2.2 Environmental variables
We collected 29 environmental variables for this study, which can be roughly divided into four categories: bioclimate, topography, habitat, and human impact, source of all variables show in Table S2. In order to avoid model overfitting, all collected variables were processed by the band collection statistics function in ArcGIS 10.2. Remove the environment variables with high correlation (r>0.7). After screening, Bio3, Bio5, Bio9, Bio14, Bio15, slope, Veg, NDVI, POP, GDP and HI were used for this study. Bio3, Bio5, Bio9, Bio14, Bio15, slope, Veg, NDVI, POP, GDP and HI. Since variables other than future bioclimatic factors cannot be predicted and change relatively little over a short period of time, we assume that they remain constant and use their projections of the future distribution of B. gargarizans(Yang et al., 2020).
In order to reduce the uncertainty of model prediction, we chose four internationally recognized GCMs (BCC-CSM1-1, HadGEM2-ES, IPSL-CM5A-LR, MIROC5) and four RCPs (RCP2.6, 4.5, 6.0 and 8.0) published in the IPCC fifth Assessment report for prediction in 2050 and 2070.
2. 3 Model optimization and parameter setting
To establish the most stable and reliable model, we choose model calibration, which includes regularization multiplier (0.5 to 6, interval 0.5) and combination of basic feature classes and from one different sets of all layers. The best parameters were selected based on statistical significance (partial ROC), predictive ability (miss rate E=5%) and complexity level (AICc) by useing R package kuenm (Cobos, Peterson, Barve, & Osorio-Olvera, 2019). After the model parameters were determined, chose the data output format of mathematical logic to ensure that the output data was between 0 and 1. And the jackknife method was selected to evaluate the contribution of environmental variables. Then the distribution data ofB.gargarizans were randomly divided into two groups: 25% were randomly selected as the test set, and the remaining 75% were used as the training set, repeat Bootstrap replicates 10 times in MaxEnt.3.4.1. Other model parameters are selected by system default. Finally, the predictive performance of the model was validated by using the area under the receiver-operating characteristic (ROC) curve (AUC) and the mean omission error. The higher the AUC value, the more accurate the model performance. Normally, AUC > 0.9 represents excellent prediction performance of the model. The ensemble threshold of model was calculated according to the Maximizing Sensitivity and Specificity (MSS) by using the dismo package in R. The MSS method is commonly used in presence only kind of occurrence data(Hijmans, Phillips, Leathwick, & Elith, 2017; Liu, White, Newell, & Pearson, 2013).
2.4Distribution barycenter migration under different scenarios of future global climate change
We divided the study area into small grids of 0.1°×0.1° to evaluate the distribution barycenter migration of B. gargarizans in China. Assuming that our research area is composed of N small grids, the proportion of the kth grid is GK=PK×SK, and PK represents the probability of the species appearing in the kth small grid, and SKrepresents the area of the kth small grid. The coordinates of each grid are obtained by Arcgis10.2. The following formula is used to obtain the coordinates of the barycenter:
X= \(\frac{\sum_{k=1}^{k}G_{k}X_{k}}{\sum_{k=1}^{k}G_{k}}\); Y=\(\frac{\sum_{k=1}^{k}G_{k}Y_{k}}{\sum_{k=1}^{k}G_{k}}\)
Here, X and Y represent the longitude and latitude of the kth grid, respectively(He et al., 2011).
The results of four GCMs under four RCPS were averaged to obtain Pf, using △P(△P=Pf-Pc, Pc represents the current distribution probability ofB. gargarizans ) to evaluate the change trend of B. gargarizans distribution in China in 2050 and 2070: -0.1<ΔP≤0.1 indicated that the habitat suitability changes were not obvious. 0.1<ΔP≤0.3 , 0.3<ΔP≤0.5 and 0.5<ΔP≤0.7 represented slight, moderate and severe increasein in habitat suitability, respectively. -0.3<ΔP≤-0.1, -0.5<ΔP≤-0.3 and -0.7 <ΔP≤ -0.5 represented slight, moderate and severe declines in habitat suitability, respectively.