The Potential Distribution of Juniperus rigida Sieb. et Zucc. Vary Diversely in China under the Stringent and High GHG Emission Scenarios Combined Bioclimatic, Soil, and Topographic Factors

: Global warming poses an enormous threat to particular species with shifts to their suitable habitat. Juniperus rigida Sieb. et Zucc., an endemic species to East Asia and a pioneer species in the Loess plateau region, is endangered because of the shrinking and scattered habitat threatened by climate change. For the sake of analyzing the impact of climate warming on its possible habitat, we herein projected the current and future potential habitats of J. rigida in China and comparatively analyzed the ecological habitat changes in three main distribution regions. There were 110 specimen records of J. rigida collected across China and 22 environmental datasets, including bioclimatic variables and soil and topographical factors, selected by the Pearson Correlation Coefﬁcient. The MaxEnt model based on specimen presence and environmental factors was used for projecting the potential habitats of J. rigida in China in the 2050s and the 2070s of RCP 2.6 and RCP 8.5 scenarios. The results indicated an excellence model performance with the average value of the area under curve (AUC) is 0.928. The mean temperature of the driest quarter (MTDq) and the temperature annual range (TAR) provided important contributions to the potential distribution of J. rigida . There were three main distribution areas in China, the Xinjiang region, the Loess-Inner Mongolian Plateau region, and the Changbai Mountain region. The distribution increased overall in area under RCP 2.6 and decreased for RCP 8.5. The mean altitude of the core distribution shifted upward in general under both scenarios. The Loess–Inner Mongolian Plateau region is the biggest distribution, encompassing ca. 61.39 × 10 4 km 2 (86.87 × 10 4 km 2 in China). The region threatened most by climate change is located in the Changbai Mountain distribution, with the centroid of the cord suitable habitat migrating southwest about 227.47 and 260.32 km under RCP 2.6 and RCP 8.5 by the 2070s. In summary, these ﬁndings provided a well-grounded understanding of the effect of climate change on ecological distribution and furnished theory evidence for the protection, management, and sustainable use of J. rigida .


Introduction
Greenhouse gas (GHG) emissions have led to global warming, which has now become a global issue [1]. It was projected that, by the end of the 21st century (2081-2100), relative to 1986-2005 with the reason of increasing GHG emission during the industrial age, the global average surface temperature is likely to increase by 0.3-1.7 • C under the RCP 2.6 scenario, 1.1-2.6 • C under the RCP 4.5 scenario, and could be 1.4-3.1 • C under the RCP 6.0 scenario and 2.6-4.8 • C under the RCP 8.5 scenario [2]. Changing climatic conditions have had a profound impact on the distribution range of species, placing the global ecological environment under inconceivable pressure [3,4], determining the viability and conservation

Collection and Processing of Species Distribution Data
The current occurrence records of Juniperus rigida Sieb. et Zucc. throughout China derived from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/, accessed on 20 November 2020), the Chinese Virtual Herbarium databases (CVH, http:// v5.cvh.org.cn/, accessed on 28 November 2020) and Specimen Resources Sharing Platform for Education (http://mnh.scu.edu.cn/, accessed on 5 December 2020). We obtained 557 occurrence records. The duplicates and those that lacked exact geo-coordinate record points were not taken into account. Finally, there remained 110 occurrence points that could carry out the MaxEnt model program and develop the distribution map integrated the topographic information with the ArcGIS 10.2 (Esri, Redlands, CA, USA) ( Figure 1). Those records were divided into three regions, which will be emphasized.

Collection and Selection Climate, Soil, and Topographical Data
The current period  bioclimatic variables with 30 arc seconds (approximately 1 km 2 ) spatial resolution were downloaded from the WorldClim-Global Climate Database (http://worldclim.org/, accessed on 11 October 2020). They represented annual or seasonality trends and extreme or limiting environmental factors [20] and were the most crucial variables affecting potential species distribution [21].
Future period (2050s and 2070s) bioclimatic variables with the same geographical accuracy and variable names with current bioclimatic variables were extrapolated by a global circulation model. There were four representative concentration pathways (RCPs) established in the 5th report of the IPCC (Intergovernmental Panel on Climate Change) [22], which means the future climate change was diverse and uncertain. RCP 2.6 was a lower greenhouse gas emission scenario and was paid widespread attention because it could achieve the goal of no more than 2 • C of global warming by 2100 relative to pre-industrial levels [23,24]. Holding the increase in the global average temperature below 2 • C above pre-industrial levels and striving to limit the increase to 1.5 • C above pre-industrial levels is the common goal for all countries in the world by reducing greenhouse gas emissions [25]. RCP 8.5 was a robust GHG emission scenario without additional efforts to constrain emissions, which posed a high risk of abrupt and irreversible regional-scale change in the composition, structure, and function of ecosystems [22]. Therefore, we chose RCP 2.6 and RCP 8.5 under BCC-CSM1-1 climate change modeling, an efficient global climate tool for the simulation of future climatic conditions that performs perfectly for the simulation of precipitation and temperature in the Chinese region [26,27]. There were two periods for model projection; the 2050s and 2070s represent 2041-2060 and 2061-2080, respectively.

Collection and Selection Climate, Soil, and Topographical Data
The current period  bioclimatic variables with 30 arc seconds (approximately 1 km 2 ) spatial resolution were downloaded from the WorldClim-Global Climate Database (http://worldclim.org/, accessed on 11 October 2020). They represented annual or seasonality trends and extreme or limiting environmental factors [20] and were the most crucial variables affecting potential species distribution [21].
Future period (2050s and 2070s) bioclimatic variables with the same geographical accuracy and variable names with current bioclimatic variables were extrapolated by a global circulation model. There were four representative concentration pathways (RCPs) established in the 5th report of the IPCC (Intergovernmental Panel on Climate Change) [22], which means the future climate change was diverse and uncertain. RCP 2.6 was a lower greenhouse gas emission scenario and was paid widespread attention because it could achieve the goal of no more than 2 °C of global warming by 2100 relative to preindustrial levels [23,24]. Holding the increase in the global average temperature below 2 °C above pre-industrial levels and striving to limit the increase to 1.5 °C above pre-industrial levels is the common goal for all countries in the world by reducing greenhouse gas emissions [25]. RCP 8.5 was a robust GHG emission scenario without additional efforts to constrain emissions, which posed a high risk of abrupt and irreversible regional-scale Insolation duration and relative humidity were provided by the China Meteorological Data Service Center (CMDC, http://data.cma.cn/en, accessed on 25 October 2020). These data were observed by 613 data and basic ground meteorological observation stations in China, and the mean value of each observation station was calculated during 1951~2000 to develop the relevant lay using kriging interpolation with the software of ArcGIS 10.2 (Esri, Redlands, CA, USA). The China Soil Map-Based Harmonized World Soil Database (v1.1) was derived from the Cold and Dry Area Scientific Data Center (http://westdc.westgis. ac.cn, accessed on 10 September 2020). The aspect and slope variables were generated based on the altitude dataset obtained from ISDSP (http://datamirror.csdb.cn/, accessed on 13 September 2020).
All the data were processed with the same spatial resolution as the climate data. Many of the environmental variables that are spatially correlated can result in model overfitting [16]. We used the Pearson correlation coefficient to examine the multi-collinearity of variables and removed one or more highly correlated environmental factors if p > 0.8. Ultimately, we retained 22 factors for modeling (Table 1) and converted the data to the ASCII format.

Model Establishment and Accuracy Test
The maximum entropy model (MaxEnt, version 3.4.1) is available freely online (http://www.cs.princeton.edu/wschapire/Maxent/, accessed on 4 June 2020). All the variables selected were imported into the model to project the potential distribution. The current climate data were replaced with the future climate data when the potential distribution is projected in the future climate scenario, while the soil and terrain data remain unchanged. The Akaike information criterion (AICc) visually manifest the optimal combination for two parameters, feature classes and regularization multiplier, which were of great importance to the fitting degree of MaxEnt [28]. We ran the model with default values of the two parameters, and the good performance had been verified [29].
The fitted model precision was evaluated based on the area under the receiver operating characteristic curve (AUC-ROC), with the criterion of values for AUC > 0.9 representing more adequate fit [30] by 10 fold cross-validation. The receiver operating characteristic curve (ROC) graphically displayed the performance of the model [31]. The output was in cloglog format with the probability (0~1) of the presence of J. rigida in each grid of 30 arc-seconds resolution. The maximum training sensitivity plus specificity criterion we utilized had been recognized as one of the best threshold (TH) selection methods on account of the optimized specificity and sensitivity using training data [32]. The maximum number of background points and iterations were set to 10,000 and 3000, respectively, which ensured that the model had adequate time to converge, while the other parameters defaulted according to the MaxEnt Tutorial [33].

Potential Habitat Assessment
We classified the distribution into 4 levels, unsuitable (<TH), general (TH~0.5), moderate (0.5~0.8), and core (0.8~1), and illustrated them in pictures, according to the suitable indices for each climate scenario provided by MaxEnt.
The distribution changes in the future climate scenarios were in three patterns, new, disappeared, and unchanged, compared with current habitat. To illustrate the change in future scenarios, we replaced the continuous probability value with the discrete binary and recalculated the grid value using the following formula. The grids with a probability above the threshold were assigned to 1 representing suitability habitats, and the others to 0 representing unsuitability habitats.
It can be calculated using the formula following.
where GCC was the grid code under the current climate, GCF was the grid code under the future climate. When X = 0 represented the unsuitable distribution under current and future climate, X = 1 represented the new distribution under future climate, X = 2 represented the disappeared distribution under future climate, and X = 3 represented the unchanged distribution.

The Shift of the Core Distribution Centroid
We assumed the core distribution area of Juniperus rigida Sieb. et Zucc. to a single centroid point that could be calculated with the help of the SDM tool-box, a type of pythonbased GIS software. The core distribution is the geospatially highly suitable area for J. rigida. The magnitude and direction of distribution area change were analyzed more clearly by comparing the centroid point position under the current and future climate.

Model Evaluation and Important Variables
The MaxEnt model for Juniperus rigida Sieb. et Zucc. generated excellent evaluations under the given parameter setting conditions, with the average test AUC value of 0.928 ± 0.025 for the 10 times replicate runs. Making a comprehensive analysis of the percent contribution and permutation importance of all predictive variables provided by the MaxEnt result, those variables, slope (Slo), altitude (Alt), mean temperature of driest quarter (MTDq), and temperature annual range (TAR), were the domain environmental factors with higher percent of contribution relative to other variables. Furthermore, annual mean temperature (AMT), annual precipitation (AP), temperature seasonality (TS), and precipitation of wettest month (PWm) made a greater contribution to the permutation importance. The cumulative contributions of these factors hit a value as high as 76.8% (Table 1).

The Potential Distribution of J. rigida under the Current Climate
The suitability was classified into four levels and was shown in Figure 2A. The simulated suitable habitat under the current climate scenario was mainly distributed in three regions, the Xinjiang region, the Loess-Inner Mongolian Plateau region, and the Changbai Mountain region, which were illustrated in Figure 2B-D, respectively.
The national ecological suitable area was 86.87 × 10 4 km 2 and the core suitable area was 10.13 × 10 4 km 2 accounting for 11.66% of the total potential habitat ( Table 2). The suitable distribution of the Loess-Inner Mongolia plateau region was larger than other two regions, reaching 61.39 km 2 . The core ecological suitable area in this region accounted for 13.04% of the total suitable area, and the percentage was higher than the national average.
The Xingjiang region had the smallest ecological suitable area, with the area size of 5.65 km 2 . The proportions of the moderate and core suitable habitat were 19.25% and 3.89% respectively, which were significantly lower than the national average.   Under the future climate scenario, the variation trends of the suitable habitat area and the pattern of increase-decrease were not exactly the same among the three regions. However, by the period of the 2070s, the area of the newly increased and disappeared suitable habitat in the three regions gradually decreased compared with that in the previous period.
As far as the Loess-Inner Mongolian plateau region was concerned, in the 2050s relative to the current distribution, the area of the potential habitat increased to 74.04 × 10 4 and 72.10 × 10 4 km 2 ( Table 3). The new areas were mainly distributed around the unchanged area, while the disappeared habitat is mainly in the west and northeast of the area under the two climate scenarios (Figure 3(A1,A3)). By the 2070s, the total habitat area decreased to 68.38 × 10 4 and 58.43 × 10 4 km 2 , respectively (Table 3), and the disappeared and new regions will be distributed almost randomly around the suitable regions relative to the 2050s (Figure 3(A2,A4)). For Xinjiang region, the new habitat area was in the northwest of Xinjiang in 2050s under the two scenarios (Figure 3(B1,B3)), amounting to 4.01 × 10 4 and 3.42 × 10 4 km 2 (Table 3), and the area of the ecological suitable raised 2.67 × 10 4 km 2 with the percentage of 47.16 under RCP 2.6 in 2050s relative to current distribution, 1.89 × 10 4 km 2 for RCP 8.5. The area of ecological suitable loss was 1.35 × 10 4 and 1.53 × 10 4 km 2 , respectively, projected for the south of the Tianshan mountains (Figure 3(B1,B3)). By the 2070s, the total habitat area went down to 7.90 × 10 4 and 5.6 × 10 4 km 2 , and the disappeared and new regions appeared in a sporadic pattern throughout the entire distribution area (Figure 3(B2,B4)).

Varied Spatial Shifts of Core Distribution in China under the Future Climate Scenario
In  Figure 4(B1,B2)).  Figure 4(A1,A2)).  8.5. Similarly, the elevation of the geographical centroid for core area showed an escalating trend and then task back under the two scenarios. The initial altitude was 2138 m, and in the 2050s descended to 2004 and 2099 m, respectively, afterwards ascending to 2070 m for RCP 2.6 and 2104 m for RCP 8.5 (Table 4, Figure 4(B1,B2)).
There was a difference in the Changbai mountain region compared with the other two regions in terms of migration direction and elevation change. The geographical centroid migrated southwesterly without changing direction, and the elevation raised continuously over the future two periods under the two scenarios.

Suitable Distribution in Current Climate Scenario
In this study, we calibrated the MaxEnt model and projected the distribution regions of J. rigida on the national scale under the current and future climate scenario base on the climate, soil, and topographical variables. Overall, the model showed a perfect performance under all the climate scenarios [34], with the AUC 0.928 ± 0.025.
According to the estimated result, the distribution regions of Juniperus rigida Sieb. et Zucc. concentrated on three regions under current and future climate scenarios, which were renamed Xinjiang region, Loess-Inner Mongolian plateau region, and Changbai Mountain region from west to east (Figure 2). The result was in general agreement with the previous occurrence recorded in nine provinces, including Heilongjiang, Jinlin, Liaoning, Inner Mongolian, Shaanxi, Shanxi, Gansu, Ningxia, and Hebei [35]. There was ca. 86.87 × 10 4 km 2 area for the projected suitable habitat of J. rigida in China, and the largest area was located in the Loess-Inner Mongolian plateau region, whose area was ca. 61.39 × 10 4 km 2 with 70.59% of the national ecological suitable area ( Table 2). There was a region that we should emphasize, the Xingjiang region. Although there were only three natural occurrence records ever in this region, we projected about 5.65 × 10 4 km 2 suitable distribution in this region. In the background of global climate warming, the arid climate and less rainfall coincided with its physiological characteristics [35] and provided inhabitable conditions for sprouting and growth of J. rigida in this region [36]. Furthermore, the Changbai Mountain region was also the major distribution area, which can be verified by the previous studies on species distribution to some degree [37].

Ecological Variables Influencing the Spatial Distribution Pattern of J. rigida
The geographical spatial distribution of plants is the result of long-term interaction between plants and the environment. Nevertheless, environmental factors do not affect the geographical spatial distribution equally. There were four identified variables that played a significant role in influencing the distribution pattern among all the factors that characterized the ecological habitat of Juniperus rigida Sieb. et Zucc. by running the model. The properties of J. rigida, light-loving and drought-resistant [35], made the specific requirements for slope (Slo) and altitude (Alt) in different distribution regions. The distribution elevation was 300~600 m in the Changbai mountain region, which was the lowest distribution area, and was ca. 1100~1600 m in Yinshan mountain, 1100~1300 m for the Erdos plateau, and 2000~2500 m for Helan mountain [35]. The complex topography was responsible for the diverse climate on a regional scale, and the slope affected the light intensity, evaporation, soil nutrients, and available water content to some extent.
The two factors associated to the temperature, mean temperature of the driest quarter (MTDq) and temperature annual range (TAR), also played important roles in affecting the distribution. It is stated that temperature was the dominant factor in the progress of the vegetative growth [38], flower bud differentiation [39], and seed dormancy and germination [40]. Temperature could germinate seeds and shortened the dormancy period caused by the encased in hard seed coat and oil [41]. The optimum germination temperature for Juniperus rigida Sieb. et Zucc. seeds was 15~25 • C, and the light could inspire its germination rate [19]. Global warming was accompanied by rising precipitation [42]. The effect of precipitation on the geographic distribution should not be ignored when we focus on the effect of temperature. We projected that the precipitation of the wettest month (PWm) and annual precipitation (AP) were also important factors in terms of percent contribution and permutation importance of all variables.

Prospective Change on Distribution in the Face of Climate Warming
With global warming, it was reported that species will either migrate [43] or change their physiology [44] to resist climate change. There may be a main trend of species migrating to high latitude or high elevation [45,46]. The mean altitude of the core suitable area gradually increased in Chanbai mountain regions, first up and then down in other two regions, over time, which may be a strategy for Juniperus rigida Sieb. et Zucc. to adapt to the changing environment (Table 3). In this respect, our results were completely consistent with those of previous studies [47]. The migration paths were almost the same in the same study area under the two different climate scenarios; nevertheless, they moved in different directions in different regions (Figure 4). The trajectory of core area centroid was triangular in the Loess-Inner Mongolian Plateau Region, first northeast and then southwest, and almost linear in the Xinjiang Region and Changbai Mountain region, northwest and southwest direction, respectively (Figurre 4). This phenomenon was mainly related to the complexity of climate change. The centroid of Changbai Mountain distribution migrated to low latitude that was different from the previous work [45,46]. The impact of global climate change on China was more obvious than other countries, and the temperature in Northeast of China rose faster than that in other regions [48]. All of those could clarify why the distribution centroid in Changbai Mountain migrated to lower latitudes rather than higher ones. In any region, the centroid shifted further under RCP 8.5 to RCP 2.6 scenario. RCP 8.5, a very high GHG emission scenario, led to robust global surface warming [22][23][24][25], affecting more prominently the distribution of J. rigida in China to RCP 2.6.
In addition to the migration of suitable areas, climate change also affected the area of suitable distribution. On the side, it is also affected the adaptability of species to climate change [49]. We gained the similar achievement of previous studies that showed that the suitable habitat range will increase as global warming intensity proceeds under the low concentration emissions pathway (RCP 2.6) [50], but excessive temperature raised might restrict plants' diffusion under the powerful GHG emissions pathway (RCP 8.5) [51]. The suitable area increased and then decreased in the Xinjiang region and the Loess-Inner Mongolian plateau region and decreased continuously in Chanbai Mountain region in both scenarios. This may be caused by the more significant effect of climate warming and precipitation in the northeast of China than other regions [48,52]. Comparing the effects of climate change and precipitation on China [53], the change in temperature and precipitation were marked in the following areas where the area of disappeared suitable habitat were severe, including the north part of Chanbai Mountain and the south part of Taihang Mountain and Lvliang Mountain. On the contrary, the adjacent regions of Hebei, Inner Mongolia, and Shaanxi province were flooded with new distributions in the 2050s and 2070s, in which there were mild changea in temperature and precipitation under the stringent concentration emissions pathway (RCP 2.6) and the high greenhouse emission scenario (RCP 8.5) [48].

Conclusions
There, we projected the current and future distribution areas of Juniperus rigida Sieb. et Zucc. in China using the MaxEnt model under the climatic scenario of RCP 2.6 and RCP 8.5 for two periods, the 2050s and 2070s. These results indicated that the temperaturedependent variables affected the suitable habitat more than the factors associated with the precipitation. The suitable area increased first and then decreased on the whole, and the centroid of the core habitat migrated to high elevation. These insights were remarkable in the Changbai Mountain region under the same scenario and were more obvious under RCP 8.5 in any region. Though we gained meaningful results that will provide a more effective and evidence-based conservation strategy, they were not comprehensive, and other factors should be taken into account too, for instance, adaptability to climate change, diffusivity, and so on. However, these were unfamiliar to us in this paper.