The thinned occurrence data (Table S2) and retained variables (Table S3) for each species were used to conduct ENM analysis usingBiomod2 (Thuiller et al., 2009) based on 10 models (Table S4) from the LIG (Last Interglacial, 116–129 ka) (Shackleton et al., 2020), LGM (Last Glacial Maximum) (Clark et al., 2009), current and future (2090). For future scenarios, we used two Shared Socio–economic Pathways (SSPs): SSP126 (with 2.6 W/m² by the year 2100) and SSP585 (with an additional radiative forcing of 8.5 W/m² by the year 2100), to explore the potential distribution areas in the future (2090). The modeling steps were as follows: the bootstrap resampling method was used; 75% occurrence data for calibrating the models and the rest of the occurrence data (25%) for validating the models; the running using ten replicates; model performance was considered reliable based on two metrics: the threshold–independent area under the curve (Taberlet et al., 1998) values ≥ 0.9 and the threshold–dependent true skills statistic (TSS) values ≥ 0.7 (Table S5)