Predicting range shifts of Davidia involucrata Ball. under future climate change

Abstract Understanding and predicting how species will respond to climate change is crucial for biodiversity conservation. Here, we assessed future climate change impacts on the distribution of a rare and endangered plant species, Davidia involucrate in China, using the most recent global circulation models developed in the sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC6). We assessed the potential range shifts in this species by using an ensemble of species distribution models (SDMs). The ensemble SDMs exhibited high predictive ability and suggested that the temperature annual range, annual mean temperature, and precipitation of the driest month are the most influential predictors in shaping distribution patterns of this species. The projections of the ensemble SDMs also suggested that D. involucrate is very vulnerable to future climate change, with at least one‐third of its suitable range expected to be lost in all future climate change scenarios and will shift to the northward of high‐latitude regions. Similarly, at least one‐fifth of the overlap area of the current nature reserve networks and projected suitable habitat is also expected to be lost. These findings suggest that it is of great importance to ensure that adaptive conservation management strategies are in place to mitigate the impacts of climate change on D. involucrate.

. Thus, in order to mitigate the negative effects of climate change on species, conservation strategies should be refined by modeling species distributions to identify to what extent they could be influenced by future climate change (Ramirez-Villegas et al., 2014;Thuiller et al., 2008).
In the past two decades, species distribution models (SDMs) have been widely used to assess the impacts of future climate change on species distributions and guide conservation planning (Kujala et al., 2013;Maggini et al., 2014;Wiens et al., 2009). However, SDMs may suffer from a lack of precision and portability due to the variation in covariate selections (Zhang & Zhang, 2012), type of SDM used (Hartley et al., 2006;Pearson et al., 2006;Thuiller et al., 2009;Wenger et al., 2013), and climate projections arising from different global circulation models (GCMs) and CO 2 emission scenarios (Barry & Elith, 2006;Wenger et al., 2013), which can yield misleading or inconsistent outcomes, posing challenges for decision-making . Ensemble modeling approaches, which combined a series of SDMs, can produce consensus projections that may outperform single SDMs (Araújo & New, 2007;Marmion et al., 2009) and reduce the predictive uncertainty of single algorithm by combining their predictions (Tehrani et al., 2021;Thuiller et al., 2014). By using ensemble modeling approaches, more robust projections can be produced with reasonable interpretation (Araújo & Guisan, 2006).
Davidia involucrata Baill., commonly known as dove tree or handkerchief tree, is a rare and endangered species listed in the China Plant Red Data Book under first-grade state protection (Liu et al., 2019). It is also a Tertiary relict plant endemic to China (Fu & Jin, 1992), currently ranges approximately from [98][99][100][101][102][103][104][105][106][107][108][109][110][26][27][28][29][30][31][32]including Yunnan,Guizhou,Sichuan,southern Shaanxi,southern Gansu,Chongqing,Hubei,and Hunan Provinces (Li, 1954;Liu et al., 2019;Takhtajan, 1980;Tang et al., 2017). Its populations are often found in subtropical evergreen broad-leaved forests or in mixed forests of temperate deciduous broad-leaved trees at altitudes of between1100 and 2,600 m (He et al., 2004). Owing to the highly strict ecotope and recruitment limitation (i.e., low reproduction rate and limited dispersal ability), the population age structure of D. involucrate is declining (Wang Yusheng et al., 2019). In addition, the increasing intensity of human activities (e.g., logging) has led to a sharp decrease of its remaining habitats (Wang Yu-sheng et al., 2019). Despite its threatened status, few studies have explored the vulnerability of D. involucrate to climate change (Tang et al., 2017;Wang Yu-sheng et al., 2019). By using ecological niche models (ENMs) (Boria et al., 2014;Peterson, 2006), Tang et al. (2017) projected the potential suitable habitats of this species under past, current, and future climatic conditions. In their work, the obsolete CMIP5 climate models were used to simulate future climate conditions (Tang et al., 2017). However, many studies have shown that the most recent CMIP6 climate models perform better in the simulation of future climate conditions compared to the CMIP5 climate models (Fan et al., 2020;Xin et al., 2020). Therefore, a rigorous analysis combined ensemble modeling approaches and the CMIP6 climate models investigating potential impacts of future climate change on the distribution of D. involucrate is of great urgency and significance.
In this study, we aim to (a) assess the vulnerability of D. involucrate to climate change (i.e., whether the suitable habitats of D. involucrate will suffer great lost under future climate change), and (b) evaluate the conservation effectiveness of current nature reserve networks in protecting D. involucrate under climate change. To this end, we compiled a large dataset on spatially explicit species presence records of the D. involucrate and environmental data (bioclimatic variables) covering China and subsequently used ensemble modeling approaches to project the potential suitable habitats of D. involucrate under current and future climatic conditions. According to our knowledge, our study is one of the first studies to investigate how D. involucrate will response to future climate change by using ensemble modeling approaches and the most recent CMIP6 global circulation models.

| Study area and species occurrence data
According to previous studies on the suitable habitat of D. involucrate, we chose the whole China as the study area ( (see below), we excluded the occurrence records that were not collected from this time period. Furthermore, to reduce potential errors in species' geographic locations of occurrence data, we only include species occurrence records with geographic locations. Finally, to ensure that all occurrence records were in the species' native geographic locations, we also excluded the species occurrence records collected in manual intervention areas, such as parks and experimental forests. From this process, we had 337 occurrence records available to use (Figure 1).
In order to reduce the high uncertainty in geographic coordinates of the occurrence records and minimize the sampling bias effect in the occurrence records dataset, all the 337 occurrence records were compiled at a spatial resolution of 10 × 10 km grids cell. After removing duplicate records within each gird cell, we obtained 324 presence records to model ecological niches for D. involucrate.
The 19 climate variables used in this study are usually strongly correlated (Marshall et al., 2018). To minimize multicollinearity among variables, we used Pearson's correlations and variance inflation factors (VIFs) to exclude highly correlated variables. Variables with a Pearson correlation >0.70 were considered highly correlated (Dormann et al., 2013), and a VIF >5 was used as a signal that a model had collinearity issues (Rogerson, 2010). Finally, six climate variables were selected for modeling species distributions, including annual mean temperature (BIO1), isothermality (BIO3), temperature annual range (BIO7), precipitation of the driest month (BIO14), precipitation seasonality (BIO15), and precipitation of the warmest quarter (BIO18) (Appendix S2).

| Species distribution modeling
An ensemble of species distribution models (Araújo & Guisan, 2006) was used to model potential suitable habitat for D. involucrata using the biomod2 package in the R platform (v. 4.0.4; http://cran.r-proje ct.org). We chose the ensemble modeling approach because of its ability to create a consensus of the predictions of multiple algorithms and reduce the predictive uncertainty of single algorithm (Kanagaraj et al., 2019;Thuiller et al., 2014;Zhang, Dong, et al., 2020;Zhang, Mammola, et al., 2020). Ten algorithms were considered in the ensemble model: artificial neural network (ANN; Ripley, 1996) et al., 2006), random forest (RF; Breiman, 2001), and surface range envelope (SRE; Busby, 1991). For algorithms requiring species absence records, we generated 10,000 pseudo-absence points with the equal number as the occupied grids for GAM, GBM, GLM, and RF and 10,000 background points for MAXENT by randomly sampling F I G U R E 1 Potential suitable (blue) and unsuitable (gray) habitats suitability of Davidia involucrata Baill. under current climatic conditions in China. Red points represent occurrence records of D. involucrata without replacement. To avoid model over-fitting, the selection of pseudo-absence or background points was limited to biological relevance combinations occupied by the species' range map.
To evaluate the accuracy of each algorithm, we performed cross-validation on each algorithm using bootstrap approach, where random subsets of 80% of each dataset for training data and the remaining 20% for testing algorithm performance using the area under the receiver operating characteristics curve (AUC) and true skill statistics (TSS). This procedure was repeated 10 times to create predictions independent of the training data. Algorithms with AUC >0.90 and TSS >0.80 were considered to have good predictive performance (Allouche et al., 2006;Gallien et al., 2012;Swets, 1988) and were thus kept in the final ensemble model. Based on the selected algorithms in the final ensemble model, we evaluated the relative importance of predictor variables in determining the geographical distribution of D. involucrate by calculating the change in correlation between the covariates and the response with and without the selected variable (Thuiller et al., 2015). The final ensemble model was then projected to current and future climatic conditions by using all occurrence and pseudo-absence data. Finally, these habitat suitability maps were converted to binary presence absence maps using a threshold that maximums model sensitivity plus specificity, which has been shown generally to perform well (Lawson et al., 2014;Liu et al., 2013;Thuiller et al., 2015).

| Statistical analysis
Analyses were conducted on the ensemble model map projections of binary presence absence maps. Firstly, to assess potential impacts of climate change on species ranges, following Zhang et al. (2015), we used two metrics to quantify species' vulnerability: the relative change in total area of suitable habitat (CSH) and the loss of current suitable habitat (LSH). The first metric assumes unlimited dispersal into the projected entire suitable habitats in the future time periods and can be calculated using the following equation: where AREA future and AREA current are the area of current and future suitable habitats. The second metric assumes no dispersal into the projected suitable habitats out of the current suitable habitats and can be calculated using the following equation: Secondly, to detect the direction and distance of species range shifts under future conditions, we determined the centroids of current and future binary presence absence maps using the R package "rgeos" with the "gCentroid" function.
Finally, to explore the conservation effectiveness of current nature reserve networks in protecting D. involucrate under climate change, we also calculated the area of the current and projected suitable habitat overlapped with the current nature reserve networks, respectively.

| Model performance and variable contribution
The AUC and TSS measures provided highly consistent estimates of the model performance of the 10 modeling algorithms (

| Species' range shifts under future climatic conditions
Under the current conditions, the potential suitable area for D.  (Figures 2 and 3a, b). The loss of potential suitable habitats under zero dispersal is more severe than those under global dispersal (Figures 2 and 3c,d), despite having a similar trend in predicted species range size. However, projections indicated that a small part of Gansu, Shaanxi, and Tibet provinces will probably become suitable for D. involucrate in the future (Figure 2).

| The Effectiveness of current nature reserve networks
The climate model and RCP 8.5 scenario) by the 2070s (Figure 5a,b).
Similarly, the overlap area of the current nature reserve networks and projected suitable habitat under zero dispersal is larger than those under global dispersal (Figure 5c,d).

| D ISCUSS I ON
Understanding and predicting how species will response to future climate change is crucial for biodiversity conservation (Wiens et al., 2009) and has required novel approaches with high predictive performance and low predictive uncertainty Thuiller et al., 2014).
In this study, by using ensemble SDMs and the CMIP6 GCMs, we Previous studies have reported that potential suitable habitats of D. involucrate were mainly distributed in mountainous areas with narrow annual temperature range and high precipitation (Liu et al., 2019;Su & Zhang, 1999), reflecting this species cold intolerance (Su & Zhang, 1999  It is often assumed that more complex and more up-to-date models will perform better and/or produce more robust projections than previous-generation models (USGCRP, 2017). Consistent with previous studies Thuiller et al., 2014), our results showed that the ensemble SDMs have higher predictive ability than individual SDMs. Furthermore, our future projection suggested that a large proportion of suitable habitats will be lost under future climate change, which has provided useful dialogue informing conservation strategies for D. involucrate. Our results revealed that the loss of the suitable habitat of D. involucrate would affect the conservation effectiveness of the current nature reserve networks, which does not protect the current suitable habitat of D. involucrate adequately, nor will they protect future potential suitable habitat. These findings also inform a critical initial step in implementing the adaptation planning processes for protecting D. involucrate (Yu et al., 2017;Zhang et al., 2017). On one hand, some nature reserves located in Sichuan Basin would suffer the greatest loss of suitable habitat under future climate change ( Figure 2). These nature reserves are urgent need to be reevaluated under the background of climate change (Bellard et al., 2012;Hansen et al., 2010), and coping strategies to deal with these potential threats require further in-depth study; on the other hand, a part of Sichuan, Gansu, Shaanxi, and Tibet provinces were identified to be the potential climatic refuges (i.e., unchanged and new gained suitability habitat; Ashcroft, 2010) for D. involucrate ( Figure 2). It is urgent needs to design these regions as priority for conservation and to establish new nature reserves in the currently unprotected areas (e.g., southern Shaanxi and southern Gansu) in these regions.
Despite the predictive power of their ensemble modeling of D. involucrate, an important limitation in the present study is that we assessed future habitat suitability under two extreme dispersal assumptions (i.e., no dispersal and unlimited dispersal), which ignores the realistic rates and modes of dispersal of this species (Saupe et al., 2012). These assumptions are likely inaccurate, which could lead to overestimation of suitable habitat under the unlimited dispersal assumption or underestimation under the zero-dispersal assumption (Engler & Guisan, 2009;Viana, 2017;Zanatta et al., 2020).
For instance, Engler and Guisan (2009)  Overall, our research provides fundamental knowledge for understanding the potential impacts of climate change on the distribution of D. involucrate. This study also provides useful information for comprehending vegetation changes at global scales under climate change, especially for the climate-related range shifts of Tertiary relict plants (Tang et al., 2017). However, to effectively improve the predictive power of SDM projections, we recommend incorporating diverse ecological processes, such as morphology and dispersal strategies, into the future projections.

ACK N OWLED G M ENTS
This work was supported by the funds of China West Normal University (17E067, 17E068, 416446, and 416447).

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Datasets used in this study are available online from the Dryad Digital Repository: https://doi.org/10.5061/dryad.cnp5h qc56.