The study area spans from 100°50′00″E, 27°00′00″N to 110°00′00″E, 36°00′00″N (Figure 1), in the northeast of Hengduan Mountains, a chain of mountains bordering Qinghai-Tibetan Plateau and Myanmar to its west and Sichuan Basin and Yunnan-Guizhou Plateau to its east. The annual mean temperature ranged from about 5 to 13°C in the start of the 2000s and the annual precipitation ranges from 500 to 1000 mm (Li et al., 2010). More than half of the area is covered by coniferous forests, meadows and steppes, and there are dramatic variations in topography with elevation differences over 6,000 m (Ye et al., 2015). Within the study area, the GPNP is located on the mountain edge northwest to the Sichuan Basin, covering an area of 27,134 km² (Huang et al., 2020).

We selected a set of terrestrial (and freshwater) vertebrates which fulfilled the following criteria: (1) endemic to China; (2) the conservation status was critically endangered (CR), endangered (EN) or vulnerable (VU) according to the IUCN Red List; (3) the study area covered at least 30% of the species’ current range or the species’ current range covered at least 30% of the study area; (4) occurrence data were either available, or could be geographically referenced from published maps/sightings, or the IUCN published range was available. A total of 40 species fulfilled all four criteria according to the IUCN Red List (IUCN, 2021b) and literature review (Figure 2 and Table 1) and were thus included in our study for further assessment.

The occurrence data were coordinates either directly obtained or georeferenced from published maps, or sampled (following complete spatial randomness, sampling 1,200 points per species to provide sufficient information about its presence) from published polygons of species ranges (Table 1). The Global Biodiversity Information Facility (GBIF¹) provided coordinates for six bird species and five amphibian species, using all available data due to scarcity of records. Coordinates were also found for one bird species from literature (Lu et al., 2012). Georeferenced locations were used for one mammal (Li et al., 2018, 2021). Furthermore, distribution ranges published by the IUCN were downloaded for all 40 species. As a result, there were in total 52 occurrence datasets for 40 species.

The sample size of occurrence data was set to be at least 11, as this was found to be the minimum needed to generate accurate distributions for species of low prevalence (i.e., 5% of raster cells where the species is present) (van Proosdij et al., 2016). The directly obtained and georeferenced coordinates for studied species were all below 5% prevalence in the study area, which was as expected considering they were both endemic and threatened.

Previous studies revealed important environmental variables for 23 of the 40 species and distribution range predictions for 24 species which were, when possible, included in this study (Supplementary Appendix I). For the other species, no earlier study was found. In spite of the fact that for species with a narrow niche, microhabitat conditions could be of vital importance [i.e., water velocity for Scutiger liupanensis (Zuo et al., 2017)], this study focused on generic bioclimatic, land-use, topographic and lithological variables due to constraints on data accessibility. In addition, including microhabitat variables would require a finer scale which was not possibly due to unavailability of fine scale bioclimatic and land-use variables. 19 bioclimatic variables were obtained from WorldClim version 2², 14 land-use variables were obtained from Land-Use Harmonization² (LUH2³), three topographic factors were calculated based on a digital elevation model downloaded from DIVA-GIS⁴, namely slope, aspect, and normalized topographic position index (Normalized TPI). Normalized topographic index measures the difference between elevation at the central point and the average elevation within a predetermined radius surrounding it (de Reu et al., 2013). Normalized TPI measures the topographic position as a fraction of local relief normalized to the standard deviation of the elevation (de Reu et al., 2013). Another three topographic factors were computed from inland water (rivers and lakes) and road maps obtained from DIVA-GIS and converted to distance maps, and lithology was sourced from literature (Hartmann and Moosdorf, 2012a; Table 2). Lithology describes the geochemical, mineralogical, and physical properties of rocks (Hartmann and Moosdorf, 2012b). It is a key factor in fields like landscape evolution, river chemical composition, matter supply to ecosystems, etc. (Hartmann and Moosdorf, 2012b) and was therefore included as a potential predictor variable in the models. The lithological variable contained 13 rock types in the study area ranked by ascending resistance to weathering and erosion (Dürr et al., 2005). It was used as a continuous variable representing the gradient of rocks’ sensitivity to weathering and erosion.

To identify the current (near historic) distribution ranges, we used bioclimatic and land-use variables from the period of 1970–2000, both available at a resolution of 30 arcseconds (approximately 0.8 km by 0.8 km in the study area). The future period under consideration was 2081–2100 for bioclimatic and land-use predictor variables, at a resolution of 2.5 arcminute (approximately 4 km by 4 km in the study area) and 30 arcseconds respectively. Unfortunately, 30 arcsecond data were not available for future bioclimatic variables. We therefore downscaled the data from 2.5 arcminutes to 30 arcseconds to match the resolution of the other variables. Because the land-use projections were available for each year in each period, we averaged the grid-cell values across years to match the averages as provided for the same periods for the bioclimatic variables. Topographic (available at a resolution of 30 arcseconds) and lithological variables (available at a resolution of 30 arcminutes) were kept constant between current and future periods.

Three scenarios of bioclimatic variables were used from the general circulation model The Beijing Climate Center Climate System Model (BCC-CSM2-MR), since its focus is on East Asia (Wu et al., 2019), namely the Shared Socioeconomic Pathways (SSP) 2–4.5, SSP3–7.0 and SSP5–8.5 (IPCC, 2021). Under SSP2–4.5, the most optimistic scenario in this study, CO² emission will be around the current level until 2050 after which it declines to net zero by 2100. Compared to 1850–1900, the global annual mean temperature is estimated to be between 2.1 to 3.5°C higher in 2081–2100 (IPCC, 2021). The average annual global land precipitation is projected to increase with 1.5–8% relative to 1995–2014, with stronger seasonality in more areas (IPCC, 2021). An important feature in global land-use change is the initial decrease in forest with about 43 million ha from 2000 to 2050, with a subsequent increase in forest with about 331 million ha from 2050 to 2100, resulting in a positive forest cover gain by 2100 (Hurtt et al., 2020). Under SSP3–7.0, CO² emission will double by 2100 and global annual mean temperature in the years 2081–2100 is estimated to be about 2.8 to 4.6°C higher compared to the years 1850–1900 (IPCC, 2021). The SSP3–7.0 scenario features the highest rate of net land transition globally, with a strong expansion of crop and pasture land and large-scale deforestation (Hurtt et al., 2020). Under the SSP5–8.5, the worst emission scenario in this study, CO² emission will triple by 2075 and the global annual mean temperature is estimated to be 3.3 to 5.7°C higher by the end of this century compared to 1850–1900 (IPCC, 2021). Precipitation is projected to increase with 1–13% by 2081–2100 relative to 1995–2014 with stronger seasonality in more areas (IPCC, 2021). This scenario predicts a strong expansion of cropland on the expense of pasture and forest land and large increases in irrigated area (Hurtt et al., 2020).

Data analysis was conducted in three phases (Figure 3). In the pre-processing phase, occurrence data were set to coordinate system (WSG84). For the 12 species for which we were able to obtain GBIF data or (georeferenced) occurrences from literatures, the spatial extents for generating the random pseudo-absences were the intersections of the study area and the concave shapes which enclosed all sample points, calculated in R studio (RStudio Team, 2021) using packages ‘sf’ (Pebesma, 2018) and ‘concaveman’ (Gombin et al., 2020). The spatial extents for pseudo-absence data that were extracted from IUCN ranges for all 40 species were the intersections of the study area and the occurrence ranges plus buffer zones which were calculated in the ‘rgeos’ package (Bivand and Rundel, 2021). The width of the buffer zone equaled the diameter of a circle which had the same surface area as the occurrence range. The sizes of the buffer zones were therefore species-dependent. They were designed to create adequate contrast between presence and pseudo-absence locations, i.e., areas in which species may be present and in which pseudo-absences are meaningful. This allows Maxent to better determine variable importance and relation with species presence, while avoiding AUC inflation (Barve et al., 2011; Barbet-Massin et al., 2012).

After the occurrence and random pseudo-presence data were collected, all coordinates (those obtained from the IUCN ranges as well as those from GBIF data and [georeferenced] occurrences from literatures) were cleaned by removing duplicates and records with little accuracy (i.e., fewer than two decimal places), and thinned to a 30 arcsecond resolution grid. The remaining number of occurrences per species is the sample size. Afterward, the pseudo-absence points were sampled within the above-described spatial extents and in the same sample size as the occurrences for each species.

Environmental variables were first set to the same temporal (Current = 1970–2000, future = 2081–2100) and, species specific spatial extents (see above), and then standardized to zero mean and unit variance. Topographical predictor variables were calculated in R studio (RStudio Team, 2021) using packages ‘raster’ (Hijmans, 2022), ‘spatialEco’ (Evans, 2021), and ‘sf’ (Pebesma, 2018). To shorten computation time, distances to road, river and lake were calculated in 2.5 arcminute resolution, which was aggregated and then rescaled to 30 arcsecond via bilinear interpolation after the calculation. Afterward, to reduce the multicollinearity among predictor variables, a number of environmental variables were selected per species within its spatial extent. Variable selection was based on the variance inflation factor (VIF) they expressed. The predictor variable that had the highest VIF was chosen from a cluster of correlated predictor variables (correlation coefficient > | 0.5|), as variable representing the group of correlated variables. The remaining level of multicollinearity between kept predictors was assess using VIF as calculated with the ‘usdm’ package (Naimi et al., 2014) in R studio (RStudio Team, 2021). The choice of correlation coefficient | 0.5| as threshold for predictor variable clustering was motivated by otherwise serious multicollinearity issues which would be encountered in the processing phase and the resulted indeterministic models.

In the processing phase, the maximum entropy algorithm (Maxent) (Phillips et al., 2006) was fitted to relate species occurrences to a selection of environmental variables, where we used the implementation of the ‘sdm’ package (Naimi and Araújo, 2016) in R studio (RStudio Team, 2021). Compared to other species distribution modeling algorithms, Maxent has shown good performance across sample sizes including small data sets (Wisz et al., 2008). Ten-fold cross-validation was used for occurrence datasets with more than 100 occurrence points (Naimi et al., 2014). For samples that were smaller than 100 points (in practice less than 77 points), leave-one-out cross-validation was used since ten-fold cross validation is less appropriate for datasets with limited occurrences (Pearson et al., 2007). Using leave-one-out cross-validation maximizes the use of limited data for training and equalizes the contribution of each point in testing.

For model fitting we used the ‘sdm’ function from the ‘sdm’ package (Naimi and Araújo, 2016) in R studio (RStudio Team, 2021). The model fitting started with all predictor variables. Then, variable importance and exclusion was evaluated based on Pearson correlation between the fitted values and the predictions when the focal variable was randomly permutated 5 times using the ‘getVarImp’ function from the ‘sdm’ package in R (Thuiller et al., 2009; Naimi and Araújo, 2016). The lower the correlation, and thus the higher “1-correlation,” the more important the variable was to the model (Thuiller et al., 2009; Naimi and Araújo, 2016). Changing the values of an independent variable and assessing the sensitivity of the predictions using a correlation based variable importance metric may highlight more transferable variables, which are important for the purpose of distribution projection under future environmental conditions (Phillips, 2017). Using a large range of variables that may show high collinearity in models requires a priori variable selection (Braunisch et al., 2013; Muscatello et al., 2021). We therefore implemented the model fitting process as an iteration which removed the least relevant variable(s) from the fitted model and then repeated model fitting until only two predictor variables were left in the model following methods proposed by Zeng et al. (2016). This process therefore yielded a series of models with a different number of explanatory variables, of which the first model always included all predictor variables and the last model only included two predictor variables. The models were subsequently assessed in the post-processing phases to determine the best model for each species. The model fitting and modification process was repeated three times with different random realizations to permutate variable importance for each occurrence dataset. We repeated the model fitting and modification process to assess the robustness of the results, assuming robustness if the three series of models that were generated produced similar results. As a result, for each occurrence dataset, the processing phase produced three sets of fitted models.

Model fitting was proceeded by model selection, prediction of current and future species distributions, and computation of the fraction of species ranges that were covered by the GPNP. Model selection was based on two criteria. First, the three fitted models from three model modification runs had to consist of the same explanatory variables. Secondly, the AUC score evaluated on the test datasets and the number of predictors in a model had to fulfill the requirement of model parsimony. In practice, models with AUC value +0.2 per additional predictor were selected, i.e., if, for a particular species, model 1 had 4 predictor variables and an AUC of 0.80 and model 2 had 5 predictor variables and an AUC of 0.82, we would select model 2 to be the better model. If for a species all models had AUC values <0.7, the projections made by this model were considered with low confidence. Because predictor variable elimination in the pre-processing phase was likely to capture the major variances in the environmental conditions but not necessarily the required niche by the species, models for species with small occurrence record datasets (<100) were indeterministic where each model modification process generated models with different predictors. For this reason, model selection was based on the consistency of models besides AUC scores and the number of predictors.

Finally, after the best model was selected for each species (given the above-listed criteria), probability maps of species presence were generated for the present time and for the future under three climate and land-use scenarios. For each species we thus had three probability maps for the present time and nine probability maps for the future (three times three scenarios). These maps were converted to binomial habitat suitability maps with a threshold value specified per model [i.e., Maximized sum of sensitivity and specificity (Nenzén and Araújo, 2011; Liu et al., 2016)]. We therefore generated three presence-absence maps from three simulation runs per species. Then, these three binary maps were overlayed and the intersection became the presence in a combined presence-absence map for each species. The approach we took was that if a grid-cell had a value of 1 in all three probability maps, that particular grid-cell would receive a 1 in the final suitability map. If on the other hand it only had a 1 in two out of three probability maps, it would receive a 0 in the final suitability map. A Pairwise Wilcoxon Rank Sum Test in R studio (RStudio Team, 2021) was used to assess if species range sizes changed significantly. Afterward, 40 binomial maps were summed up to produce a set of congruent suitability maps for the present and under each future scenario, where the value of a grid cell can range between 0 to 40, indicating the number of species which may occur in each cell. The current and future maps were compared to assess species distribution range changes and changes in endemism hotspots. Finally, the congruent suitability maps were overlayed with the GPNP to assess conservation coverage for now and for the future. QGIS (QGIS.org, 2022) and R studio (RStudio Team, 2021) were used for analyses.

The remaining environmental variables kept for prediction all had a VIF < 2, thus the kept predictors had very low levels of multicollinearity. Out of the 40 selected models, 31 models generated species presence-absence maps with a high confidence level (average AUC ≥ 0.7 and comparable projections across three simulation runs). The other 9 models generated species distribution projections with a low confidence level either because of low average AUC scores (<0.7) and/or because of less similar projections from three simulation runs (see Supplementary Appendix II for model per species and model statistics).

A range of variables were important in explaining the distribution ranges of the species assessed and the most important explanatory variable was generally species specific (see Supplementary Appendix II). Regarding all species together, the explanatory variable that most frequently explaining their geographic distribution was a land-use variable: the mean age of secondary land. It was one of the most important explanatory variables for 8 out of 40 species. The mean age of secondary land is predicted to decrease in most areas in and around the Sichuan Basin except for in the north where it is predicted to increase. The precipitation of the driest quarter was one of the most important explanatory variables for 7 out of 40 species. The precipitation of the driest quarter is predicted to increase in the southeast of the study area in future and decrease in the rest of the region. The species distributions mostly responded to both variables unimodally with a species-specific threshold. Furthermore, distance to the nearest lake (negative relationship with area suitability) and lithology also both appeared as in the most important explanatory variables for 6 species.

For most amphibians (6 out of 24), the mean age of secondary land was one of the most important explanatory predictor variables in which the response was usually unimodally and varied per species. For most birds (3 out of 13) the precipitation of the driest quarter was one of the most important predictor variables in which suitability increased with an increasing amount of precipitation during the driest quarter of the year until a species specific optimum after which suitability declined again. Mammals’ distributions were predicted by different environmental variables per species. Lithology was an important predictor for the distribution ranges of both reptiles; the suitability for both species to occur was highest in areas with rock types that were medium resistant to weathering and erosion. Precipitation seasonality was an important predictor variable for the Sichuan rat snake (Euprepiophis perlacea), in which probability of its occurrence first increased with increasing precipitation seasonality until a peak at 0.75 (coefficient of variation) after which it declined. Distance to the nearest lake was important for the Wa Shan Keelback (Hebius metusium) in which the probability of its occurrence declined with an increasing distance to the nearest lake.

At present, 37 out of 40 threatened endemic vertebrates in the region are present in the GPNP. However, climate and land use change is predicted to significantly alter the range size of the species assessed in future (Pairwise Wilcoxon rank sum test p < 0.001 for all scenarios [SSP2–4.5, SSP3–7.0, and SSP5–8.5]). More than half of the evaluated species (23 out of 40 species, 57.5%) is projected to lose between 80 to 100% of their current range under all three future scenarios. 78.0% of evaluated birds, more than half of the evaluated amphibians, and half of the evaluated mammals and reptiles are predicted to lose more than 80% of their suitable ranges (Supplementary Appendix III). Four species (10.0%) are predicted to lose between 15 and 80% of their ranges (Figure 4 and Supplementary Appendix III). Depending on the climate change scenario either only 17 or only 18 of these species can still occur in the GPNP. Three species are expected to maintain their current ranges (7.5%, i.e., range changes between −15 and 15%), and six species (15.0%, all amphibians) are projected to expand their distribution ranges (i.e., range changes >15%).

The impact of climate change on the extent of range size change of the species assessed is generally consistent across the three different future scenarios (Pairwise Wilcoxon rank sum test p > 0.950).

Since it is predicted that most species will contract their ranges in future, the hotspots of threatened endemic species in the study area, and in the GPNP, will also contract. Currently, the hotspots are located in the mountain ranges to the northwest and the southwest of the Sichuan basin, where a maximum of 18 evaluated species can co-occur (Figure 5A). By the end of this century, no region in the study area would possess the climatic and land-use conditions to support more than 10 threatened endemic vertebrates to occur simultaneously (Figures 5B–D). Regarding the GPNP, about 70% of the area is currently suitable to simultaneously host five to ten of the species included in our study (Figure 6A). In 2081–2100 this number is predicted to have been reduced to one to three species regardless the future scenario (Figures 6B–D).