Interspecific Variance of Suitable Habitat Changes for Four Alpine Rhododendron Species under Climate Change: Implications for Their Reintroductions for

: Rapid temperature changes in mountain ecosystems pose a great threat to alpine plant species and communities. Rhododendron species, as the major component of alpine and sub ‐ alpine vegetation, have been demonstrated to be sensitive to climate changes. Therefore, understanding how alpine Rhododendron species spread to new habitats and how their geographical distribution range shifts is crucial for predicting their response to global climate change and for facilitating species conservation and reintroduction. In this study, we applied MaxEnt modeling and integrated climate, topography, and soil variables in three periods under three climate change scenarios to predict the suitable habitat for four Rhododendron species in China. We measured the potential distribution change in each species using the change ratio and the direction of centroid shifts. The predicted results showed that (1) the threatened species R. protistum would have a maximum decrease of 85.84% in its distribution range in the 2070s under RCP 8.5, and R. rex subsp. rex as a threatened species would experience a distribution range expansion (6.62–43.10%) under all of the three climate change scenarios in the 2070s. (2) R. praestans would experience a reduction in its distribution range (7.82–28.34%) under all of the three climate change scenarios in the 2070s. (3) The four Rhododendron species would be moved to high latitudes in the north ‐ westward direction as a whole in the future, especially the two threatened species R. protistum and R. rex subsp. rex . (4) Aside from climate variables, soil factors also exert an important influence on the distribution of Rhododendron species. This study revealed the species ‐ specific response of Rhododendron species to climate change. The results can not only provide novel insights into conservation strategies of Rhododendron species, but also propose a valuable method for the habitat selection during the reintroduction of endangered species. under Change: Implications for


Introduction
Human-induced global changes are the most severe phenomena since the Last Glacial Maximum. Rapid global climate warming has led to distribution pattern changes for numerous species since the original habitats of these species may no longer be appropriate for their growth and survival [1]. Climate change plays a driving role in biodiversity loss, changes in the spatial patterns of species, and threatened species' survival, and it can increase the risk of extinction for endangered plants. Numerous studies have shown that the habitat range of most species is decreasing [2,3], and the geographic distributions of these species are expected to move towards high altitudes and high elevations to expand to new favorable areas under climate warming [1]. Moreover, growing evidence shows that the speed of warming is amplified by elevation under the background of global warming, which means that temperatures change more rapidly at high elevations than at low elevations [4]. This rapid change in temperature poses a great threat to alpine plant species and communities [5,6]. Therefore, understanding how alpine species spread to new habitats and how the geographical distribution range changes is crucial for predicting their response to global climate change and facilitating species conservation and reintroduction [7].
Rhododendron L., which comprises about 1025 species, is one of the largest woody plant genera. This genus is widely distributed from the northern temperate region throughout southeast Asia to northeastern Australia [8,9]. About 571 Rhododendron species exist in China, and 70% of them are endemic in the southwestern region of China. Therefore, the mountain land of Southwest China is considered the center of origin or evolution of modern Rhododendron [8]. However, with the ever-increasing impact of human activities and global climate change, the survival of Rhododendron resources in China is facing a serious threat [10]. In addition, as the major component of alpine and sub-alpine vegetation, the Rhododendron species have been demonstrated to be sensitive to climate change [11][12][13]. Yu et al. (2019) [12] applied the presence-only ecological niche model MaxEnt and predicted the distribution of 10 narrow-ranging and 10 wide-ranging Rhododendron species. The authors discovered the negative effect of climatic and land use changes on the distribution of species. Lu et al. (2021) [13] assessed the effect of protected areas and tourist attractions on the current and future suitable habitat distributions of seven species of the genus Rhododendron in Northeast China and provided conservation suggestions.
Species distribution models (SDMs), as an effective tool in predicting species' suitable habitats under climate change, play important roles in studying the processes of species' ecological evolution and in conservation planning [13,14]. On the basis of species occurrence data and environmental variables, SDMs can predict the potential distribution of species at different spatial and temporal scales, especially under scenarios of global climate change [15]. At present, SDMs based on different theoretical and analysis methods, such as surface range envelope (BIOCLIM), artificial neural network (ANN), maximum entropy (Maxent), random forest (RF), generalized linear model (GLM), and generalized additive model (GAM), are widely used in ecology and biogeography. The Maxent model is widely recognized as one of the most effective tools for predicting species' suitable habitats and changes in future distributions due to its advantages. In particular, it requires only species presence data, can deal effectively with limited occurrence data and small sample sizes, and can use continuous and categorical environmental data as input variables. Numerous studies have predicted changes in animal and plant distributions in response to contemporary climate changes using the Maxent model [7,[16][17][18][19][20][21].
Reintroduction is one of the rescue methods of plant conservation, especially for endangered species threatened by habitat destruction and loss under global changes. Predicting habitat suitability and the potential habitat distribution range is critical to the reintroduction of endangered plants [22,23]. Here, we integrate climate and soil variables to predict the suitable habitat for four Rhododendron species in China. We hypothesize that the suitable range for the four Rhododendron species would contract and shift to high altitudes in the future since the original distribution areas may become climatically unsuitable under climate change omission scenarios. Furthermore, we predict that interspecific variance would occur in the future. The results of this study can provide novel insights into the creation of conservation strategies and habitat selection during the reintroduction of Rhododendron species.  [8]. Rhododendron protistum is listed as a critically endangered species in the China Plant Red Data Book [24]. Four Rhododendron species are evergreen trees that are distributed in forest ecosystems in the southwest plateau mountains in China.

Species Occurrence Data
Occurrence data were obtained from published studies, the China Digital Herbarium (http://www.cvh.ac.cn/, accessed on 1 October 2021), and GBIF (https://www.gbif.org/, accessed on 1 October 2021). In accordance with the biological and ecological characteristics of these species, we removed unreasonable sites to ensure geographic accuracy. Each sample point was at least 2 km apart to reduce the effect of spatial autocorrelation. Then, we obtained 24, 22, 15, and 39 effective localities of Rhododendron protistum, Rhododendron rex subsp. rex, Rhododendron praestans, and Rhododendron sinogrande in Southwest China for the subsequent model construction and analysis ( Figure  1).

Environment Data
Climate variables are the major predictors in the MaxEnt model since they are the primary factors regulating species' geographic distributions. However, recent studies have found that other factors, including soil types and characteristics, land use changes, and vegetation types, have a significant influence on distribution predictions. These factors are often correlated with climate change [15,16]. Therefore, combining climate and other factors that affect species' biological characteristics in MaxEnt model building can reveal a suitable habitat under the scenarios of global climate change. This suitable habitat can be used for species reintroduction and rehabilitation.
The environment data included 19 climate variables, 3 topography variables, and 16 soil variables. Climate data (Bio1-Bio19) on current climate conditions , the 2050s (average for 2041-2060), and the 2070s (average for 2061-2080) were downloaded from the Worldclim Database (https://www.worldclim.org/, accessed on 1 October 2021) at a resolution of 30 arc-seconds (~1 km 2 ). For future climatic projection, we used the Beijing Climate Center Climate System Model (BCC-CSM1-1), which is a widely utilized global climate model (GCM) in China provided by IPCC5 [25]. The representative concentration pathways (RCPs) form a set of greenhouse gas concentration and emission pathways which are widely used to research the impacts and potential policy of the responses to climate change. For this GCM, we selected three RCPs of emission scenarios from low to high, namely, RCP 2.6, RCP 4.5, and RCP 8.5. Furthermore, we used topsoil variables at depth intervals of 0-30 cm in Maxent model building for precise suitable habitat prediction. The topsoil variable data were derived from the Harmonized World Soil Database (http://www.iiasa.ac.at/web/home/research/researchPrograms/water/HWSD.html, accessed on 1 October 2021) at the same resolution (1 km 2 ). Digital elevation model (DEM) data were obtained from the SRTM Database (http://srtm.csi.cgiar.org/, accessed on 1 October 2021) with a 90 m resolution. Meanwhile, the slope and aspect were extracted from DEM data in ArcGIS 10.2 and resampled into a 1 km spatial resolution.

Study Area and Environment Variables Selection
To approach a more realistic prediction, we delimited the study area of each Rhododendron species by drawing a minimum convex polygon (MCP) with a buffer distance of 1 degree (~111 km) around the localities and obtained the environment factors within this region by the Wallace package in R. All of the analyses include variables and the selections were accomplished within the buffer areas. In addition, we explored the potential distribution of four Rhododendron species in Southwest China.
The collinearity among environmental factors causes overfitting and uncertainty in the predicted results [26]. To establish a good performance model with few variables for four Rhododendron species, first, we used 38 variables to prebuild the MaxEnt model three times in succession and eliminated the environment factors with no contribution or importance to the models. Moreover, we eliminated four variables that combined temperature and precipitation since they have artificial discontinuities [27,28]. Therefore, we performed a correlation analysis by the Band Collection Statistics tool of ArcGIS for the remaining environment variables, except for the categorical factor. In addition, we retained one factor which has a high contribution from each set of highly cross-correlated variables (R 2 > 0.8) for further modeling (Tables 1, 2 and A1-A4).

Model Evaluation and Model Selection
When selecting the model parameters, we used models with regularization multiplier (RM) values ranging from 0.5 to 4.5 (increments of 0.5) and with six feature class (FC) combinations (L, LQ, H, LQH, and LQHP, with L = linear, Q = quadratic, H = hinge, and p = product). In accordance with the lowest delta, the AICc score and the predictor variables (RM and FC) were used in the final model. All of the analyses were performed using the Wallace package in R [29][30][31].
We inputted the species and environment data into MaxEnt 3.4.1 [32]. The feature combination and regularization multiplier were set based on Table 3. For each run, we set the output format as "Cloglog" and the output file type as "asc." The model was estimated using 75% of the occurrence data for model calibration, and 25% was used for model testing. The replicate was set as "10" with a replicated run type as "Crossvalidate" to reduce uncertainty. The remaining parameters adopted default settings. The average of 10 predictions was taken as the result. The area under the ROC curve (AUC) and the true skill statistic (TSS) are the two metrics that we used to measure model performance. The range of the AUC is from 0.5 to 1 and TSS ranges between −1 and +1. The closer the assessment value is to 1, the better the model performance [33,34].

Geospatial Analyses
To estimate the similarity of environment conditions between the buffered background area and projection area and to know where the climates are novel, we calculated the multivariate environmental similarity surfaces (MESS) following the guideline from Elith et al. [35]. The novelty of environment (i.e., extrapolation, negative values in MESS result) indicated where the prediction areas are questionable and extrapolated in ecologically unrealistic ways [35][36][37]. Therefore, we masked out these extrapolation areas (negative values in the MESS map output). As a result, a predictive map with four classes of habitat suitability for the four Rhododendron species was defined using the reclass tool in ArcGIS 10.2. The four classes were unsuitable (<0.10), low suitability (0.10-0.33), moderate suitability (0.33-0.67), and high suitability (>0.67). To measure the predicted distribution changes and centroid shifts for each species, we projected the binary species distribution models (SDMs) onto WGS_1984_UTM_Zone_47N projection in ArcGIS 10.2. Then, we calculated the change ratio between the current and future distribution and centroid change direction vectors of the current and future binary SDMs using the SDM Toolbox [38].

Model Selection and Evaluation
In accordance with the lowest delta.AICc, we obtained optimal parameters for each Rhododendron species (Table 3). On the basis of these parameters, current suitable habitat distributions of R. protistum, R. rex subsp. rex, R. praestans, and R. sinogrande in the current period were constructed. The AUC values for the models of the four Rhododendron species were between 0.920 and 0.988, and the TSS ranged from 0.66 to 0.94 (Table 3).

Current Habitat Distribution and Dominant Environment Variables
After masking out the outlier areas which have negative values in MESS ( Figure A1, Table A5), about 10% to 20% of projection grid cells were predicted as a suitable potential area under current climatic conditions. The current potential distribution prediction of the four Rhododendron species (Figure 2) indicated that R. protistum is mainly distributed in the Gaoligong Mountain of northwestern Yunnan Province in China. In addition, Tibet bordering northern Myanmar has a small suitable area. R. rex subsp. rex has a wide, highly suitable distribution range and is mainly distributed in northeastern Yunnan Province and southwestern Sichuan Province. R. praestans has a major current habitat suitability in northwestern Yunnan Province and southeastern Tibet. Moreover, the Hengduan Mountainous Region is a potential current habitat. The current habitat distribution range of R. sinogrande is similar to R. protistum. In other words, it is mainly distributed in western Yunnan Province and southeastern Tibet (Figure 2). The current distribution areas of R. protistum, R. rex subsp. rex, R. praestans, and R. sinogrande are 10,440.58, 178,394.88, 132,745.11, and 38,819.5 km 2 , respectively. The jackknife test on the variables' contributions to Maxent revealed the influence of the environmental factors on the spatial distribution of the species. High values of the regularized training gain indicate a highly significant contribution of variables [21]. The results of the jackknife test on the regularized training gain for the four Rhododendron species are shown in Figure A2. For R. protistum, Bio13 provided very high gains (>1.0) when utilized independently, indicating that the precipitation of the wettest month is the dominant factor in habitat conditions. Bio15, Bio17, Bio3, Bio4, T_CLAY, T_GRAVEL, and T_OC provided a moderate gain when used independently, suggesting that they can affect the habitat suitability distribution of R. protistum. For R. rex subsp. rex, only Bio4 and soil type (T_USDA_TEX_CLASS) provided a moderate gain when used independently. Meanwhile, Bio15, Bio2, Bio7, and the slope provided moderate gains for R. praestans when used independently. Bio2 and Bio4 provided high gains (>1.0) for R. sinogrande, and Bio13, Bio17, Bio3, Slope, T_CLAY, and T_USDA_TEX_CLASS provided moderate gains when used independently.

Projected Changes in Habitat Suitability for Future Periods
Prediction of the potential distribution of four Rhododendron species in China under different future climatic scenarios (RCP 2.6, RCP 4.5, and RCP 8.5) in the 2050s and 2070s are shown in Figure 3. Compared with the current high-suitability habitat distribution range (Figure 4) in the 2070s, R. protistum will show a decreasing-increasing-decreasing trend with intensified emission scenarios, and it will have a maximum decrease of 85.84% in the 2070s under RCP 8.5. R. rex subsp. rex would experience a distribution range expansion under all of the three scenarios in the 2070s. However, R. praestans will show an opposite performance as R. rex subsp. rex. In addition, the distribution range of R. sinogrande will decrease under RCP 2.6 and expand under RCP 4.5 and RCP 8.5.  After a comparison of the high-suitability habitat change ratio of each species between two periods (current-2050s and 2050s-2070s) under RCP 2.6, RCP 4.5, and RCP 8.5 ( Figure 5), we found that the high-suitability habitat of R. protistum and R. rex subsp. rex will increase in the two periods under RCP 4.5. However, R. protistum will exhibit a trend of shrinkage-expansion and expansion-shrinkage under RCP 2.6 and RCP 8.5, respectively. The high-suitability habitat will decrease sharply in the 2050s-2070s under RCP 8.5. Meanwhile, the high-suitability habitat of R. rex subsp. rex will increase by at most 44.38% within the current-2050s under RCP 2.6. The high-suitability habitat of R. praestans will increase first, then decrease under RCP 2.6 and RCP 4.5, and it will show a persistent shrinkage trend under RCP 8.5. On the contrary, the high-suitability habitat of R. sinogrande will decrease first, then increase under RCP 2.6 and RCP 4.5.

Shifts of the Centroids of the High-Suitability Habitat in the Future
The overall trend of the centroid change of R. protistum will continue to be at high latitudes in the period current-2050s and 2050s-2070s under the three scenarios ( Figure  6a), indicating that the climate sensitivity of the species will survive. R. rex subsp. rex will migrate to the northwest of the high latitude under RCP 4.5 and RCP 8.5, but under RCP 2.6, its centroid will move to low and then to high latitude (Figure 6b). For the centroid shift of R. praestans, it will move to high latitude, except for the period current-2050s under RCP 8.5 (Figure 6c). However, R. sinogrande is different from the former species, its centroid will change to high latitude first and then to low latitude under RCP 2.6 and RCP 4.5 (Figure 6d).

Discussion
MAXENT is reliable in predicting the potential geographic distributions of species [19][20][21]23]. Ardestani reported that [22] the species with a wide distribution range tend to have low AUC values. Our results on the four Rhododendron species' AUC values support this viewpoint. From small to large, the current habitat distribution ranges of the four focused species are in the order of R. protistum, R. sinogrande, R. praestans, and R. rex subsp. rex, and the AUC values are 0.996, 0.975, 0.922, and 0.918, respectively. In other words, R. protistum has the smallest distribution range with the largest AUC value (AUC = 0.996), and R. rex subsp. rex has the widest distribution range and the smallest AUC value.
The niche theory suggests that the species with a limited distribution range size often has a narrow niche breadth and low tolerance and adaptability to changeable environments [39]. The species with a wide distribution range tend to have a broad niche breadth, which allows them to adapt to climate change [13]. In the present study, the geographical ranges of R. protistum and R. sinogrande were locally distributed. These species have strict requirements for habitat environments. In other words, the predicted probability of presence will be greater than 0.632 if precipitation or other environment variables reach a certain value [40]. Moreover, the distribution sizes of R. rex subsp. rex and R. praestans were larger than the former. From the jackknife test on regularized training gain, we found that R. rex subsp. rex and R. praestans have a broad range of environment requirements. However, the suitable habitat proportion of limited-distribution plants R. protistum and R. sinogrande did not decrease under different future climate conditions according to the niche theory. On the contrary, their suitable habitat area will increase in the future. This prediction is consistent with the results of Yu et al. (2019) [12] and Lu et al. (2021) [13], in which the narrow-ranging species may experience range expansion under climate change. It is possible that the species occupying a warm climatological niche will benefit from climate change [41]. In addition, the distribution range of R. praestans was wider than the former, and its habitat size will decrease in the future. This result is consistent with the result of Yu et al. (2019) [12], in which the wide-ranging Rhododendron species will be adversely affected by environment changes. However, the distribution shift trend of R. rex subsp. rex was contrary to R. praestans. These species-specific responses of Rhododendron species imply specific conservation and reintroduction strategies for different species.
In general, the species would tend to distribute to high latitudes or elevations in the future to achieve a suitable habitat under continuous climate warming. Previous studies have demonstrated that numerous species change their geographic distribution range towards high latitudes and altitudes [1,16,17]. However, an increasing number of reports have found other types of range changes, such as east-west directions across longitudes or towards low elevations and tropical latitudes [42]. In this study, the four Rhododendron species were predicted to move towards high latitudes in the north-westward direction as a whole, especially the two threatened species R. protistum and R. rex subsp. rex. These results are in line with the universal migration rule of alpine plant species. However, the different species showed different movement trends (towards low latitudes or east-west direction across longitudes) under the different climate scenarios. For instance, the centroid shifts of the highly suitable habitat of R. sinogrande were predicted to move to high latitudes and to low latitudes in RCP 2.6 and RCP 4.5, respectively. This result implies that the species R. sinogrande has stronger adaptability to climate warming than the other species.
An organism's habitat is the combination of the space it occupies and all of the ecological factors in that space, including the abiotic environment and the biotic environment, which are necessary for the survival of the other organisms [21]. Therefore, before conducting the predictions, we constructed a buffer area around the species occurrence localities to delimit the study region. This can greatly improve the accuracy and truthfulness of the predictions since the area includes environments that are accessible to the species, given its limitations and the configuration of barriers [43]. If possible, we could base this on a biological justification of the factor, physiological information from laboratory experiments, and natural history information when selecting the predictor variables to make the predictions realistic [44]. Previous studies have reported that aside from climate variables, other environmental factors, such as soil, topography and land use, vegetation dynamics, and inter-specific competition, affect the response of Rhododendron species to climate change and their future distribution [11,45,46]. Feng et al. (2020) [16] built ecological niche models using climate and soil variables when predicting the habitat distribution of Camptotheca acuminate. The authors found that soil factors can accurately limit the distribution range of species on the basis of climate factors, and the accuracy of climate and soil models for C. acuminata is high. Abdelaal [23] suggested that combining climate variables with soil variables can predict the distribution of species accurately. In the present study, climate, topographic, and soil variables were combined to explore the potential habitat distribution shift in current and future periods of four alpine Rhododendron species under rapid climate change. The results of the regularized training gain of each species showed that soil factors affect the distribution of Rhododendron species. T_CLAY, T_GRAVEL, and T_OC are important factors in predicting the distribution of R. protistum. T_USDA_TEX_CLASS has a moderate influence on R. rex subsp. rex. T_CLAY and T_USDA_TEX_CLASS can affect the habitat distribution of R. sinogrande. The clay content (T_CLAY) provides a high soil organic carbon, and the gravel content in the topsoil (T_GRAVEL) affects the respiration and water absorption of the roots. The content of topsoil organic carbon (T_OC) is a significant indicator in evaluating soil fertility. Moreover, the soil type (T_USDA_TEX_CLASS) can reflect the soil physical properties. Overall, we reasonably speculate that the nutrient-holding and water-holding capacities of soil significantly affect the distribution of Rhododendron species.
Rhododendron species have high ecological and ornamental values, but 70% of the plants in this group are classified as vulnerable, threatened, endangered or critically endangered [9,10]. Therefore, conservation recommendations must be made in accordance with the characteristics of the species, especially under future climate change scenarios. SDMs are easy-to-use tools for selecting the probable distribution and areas suitable for the reintroduction of endangered species [20,21,23]. Moreover, our results indicate that aside from climate factors, the contents of clay, gravel, and organic carbon in the topsoil have a significant influence on the endangered plant R. protistum. This result implies that the conservation and future reintroduction of R. protistum should consider not only the climate factors, but also the characteristics of topsoil in the micro-habitat. In addition, we found that R. protistum will face a severe contraction across its highly suitable distribution area in the 2070s under the RCP 8.5 climate change scenario. Notably, the remaining habitat area is still concentrated in the national natural conservation zone in the Gaoligong Mountain. Therefore, we reasonably speculate that the survival of R. protistum depends on the specific environment or local microhabitats, and R. protistum may have nearly no adaptive capacity to climate change. We recommend protecting R. protistum in situ where its preference is. In addition, we recommend that the habitat selection of future reintroduction and population recovery of R. protistum should focus on the Gaoligong Mountain. R. rex subsp. rex is an endangered plant endemic to China. However, the present study showed that this species will experience distribution range expansion under all of the three scenarios in the 2070s. This finding indicates that R. rex subsp. rex may have a certain adaptive capacity to future global warming. Meanwhile, a previous study has reported that the main reason for the endangerment of R. rex subsp. rex is exotic anthropogenic disturbance [8,47]. Furthermore, the present study found that the distribution of R. rex subsp. rex is mainly affected by the temperature seasonality standard deviation and soil type. Therefore, we suggest that R. rex subsp. rex should be given priority by building a conservation area to protect it from high-frequency human disturbances and implement reintroduction strategies for its future suitable habitat. Although R. praestans and R. sinogrande are not listed as threatened plants in China, the distribution size of R. praestans will decrease by 28.34% in the future. Therefore, we suggest that comprehensive studies and additional other alpine Rhododendron species should be given urgent attention under the scenario of climate change.

Conclusions
Global climate change causes distribution changes in organisms, which significantly affect biodiversity conservation and ecosystem safety. Therefore, predicting the trends of distribution changes of plants can provide a basis for plant conservation and habitat selection in species' reintroduction strategies. The present study applied MaxEnt modeling and integrated climate, topography, and soil variables in three periods under three climate change scenarios to predict the suitable habitat for four Rhododendron species in China. We revealed species-specific responses to climate change and the influence of topsoil factors on the distribution of Rhododendron species. Moreover, we presented corresponding conservation strategies for different species.

Conflicts of Interest:
The authors declare no conflict of interest.