Modeling the Effect of Climate Change on the Potential Distribution of Qinghai Spruce ( Picea crassifolia Kom . ) in Qilian Mountains

Qinghai spruce forests play a key role in water conservation in the dry region of northwest China. So, it is necessary to understand the impacts of climate change on the species to implement adaptation strategies. Based on the four-emission scenario (i.e., RCP2.6 (Representative Concentration Pathway), RCP4.5, RCP6.0 and RCP8.5) set by the Intergovernmental Panel on Climate Change (IPCC) fifth assessment report, in the study, we predicted the potential distribution of Qinghai spruce (Picea crassifolia Kom.) under current and future scenarios using a maximum entropy (Maxent) model. Seven variables, selected from 22 variables according to correlation analysis combining with their contribution rates to the distribution, are used to simulate the potential distribution of the species under current and future scenarios. Simulated results are validated by area under the operating characteristic curve (AUC). Results demonstrate that elevation, mean temperature of wettest quarter, annual mean temperature, and mean diurnal range are more important in dominating the potential distribution of Qinghai spruce. Ratios of the suitable area to the total study area are 34.3% in current climate condition, 34% in RCP2.6, 33.9% in RCP4.5, 33.8% in RCP6.0, and 30.5% in RCP8.5, respectively. The warmer the climate condition is, the more area of higher suitable classification is changed to that of lower suitable classification. The ratios of real distribution area in simulated unsuitable class to the real distribution area change from 4.3% (60.7 km2) in the current climate to 13% (185 km2) in RCP8.5, suggesting that the real distribution area may decrease in the future. We conclude that there is a negative effect of climate change on the distribution of Qinghai spruce forest. The result can help decision-makers to draft adaptation countermeasures based on climate change.


Introduction
Global warming is an indisputable fact, and the global annual temperature has increased 0.85 • C between 1880 and 2012 and will continuously rise in the future [1].Climate warming has been widely proved to have a great influence on community composition and structure, vegetation pattern, and ecosystem functions, especially, on species' distribution [2][3][4][5][6][7][8][9].Predicting the potential species distribution is significant to take safe and effective countermeasures to reduce the ecological risk happening under climate change [10][11][12].
Species distribution model (SDM) is a useful tool to understand the effect of climatic change on species' distribution, which was developed by combining present data of a species with relevant environment variables [13][14][15][16][17][18][19].Many SDMs have been developed for predicting the potential distribution of species under climate change, such as the Genetic Algorithm for Rule-set Prediction (GARP) [20,21], Surface Range Envelope (SRE; usually called BIOCLIM) [22], Random Forests (RFs) [23], Ecological Niche Factor Analysis (ENFA) [24], and Maximum Entropy (Maxent) [25].They have been applied to predict the potential distribution of the selected species under current as well as the future scenarios in many places [18,26,27].
Among these models, Maxent has been widely used in the world due to its advantages: (1) needing only species occurrence data, that makes data easily collected; (2) using both continuous and categorical environmental data at the same time; (3) generating a continuous probabilistic output, that is apt to classify species suitability; (4) small sample size can meet its demand [28][29][30][31][32][33], that reducing laborious jobs in data collection; (5) facilitating model interpretation [10,25,34,35].Compared to other algorithms such as GARP, BIOCLIM, and DOMAIN, Maxent has been found to be consistently better in its prediction performance [25,36].It has been extensively used in the studies on the prediction of species distribution [18,26,27,37,38].
Qilian Mountains located in the northwest China and surrounded by desert and Gobi, is a crucial water source in the northwest China because of the higher precipitation.Many inland rivers, such as Shiyanghe River, Heihe River and Shulehe River, are originated from Qilian Mountains [39].Qinghai spruce (Picea crassifolia Kom.), the dominant tree species in the Qilian Mountains, is distributed in the north-facing slope.It is characterized by a small population, a limited geographic range, and habitat specialization [40].The population with these characteristics has remarkable sensitivity to climate change [41,42].Qinghai spruce forest plays an important role in ecological service functions, especially water conservation.Given this, Qinghai spruce becomes the protected object in Qilian Mountains National Natural Reserve.One widespread hypothesis is that global warming will change distribution ranges of species.Understanding changed distribution ranges of Qinghai spruce is meaningful for maintaining ecological services and protection strategies.
The objectives in this study are to: (1) predict the potential distribution of the species under the current situation; (2) forecast its suitability areas under four future climate scenarios, and (3) evaluate the effects of climate change on the distribution of Qinghai spruce.Identifying the shift ranges of suitable area under future scenarios is the innovation in the study.

Study Area
Qilian Mountains selected as our study area is situated at the northeast edge of the Tibet Plateau, between 99 • 25 -103 • 28 E and 36 • 45 -39 • 36 N (Figure 1).The elevation ranges from 2000 to 6000 m.In the southeast, the Southeast Asian monsoon is the influencing weather system, and in the northwest the westerlies bring relatively dry air from Central Asia.Two atmospheric circulation systems together with high elevation differences result in a significant precipitation gradient.Observations at various meteorological stations (for the period 1960-2009) present that horizontally, higher mean annual precipitation (525.9 mm) appeared in the southeast part and lower mean annual precipitation (293.1 mm) occurred in the northwest part, vertically, the mean annual precipitation ranged from 250 mm in lowlands to 700 mm in high elevations.Approximately 88% of precipitation occurs from May to September.The mean annual temperature is 6 °C in lowlands and −10 °C in high elevations [43].The combination of precipitation and temperature results in horizontal and vertical differentiations of vegetation and soil.Vegetation types in the study area are displayed in spatial sequences (from low to high elevations): desert, desert steppe, steppe, steppe-forest, sub-alpine shrubby meadow and alpine frost-action barren zone.Qinghai spruce forests distributed from 2500 m to 3500 m are in the steppe-forest zone [44].
steppe-forest, sub-alpine shrubby meadow and alpine frost-action barren zone.Qinghai spruce forests distributed from 2500 m to 3500 m are in the steppe-forest zone [44].

Qinghai Spruce Distribution Data
The date on the locations of Qinghai spruce was obtained from field survey conducted from 2015 to 2017 in the study area.A total of 178 species occurrence were recorded.To minimize the spatial autocorrelation of the species occurrence points and reduce its impact on the model's results, the sampling sites were selected based on different environmental conditions and the distance between two sampling points was more than 5km.The real distribution map of forests was collected from Gansu Qilian Mountains National Nature Reserve administration.We extracted the distribution of Qinghai spruce forests after digitalizing the forests map.

Environmental Data
Environmental variables included three terrain variables and 19 bioclimatic variables.Terrain variables like elevation, aspect, slope with a 90m spatial resolution were downloaded from the USGS website (www.srtm.usgs.gov),which were re-sampled into 1 km spatial resolution using the nearest neighbor re-sampling technique in ArcGIS10.3.The 19 bioclimatic variables for the current and future scenarios with a 1 km spatial resolution were downloaded from the WorldClim data set (http://www.worldclim.org/bioclim).The dataset was generated by interpolation of observed weather data using a thin-plate smoothing spline during the period of 1950-2000.The interpolation method considers latitude, longitude and elevation as independent variables [45,46].The future climate scenarios were present by 2070s data (the average data for 2060-2080) modeled by the Community Climate System Model version 4 (CCSM4) representing four future greenhouse gases concentration trajectories (i.e., RCP2.6 (Representative Concentration Pathway), RCP4.5, RCP6.0, and RCP8.5).The four scenarios described the total radiative forcing values in 2100 would have reached 2.6 W/m 2 , 4.5 W/m 2 , 6.0 W/m 2 , and 8.5 W/m 2 over the value in preindustrial period [47].The 19 variables downloaded were on a global scale, so they were extracted for the study area by GIS tool.
In order to avoid the cross-correlation within 19 bioclimatic variables and three terrain variables, first, we assessed the contribution rate of each variable by the jackknife test of Maxent model, then we conducted the multi-collinearity test using Pearson's correlation coefficient in IBM-SPSS statistical software (Statistical Product and Service Solutions, Version 20.0, SPSS Inc., Chicago, IL, USA), If two variables have a Pearson correlation > |0.8|,only one with a high contribution rate was selected [48].

Qinghai Spruce Distribution Data
The date on the locations of Qinghai spruce was obtained from field survey conducted from 2015 to 2017 in the study area.A total of 178 species occurrence were recorded.To minimize the spatial autocorrelation of the species occurrence points and reduce its impact on the model's results, the sampling sites were selected based on different environmental conditions and the distance between two sampling points was more than 5 km.The real distribution map of forests was collected from Gansu Qilian Mountains National Nature Reserve administration.We extracted the distribution of Qinghai spruce forests after digitalizing the forests map.

Environmental Data
Environmental variables included three terrain variables and 19 bioclimatic variables.Terrain variables like elevation, aspect, slope with a 90m spatial resolution were downloaded from the USGS website (www.srtm.usgs.gov),which were re-sampled into 1 km spatial resolution using the nearest neighbor re-sampling technique in ArcGIS10.3.The 19 bioclimatic variables for the current and future scenarios with a 1 km spatial resolution were downloaded from the WorldClim data set (http://www.worldclim.org/bioclim).The dataset was generated by interpolation of observed weather data using a thin-plate smoothing spline during the period of 1950-2000.The interpolation method considers latitude, longitude and elevation as independent variables [45,46].The future climate scenarios were present by 2070s data (the average data for 2060-2080) modeled by the Community Climate System Model version 4 (CCSM4) representing four future greenhouse gases concentration trajectories (i.e., RCP2.6 (Representative Concentration Pathway), RCP4.5, RCP6.0, and RCP8.5).The four scenarios described the total radiative forcing values in 2100 would have reached 2.6 W/m 2 , 4.5 W/m 2 , 6.0 W/m 2 , and 8.5 W/m 2 over the value in preindustrial period [47].The 19 variables downloaded were on a global scale, so they were extracted for the study area by GIS tool.
In order to avoid the cross-correlation within 19 bioclimatic variables and three terrain variables, first, we assessed the contribution rate of each variable by the jackknife test of Maxent model, then we conducted the multi-collinearity test using Pearson's correlation coefficient in IBM-SPSS statistical software (Statistical Product and Service Solutions, Version 20.0, SPSS Inc., Chicago, IL, USA), If two variables have a Pearson correlation > |0.8|,only one with a high contribution rate was selected [48].

Ecological Niche Modeling
We used the Maxent model to predict the current and future potential distribution of Qinghai spruce in the Qilian Mountains.The Maxent model is a machine learning technique that has successfully been applied to model species distribution by only occurrence data and environmental data [25].The Maxent model contains two accessional procedures; one is area under the curve (AUC) calculation.AUC values were used to evaluate the model performance, which are automatically generated using random pseudo absence background points [49], ranging from 0 to 1. Performance of the model can be divided into five groups according to AUC values: excellent (>0.9), good (0.8-0.9), accepted (0.7-0.8), bad (0.6-0.7), and invalid (<0.6) [26].The other is the jackknife test; it can be used to evaluate the importance of variables.
The species distribution point data should be saved in CSV format and contains coordinate information such as species name, longitude and latitude.All environmental variables converted to ASCII raster grids and species occurrence coordinates were converted to decimal degrees in ArcGIS 10.3.Other settings were kept as default (500 iterations, 0.00001 convergence threshold, 10,000 maximum background points).The program ran with 75% of presence locations and tested the model performance with the remaining 25% of the presence locations [25].A detailed mathematical definition of Maxent was described in Phillips et al. [25].
The output of the model is a continuous probabilistic layer, ranging from 0 to 1. Areas with higher values imply more favorable conditions for the species growth [25].We selected the minimum of the output by training presences (MTP) as a threshold or "cutoff" value for each scenario [50].Ecologically, the MTP can be interpreted to contain those cells that are predicted to be at least as suitable as those where the species was identified as present.We divided habitat suitability into four classes: unsuitability, low suitability, moderate suitability, and high suitability.

Variables Analysis
According to the contribution rate and Pearson's correlation coefficient (see Table S1), only seven variables were retained for modelling (see Table 1 in bold), including Elevation, Slope, Aspect, Bio1, Bio2, Bio3, and Bio8.Bio1 is annual mean temperature, Bio2 is mean diurnal range, that is, mean of monthly (max temp-min temp), Bio3 is isothermality that equals (Bio2/Bio7) (* 100), and Bio7 is temperature annual range that is the maximum temperature of warmest month minus minimum temperature of coldest month, bio8 is mean temperature of wettest quarter.The wettest quarter is three months with the most rainfall during a year.The driest quarter is three months with minimum precipitation during a year.In the study area, the wettest quarter includes June, July, and August, and the driest quarter is December, January, and February.The wettest quarter is consistent with the warmest season.During the wettest quarter, precipitation is in the form of rain.The variables in bold were key variables selected by their contribution rates and multi-collinearity test.RCP: (Representative Concentration Pathway).
We suppose the terrain variables would not change in the future scenarios.The contribution of the other four variables to the climate variability under the four scenarios was shown in Table 2. Bio2 and Bio3 had little contribution to the climate variability, while Bio1 had great contribution to the climate variability.Bio1 increased rapidly with the increase in the greenhouse gas concentrations, being 0.30

Accuracy of the Maxent Model
AUC values were greater than 0.98 under the current and four climatic scenarios (Table 3).According to the performance classification standard, the prediction accuracy was excellent in our study, as AUC values are more than 0.9 in both training and test.
The jackknife test indicated that Elevation, Bio8, Bio1 and Bio2 provided very high gains when they were used independently (Figure 2).According to gains being more than 1, Figure 2A showed that potential distribution of Qinghai spruce under the current was significantly affected by Bio8 with gain 1.57, elevation with gain 1.53, and Bio1 with gain 1.43.The potential distribution of Qinghai spruce in the RCP2.6 scenario was mainly contributed by Bio8, elevation, Bio1, and Bio2, their gains were 1.66, 1.55, 1.43, and 1.04, respectively (Figure 2B).The important variables under the RCP4.5 scenario were Bio8 (gain, 1.68), elevation (gain, 1.55), Bio1 (gain, 1.44), and Bio2 (gain, 1.08) (Figure 2C).The variables most constraining the potential distribution of Qinghai spruce under the RCP6.0 scenario were the same as that under RCP4.5 with the same gains except for Bio8 (gain 1.65) (Figure 2D).The potential distribution of the species under the RCP8.5 scenario was most strongly associated with Bio8, elevation, Bio1, and Bio2, their gains were 1.67, 1.55, 1.42, and 1.01, respectively (Figure 2E).On the whole, Bio8, elevation, Bio1, and Bio2 were more important variables for predicting the distribution of Qinghai spruce.

Potential Distribution of Qinghai Spruce under Different Climate Scenario
The suitability of Qinghai spruce was classified using four probability classes based on Khafaga [51] and Remya [52].The MTP threshold, or the level at which no omission errors, was detected as 0.089.That is, areas with probability <0.089 were unsuitability areas.The range between 0.089 and 1 was divided into three classes, 0.089-0.3named as low suitability, 0.3-0.6 as moderate suitability, >0.6 as high suitability.So, four classifications of suitability were obtained.
The potential distribution maps of Qinghai spruce under current climate and four scenarios were shown in Figure 4.The Qinghai spruce is distributed mainly over the middle and eastern parts of the study area.The simulated distribution area of Qinghai spruce under current climate and four scenarios are considerably larger than real distribution area (Figure 4).
As shown in Table 4, the distribution area with moderate suitability under current climate is 4250.7 km 2 , the area with low suitability is 2994.1 km 2 , and the area with high suitability is 1312.0km 2 .The total suitable area occupied 34.3% of the study area.In four scenarios, the total suitable area gradually decreased with climate warming, being 8556.8,8473.1, 8438.5, 8420.6, and 7604.2 km 2 under current, RCP2.6,RCP4.5, RCP6.0 and RCP8.5, respectively (Table 4).The warmer the climate

Potential Distribution of Qinghai Spruce under Different Climate Scenario
The suitability of Qinghai spruce was classified using four probability classes based on Khafaga [51] and Remya [52].The MTP threshold, or the level at which no omission errors, was detected as 0.089.That is, areas with probability <0.089 were unsuitability areas.The range between 0.089 and 1 was divided into three classes, 0.089-0.3named as low suitability, 0.3-0.6 as moderate suitability, >0.6 as high suitability.So, four classifications of suitability were obtained.
The potential maps of Qinghai spruce under current climate and four scenarios were shown in Figure 4.The Qinghai spruce is distributed mainly over the middle and eastern parts of the study area.The simulated distribution area of Qinghai spruce under current climate and four scenarios are considerably larger than real distribution area (Figure 4).

The Transfer Matrix of Suitable Distribution Areas under the Current and Future Scenarios
The transfer matrix of suitable distribution classes under the current and future scenarios showed that the warmer the climate condition is, the more area with higher-level suitability is changed to that with lower-level suitability (Figure 5).For example, there were 443.5 km 2 areas with As shown in Table 4, the distribution area with moderate suitability under current climate is 4250.7 km 2 , the area with low suitability is 2994.1 km 2 , and the area with high suitability is 1312.0km 2 .The total suitable area occupied 34.3% of the study area.In four scenarios, the total suitable area gradually decreased with climate warming, being 8556.8,8473.1, 8438.5, 8420.6, and 7604.2 km 2 under current, RCP2.6,RCP4.5, RCP6.0 and RCP8.5, respectively (Table 4).The warmer the climate is, the less suitable the area is.The suitable distribution area occupied 34%, 33.9%, 33.8%, 30.5% of study area under four scenarios (i.e., RCP2.6,RCP4.5, RCP6.0 and RCP8.5), respectively.

The Matrix of Suitable Distribution Areas under the Current and Future Scenarios
The transfer matrix of suitable distribution classes under the current and future scenarios showed that the warmer the climate condition is, the more area with higher-level suitability is changed to that with lower-level suitability (Figure 5).For example, there were 443.5 km 2 areas with low suitability, 58.8 km 2 areas with moderate suitability and 3.9 km 2 areas with high suitability area in current climate condition would be changed to unsuitable area in RCP2.6.There was 598.3 km 2 of area with moderate suitability and 18.0 km 2 of area with high suitability area in current climate condition that would be changed to area with low suitability in RCP2.6.There was 367.6 km 2 of area with high suitability in current climate condition that would be changed to area with moderate suitability in RCP2.6.Totally, 1490.1 km 2 of area with higher level suitability in current climate condition would be changed to area with lower level suitability in RCP2.6.There was 1108.4 km 2 of area with low suitability, 684.8 km 2 of area with moderate suitability and 149.6 km 2 of area with high suitability area in current climate condition would be changed to unsuitable area in RCP8.5.There was 1192.1 km 2 of area with moderate suitability and 191.1 km 2 of area with high suitability in current climate condition would be changed to area with low suitability in RCP8.5.There was 538.0 km 2 of area with high suitability in current climate condition would be changed to area with moderate suitability in RCP8.5.Totally, 3672.9 km 2 of area with higher level suitability in current climate condition would be changed to area with lower level suitability in RCP8.5.
Forests 2018, 9, x FOR PEER REVIEW 10 of 16 low suitability, 58.8 km 2 areas with moderate suitability and 3.9 km 2 areas with high suitability area in current climate condition would be changed to unsuitable area in RCP2.6.There was 598.3 km 2 of area with moderate suitability and 18.0 km 2 of area with high suitability area in current climate condition that would be changed to area with low suitability in RCP2.6.There was 367.6 km 2 of area with high suitability in current climate condition that would be changed to area with moderate suitability in RCP2.6.Totally, 1490.1 km 2 of area with higher level suitability in current climate condition would be changed to area with lower level suitability in RCP2.6.There was 1108.4 km 2 of area with low suitability, 684.8 km 2 of area with moderate suitability and 149.6 km 2 of area with high suitability area in current climate condition would be changed to unsuitable area in RCP8.5.There was 1192.1 km 2 of area with moderate suitability and 191.1 km 2 of area with high suitability in current climate condition would be changed to area with low suitability in RCP8.5.There was 538.0 km 2 of area with high suitability in current climate condition would be changed to area with moderate suitability in RCP8.5.Totally, 3672.9 km 2 of area with higher level suitability in current climate condition would be changed to area with lower level suitability in RCP8.5.

The Relationship Between Qinghai Spruce' Actual Distribution and Potential Distribution Under Current and Four Scenarios
According to the actual distribution area (Figure 4), there was 1422.1 km 2 of Qinghai spruce forests in the study area, mainly distributed in the middle and east parts of the reserve.Overlapping the real distribution with potential distributions of current and future scenarios, we found that the area of real distribution in unsuitable class increased from 60.7 km 2 (4.3%) in the current climate to

The Relationship Between Spruce' Actual Distribution and Potential Distribution Under Current and Four Scenarios
According to the actual distribution area (Figure 4), there was 1422.1 km 2 of Qinghai spruce forests in the study area, mainly distributed in the middle and east parts of the reserve.Overlapping the real distribution with potential distributions of current and future scenarios, we found that the area of real distribution in unsuitable class increased from 60.7 km 2 (4.3%) in the current climate to 185 km 2 (13%) in RCP 8.5 (Table 5).

Discussion
It is necessary to understand the interrelationship between the size of species geographic range and species extinction risk under global climate change scenarios [53,54].Which is crucial to set up management strategies for the habitat conservation and sustainability of species in the future [55,56].For the purpose, we must identify the distribution of species under the current and climate change scenarios.SDM is one of the important tools for determining species distribution [57].Many SDMs have been developed, among them the Maxent model performs well with its advantages mentioned in introduction to predict species' distribution ranges at macro-scales [58,59].In this study, we used the Maxent model to predict Qinghai spruce distribution ranges on the current and future scenarios and presented an assessment about habitat suitability of the species, an extremely important ecological and hydrological tree species in the Qilian Mountains.
Model accuracy assessment is an essential step to ensure that the species distribution models (SDMs) reflect species-habitat relationships [60].Generally, AUC is selected as an assessment standard.Models with a value larger than 0.7 exhibit good performance [26,61].In our study, AUC values were greater than 0.964 for Qinghai spruce under the current and four climatic change scenarios, indicating that the model has excellent performance for simulating the distribution of the species.
Our results showed that the variables related to air temperature (Bio1, Bio8) and elevation were dominant variables in determining the suitable habitat for Qinghai spruce.Previous studies demonstrated that distribution of the species was affected by temperature in July, annual precipitation and aspect [44].Xu et al [39] indicated mean temperature of the warmest quarter, precipitation of the wettest quarter, annual solar radiation, and topographic wetness index were important variables.High suitability of Qinghai spruce was the optimal combination of average annual precipitation, average air temperature in July, and solar radiation, being 380 mm, 11.5 • C, and 2 × 10 3 kWh/m 2 , respectively [62].Therefore, the temperature factor was the dominant variable that determined the distribution of Qinghai spruce.Qinghai spruce was mainly distributed on the shady slope.Increased temperatures may lead to a decrease in the distribution of the species in low altitude and an increase in the growing season.Under different climate change scenarios, the increase of the temperature in the study area was much greater than that in precipitation, which may lead to drought in the study area.So, it is not conducive to the growth of spruce.It should be noted that Qilian Mountains regions are dominated by westerly jet, that is, so-called "silk-road" pattern [63].Song et al. [64] and Song and Zhou [65] revealed the pattern influencing the air temperature.Therefore, the "silk-road" pattern should be further considered to study the effect of temperature on the distribution of Qinghai spruce.
Our research found that suitable area would decrease under four climate change scenarios (RCP2.6,RCP4.5, RCP6.0 and RCP8.5).This is consistent with results of previous studies, which suggested that the global warming has negative effects on species distribution ranges [7,66].The suitable area decreased under RCP8.5 is larger than the other three scenarios due to different changes of temperature under four scenarios.RCP2.6 is the lowest carbon emission scenario, RCP4.5 and RCP6.0 are medium carbon emission scenarios, whereas, RCP8.5 is an extreme carbon emission scenario.The increase in temperature was much greater under RCP8.5 than under the other three scenarios (Table 2).By contrast, Xu et al. [67] reported that the potential distribution area of Qinghai spruce increased under climate warming.The difference may be attributed to different selected climate models.Previous studies demonstrated that different models and different realizations in the same model may show quite different historical temperature evolution [68] There are 19 general circulation models (GCMs) in WorldClim, in these studies, we selected the CCSM4 model because it was widely adopted to simulate future climate change projection [69,70].However, the comparative analysis of climate model performance in Qilian Mountains is needed in further research, which could provide the criterion for selecting models.In this way, we can use the average value of climatic variables obtained from GCMs to run Maxent model in the Qilian Mountains, and thus avoiding the uncertainty of future climate projections.In addition, input variables' selection is also a key step.Different results about the effect of climate warming on Qinghai spruce distribution in the Qilian Mountains may be caused by input variables.The study of Xu et al. [67] had no screening variables, and the sampling points were mostly concentrated in the central and eastern regions.The multi-collinearity problem existed among variables and the autocorrelation problem of the sampling points was not considered.This may be the reason why our results are different from theirs.
Many variables have multi-collinearity that results in an overfitting simulation of species distribution [71].Therefore, it is necessary for screening variables to run SDMs.There are various selection techniques such as expert knowledge method [27,62], pair-wise correlation analyses [72], X 2 test [73], PCA analysis [74], and jackknife analyses [75][76][77][78].In our study, we used jackknife test and Pearson's correlation coefficient to select variables.Variables were selected with the higher contribution and least correlation (r < 0.8).Comparative study of variable selection methods will be conducted further for promoting the accuracy of SDM.
Our results showed that the predicted potential distribution under current and future scenarios was larger than realized distribution.This is possibly because the Maxent model predicted the species fundamental niche rather than realized niche [27,31].In reality, a species might have failed to disperse due to geographic barriers, human disturbance or associated species' competition [31,42,79,80].In our study area, realized distribution of the species had suffered from intensive deforestation.It was reported that some area of Qinghai spruce forest has been converted into grasslands [44] in the middle part of the study area and has been transformed to farmland in the east [39].The comparative analysis of potential species distribution under current climate and the realized species distribution could identify human activities and provide information about prior areas to restore the species.

Conclusions
In this study, we successfully modeled the potential distribution of Qinghai spruce for the current and future climate change scenarios.Our results verify the hypothesis that global warming will change distribution ranges of species, which is crucial for understanding spatiotemporal dynamics of Qinghai spruce under climate change scenarios.The current suitable areas predicted should have two measures to protect the distribution area of the species.One is conservation prioritization measure; the other is restoration measure.Conservation prioritization area is the predicted suitable area which has real distribution.The restoration area is the place where the predicted suitable area which has no real distribution, especially the high suitable area in the middle and eastern part of the Qilian Mountains.Our results could be used to provide reliable information on devising adaptive responses for the sustainable management of the species.

Figure 1 .
Figure 1.Location of the study area in Gansu Province, China, and its digital elevation model (DEM).

Figure 1 .
Figure 1.Location of the study area in Gansu Province, China, and its digital elevation model (DEM).

Forests 2018, 9 , 16 Figure 2 .
Figure 2. (A) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under current climate condition; (B) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under climate change scenarios RCP2.6; (C) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under climate change scenarios RCP4.5; (D) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under climate change scenarios RCP6.0; (E) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under climate change scenarios RCP8.5.(The regularized training gain describes how much better the Maxent distribution fits the presence data compared to a uniform

Figure 2 .
Figure 2. (A) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under current climate condition; (B) The results of the jackknife test of variables' contribution in modeling Qinghai spruce's habitat distribution under climate change scenarios RCP2.6;

Figure 3 16 Figure 3 .
Figure 3 showed that response curves of seven environmental variables to the potential distribution of Qinghai spruce.From Figure 3 we could see: the highest suitability for Qinghai spruce occurred in the areas where elevation was from 2650 to 3050 m, annual mean temperature (Bio1) was between 1.5 and 2.5 • C, mean temperature of wettest quarter (Bio8) between 11.0 and 12.5 • C, Mean diurnal range (Bio2) from 11.0 to 12.1 • C, Isothermality (Bio3) from 31.8 to 33.0, aspect from 310 to 40 • , slope between 4.8 to 30 • .Forests 2018, 9, x FOR PEER REVIEW 8 of 16

Forests 2018, 9 , 16 Figure 4 .
Figure 4.The actual distribution and potential of Qinghai spruce under current climate and four scenarios.

Figure 4 .
Figure 4.The actual distribution and potential of Qinghai spruce under current climate and four scenarios.

Table 1 .
22Environmental variables downloaded and variables contribution rate (%) to the species distribution under each scenario.

Table 2 .
Climatic variables selected in this study and their percent contribution to each scenario.

Table 3 .
Results of receiver operating characteristic (ROC) analysis under current climate and four future scenarios.

Table 4 .
Habitat suitability classes area of Qinghai spruce.

Table 4 .
Habitat suitability classes area of Qinghai spruce.

Table 5 .
The statistic area of Qinghai spruce' actual distribution in four classes under current and four scenarios.