Response of distribution patterns of two closely related species in Taxus genus to climate change since last inter‐glacial

Abstract Climate change affects the species spatio‐temporal distribution deeply. However, how climate affects the spatio‐temporal distribution pattern of related species on the large scale remains largely unclear. Here, we selected two closely related species in Taxus genus Taxus chinensis and Taxus mairei to explore their distribution pattern. Four environmental variables were employed to simulate the distribution patterns using the optimized Maxent model. The results showed that the highly suitable area of T. chinensis and T. mairei in current period was 1.616 × 105 km2 and 3.093 × 105 km2, respectively. The distribution area of T. chinensis was smaller than that of T. mairei in different periods. Comparison of different periods shown that the distribution area of the two species was almost in stasis from LIG to the future periods. Temperature and precipitation were the main climate factors that determined the potential distribution of the two species. The centroids of T. chinensis and T. mairei were in Sichuan and Hunan provinces in current period, respectively. In the future, the centroid migration direction of the two species would shift towards northeast. Our results revealed that the average elevation distribution of T. chinensis was higher than that of T. mairei. This study sheds new insights into the habitat preference and limiting environment factors of the two related species and provides a valuable reference for the conservation of these two threatened species.

. These migrations will produce new species combinations and species interactions (Guo et al., 2021;Liu et al., 2014;Wang et al., 2020). Meanwhile, migrations will pose a threat of local extinction for many species and/or accelerate the reproduction of some species (Dullinger et al., 2012;Elsen & Tingley, 2015;Faleiro et al., 2018;Hulme, 2017;Wiens, 2016;Zhang et al., 2019). Therefore, it is crucial to understand how climate change alters the distribution of species. Additionally, ample evidence has shown that land-use changes (Fischer et al., 2008;Guo et al., 2019;Li et al., 2020;Madella et al., 2021;Ru, 2006), topography alterations (Elsen et al., 2020;Keppel et al., 2017;Wang et al., 2019), and human behaviors (Gallardo et al., 2015;Pecl et al., 2017) also play a non-negligible role in regulating the effect of environmental changes on species distribution. It is essential to assess the species' habitat by adopting integrated variables and take effective measures to protect the ecological systems.
Species distribution models (SDMs) aim at predicting habitat suitability by integrating species distribution data and environmental data. SDMs have extensively been used to hind-cast and predict future species distribution range change, to assess the impact of species invasion, to reveal niche conservatism, to provide guidance for species reintroduction site selection and conservation strategy formulation (Elith & Leathwick, 2009;Lenoir et al., 2009;Thuiller et al., 2019). Various approaches have been developed to construct SDMs such as Bioclim (Booth et al., 2014), Generalized Linear Models (Guisan et al., 2002), Random Forest (Breiman, 2001), and Maximum Entropy (MaxEnt) (Phillips et al., 2006). Here, we adopted the top-performing Maximum Entropy (MaxEnt) approach which has been widely used in SDMs due to its simple clear graphical interface, high prediction accuracy, and its easy-to-understand output (Elith et al., 2011;Lissovsky & Dudov, 2021;Phillips et al., 2006;Phillips & Dudík, 2008). However, latest research has shown that MaxEnt is prone to over-fitting, resulting in low model transfer ability, which seriously affects its application to various fields such as the invasion biology and phylogeography (Jiménez-Valverde et al., 2011;Syfert et al., 2013;Zhu & Qiao, 2016). Model complexity of MaxEnt is mainly affected by the 4 parameters, namely, background data, feature class (FC), regularization multiplier (RM), and sampling bias (Merow et al., 2013). Recently, Muscarella et al. (2014) have developed the ENMeval package to perform automated tuning and evaluations of species distribution models. This species-specific tuning in Maxent settings can avoid over-fitting in niche models and improve predictive ability.
Taxus is the largest and most widely distributed species in Taxaceae (Fu et al., 1999). In the 1980s, Taxus attracted great attention, meanwhile suffering huge damage since the Taxol extracted from it was found to be one of the most popular natural anticancer materials . Due to over-exploitation and anthropogenic disturbances, Taxus population diminished sharply and become fragmented. At present, it is at a high risk of extinction (Liu et al., 2011Yu et al., 2014). Moreover, its biological properties such as low pollination rate, long seed dormancy, and weak competitive ability of seedlings also cause its current endangered status Liu, Feng, et al., 2019;.
Internationally, the three Taxus species (Taxus wallichiana, Taxus contorta, and Taxus chinensis) are listed as endangered species (EN, Thomas et al., 2013;Thomas, 2011;Thomas & Farjon, 2011), and Taxus mairei as Vulnerable (VU, Yang et al., 2013), and Taxus cuspidata as the least concern (LC, Katsuki & Luscombe, 2013) 1999). It is worth mentioning that the two closely related Taxus species Taxus chinensis and Taxus mairei are mainly distributed in the Sino-Himalayan forest floristic subkingdom and the Sino-Japanese forest floristic subkingdom (Wu & Wu, 1996). Taxus chinensis is endemic to China, and it is mainly distributed in the mountains around Sichuan Basin. Taxus mairei has a wide distribution in the south of the Yangtze River in China and other countries of South Asia and Southeast Asia. Generally speaking, Taxus mairei is usually at lower elevation than Taxus chinensis (Fu et al., 1999). However, the morphological characteristics of the two species are almost the same, and thus it is difficult to distinguish them (Farjon, 2017). As an ancient and long-lived tree species, the two species experience Quaternary glaciation Möller et al., 2020). However, how the two species respond to elevation change and climatic oscillation in the Quaternary remains to be further explored .
With the increasing energy consumption, global warming is one of the main challenges in the 21st century (Durán-Martín et al., 2019).
The distribution change of these two species in response to global warming is largely unknown. Hence, the mappings of suitable habitats and predictions of the impacts of climate change are vital for habitat protection and the sustainable development of these two species.
In this study, we integrated optimized species distribution models (SDMs) and geographical information system (GIS) software to analyze the two species distribution pattern in response to climate change with the aims to (1) uncover the dominant environmental factors in their niche differentiation; (2) determine elevation differences of the two species since last inter-glacial (LIG) periods; (3) reveal the conservation implications for the species. Overall, this study will deepen our understanding of their evolutionary history and provide some useful guidelines for the conservation of these two threatened species.

| Species occurrence data
Organism photographs in the fruiting stage (August-December) of T. chinensis and T. mairei were shown in Figure 1. Species occurrence data of T. chinensis and T. mairei were collected from the fieldwork, previous studies (Liu et al., 2011, the Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn/) and the Plant Photo Bank of China (PPBC, http://www.plant photo.cn/). (Number of occurrence data for each species from each dataset was shown in Table S1.) Due to Taxus chinensis and Taxus mairei were mainly distributed in China, and the aim of the study is to infer the potential distribution area of the two species in China; thus, we do not consider the Global Biodiversity Information Facility (GBIF) database in the current study.
Then, the data with obvious geographical coordinates errors were removed by the ArcGIS 10.4. In addition, the duplicate data were removed to ensure only one record in the 2.5′ × 2.5′ grid by the "spThin" package. Finally, a total of 63 sampling points for T. chinensis, and 140 sampling points for T. mairei were retained ( Figure S1; Table S1).

| Predictor variables
The relevance and completeness of variables are key components for constructing SDMs (Elith & Leathwick, 2009;Guo et al., 2017;Zimmermann et al., 2010) Figure S2). Three topographical variables, including the elevation obtained from the WorldClim database (https://world clim. org), and the slope and aspect obtained from Digital Elevation Model using the 3D analyst tools in the software ArcGIS 10.4. Five soil variables were downloaded from the Harmonized World Soil Database (HWSD, https://www.fao.org/soils -porta l/data-hub/soil-maps-anddatab ases/harmo nized -world -soil-datab ase-v12/en/) based on previous studies Li, Zhang, & Griffith, 2021;Li, Zhang, Zhu, et al., 2021;Ru, 2006), and the correlation coefficient between the five soil variables was less than 0.7 ( Figure S3). Human interference index (HII) was downloaded from Socioeconomic Data and Applications Center (SEDAC, http://sedac.ciesin.colum bia.edu) (Gallardo et al., 2015;Madella et al., 2021). Ultimately, 15 ecological variables were chosen for each species for further analysis (Table 1). CCSM model is a coupled climate model for simulating the earth's climate system, whereas RCP4.5 provides a platform for climate models to explore the climate system response to stabilizing the anthropogenic components of radiative forcing (Thomson et al., 2011).

Furthermore, bioclimatic variables of
However, there are no available data on the other three types of environmental variables (topographical variable, soil variables, and human interference index) for the past and future periods corresponding to the same periods. Since the two species are mainly distributed in the mountain areas, they are relatively less influenced by topograghy, soil, and human behaviors. Therefore, these three types of environment variables were assumed to be constant, as reported in previous studies (Evans et al., 2020;Lv et al., 2021;Zhang et al., 2019).

| Modeling procedures
We conservatively defined the background range as an area surrounding our occurrence location by buffering a bounding box. We set the buffer distance as 7 degrees (about 779 km) based on the dispersion distance of the Taxus pollens and seeds (Brown et al., 2017;Li et al., 2015;Li, Zhang, & Griffith, 2021;Li, Zhang, Zhu, et al., 2021).
According to the relevant parameter ahead, Maxent v3.4.1 was used to investigate the effects of past and present climatic conditions on T. chinensis and T. mairei (Phillips et al., 2006). The 75% of species occurrence data and the remaining 25% were used as training data and testing data for model validation, respectively. The 10 replications and 5000 bootstrap iterations were set, and other parameters had the default settings. Model performance was evaluated by the area under the receiver operating characteristics curve (AUC). The AUC ranged from 0 to 1. The model with AUC value of more than 0.9 was considered as excellent (Araújo et al., 2005;Elith et al., 2006).

| Distribution area change and centroid transition
Based on the results of the species distribution modeling, SDMtoolbox (Brown et al., 2017) in ArcGIS v10.4 was used to evaluate the changes in the area during different periods and the centroid shift for the two Taxus species. We cross-checked the changes of the highly suitable area to identify regions as (i) expansion, (ii) unchanged, and (iii) contraction relative to the previous periods . Centroid shift concentrated the species distribution on an independent central point and created a vector file to depict the magnitude and direction of changes over time (Cong et al., 2020;Hu et al., 2015).

| Model accuracy and contributions of predictor variables
Species distribution modeling was constructed for each species to predict its geographical distributions at present, in the past and future. Based on the results of ENMeval, the optimal combination (Delta-AICc = 0) ( Figure S4) of RM/FC for T. chinensis and T. mairei was 0.5/LQ and 1/LQH, respectively. These parameters combination could avoid model's over-fitting and improve its prediction ability. Our data showed that the area under the receiver operating characteristics curve (AUC) obtained from all the models was larger than 0.9 (AUC > 0.9), indicating the robustness and reliability of predictions of our models ( Figure S5). were presented in Figure S6. The kernel density clearly showed that the two Taxus species had different temperature preferences, and the suitable temperature was about −10 to 15°C for T. chinensis and −5 to 20°C for T. mairei with an average temperature of 1.79 and 5.60°C, respectively (Figure 2; Figure S6).

| Distribution area during different periods
The potential distributions of the two species during different climatic periods were presented in Figure 3. The suitable area

| Centroid migration
In different periods, the centroid of T. chinensis was distributed in the adjacent areas of Sichuan province and Chongqing city, while T. mairei was distributed in Hunan province ( Figure 5) Figure 5).
The migration distance of each species between the two adjacent periods was shown in Table 3. For T. chinensis, the largest migration distance was 33.53 km from the LIG to the LGM periods. For

T. mairei, the largest migration distance was from the LGM to the
Holocene periods with the migration distance of 80.27 km.

| Elevational differences
To reveal the elevation difference between the two species, we cal-  indicate that Taxus is shade-tolerant species, and that it prefers to grow along the river (Liu et al., 2013;Song, 2013). In this study, the mean temperature of the coldest quarter (bio11) was found to be the most important and most contributing factor for T. chinensis. Thus the temperature was the main factor influencing its spatial distribution. For T. mairei, the most contribution and permutation importance factor was bio18 (precipitation in warmest quarter). Wang et al. (2019) have shown that annual precipitation (bio12) and topographical variables have a strong effect on the distribution of T. chinensis and T. mairei. Li et al. (2022), Liu, Feng, et al. (2019) and  have also supported that precipitation is the most important climatic factor that restricts the habitat distribution of the T. mairei. Poudel et al. (2012) have reported that great differences in rainfall between winter (low) and summer (high) are the determining factor responsible for the distribution of T. mairei in the east of the Himalayas in Nepal.
Overall, precipitation is the dominant factor determining the distribution of the T. mairei.

F I G U R E 5 Centroid migration routes under different climate periods for Taxus chinensis and Taxus mairei
It is worth mentioning that rooting condition is the third permutation importance factor (15.6%) for T. mairei (Table 1). Ru (2006) has shown that T. mairei prefers living in an environment with moist fertile soil and good water permeability. Owing to T. mairei lives in relatively low elevation areas with ample environment moisture; therefore, good water permeability is conducive to the growth of T. mairei.

| Changes in species distribution area
The distribution area of T. chinensis and T. mairei were almost in stasis from the LIG to the future (RCP4.5_2070). This can be understood from the following two points. First, biological traits such as limited dispersal capacity, long generation time and low rate of seed germination of Taxus (Keppel et al., 2017;Li, Zhang, & Griffith, 2021;Li, Zhang, Zhu, et al., 2021;Ru, 2006;Wang et al., 2018Wang et al., , 2019. Second, the main distribution region of the natural populations of the two species were in the mountains such as the Qinling, Nanling, and Wuyi Mountains in China (Fu et al., 1999;Li et al., 2022). T. chinensis is mainly distributed around the Sichuan Basin, while T. mairei occupies most of the southern regions of the Qinling-Daba Mountains in China. These mountains not only provide a relatively stable habitat for species but also act as the refuge (Jiang et al., 2019;Keppel et al., 2017;Ye et al., 2017;Zhao et al., 2019). This phenomenon has also been reported on T. mairei, T. wallichiana , Tsoongiodendron odorum (Hu et al., 2020), and Houttuynia cordata  and Eucalyptus grandis (Ouyang et al., 2022). Species distribution models (SDMs) are based on the niche conservatism hypothesis and niche conservatism is more prevalent than niche differentiation (Alexander, 2013;Chivers et al., 2017). We have adjusted parameters many times to test the changing tendency of the species distribution area between two adjacent periods, the end result, however, was about the same. Thus, the models was not the main reasons that lead to the overall stasis in the species distribution area of the two species.
Notably, the suitable area of the two species will shrink from current period to future period because the plant growth, development and reproduction are vulnerable to the effects of increasing global temperature (Liu, Feng, et al., 2019;. Our result is consistent with the previous findings of study of T. mairei (Li, Zhang, & Griffith, 2021;Li, Zhang, Zhu, et al., 2021). Meanwhile, the suitable area of Taxus cuspidata and Taxus wallichiana will be reduced with the rising temperature in the future (the 2050 and 2070) Su et al., 2018). This shrinking tendency is also observed in other species such as Quercus lamellosa (Guo et al., 2021) and Polyporus umbellatus .
Species may change its latitude or elevation in response to climate changes (Davis & Shaw, 2001). Previous studies have shown that species will move northward and upwards with the increasing temperature, such as Quercus lamellosa (Guo et al., 2021), Cyananthus , Quercus kerrii (Jiang et al., 2018), wild soybean (He et al., 2016), and T. wallichiana . Likewise, this study also found an upward and northward shift trend for T. mairei, which is consistent with the reports by Li, Zhang, and Griffith (2021), Li, Zhang, Zhu, et al. (2021) and Poudel et al. (2012). Instead, T. chinensis shifted downward and northward in China. This may be due to the fact that survival pressure from the higher elevation is not conducive to the growth and reproduction of T. chinensis, and that the special topography around the Sichuan Basin may provide a route for T. chinensis to migrate northward. Furthermore, Liang et al. (2018) modeled 151 representative plants in the Hengduan Mountains and its adjacent areas in China, and found that the mountain plants shifted upward with the increasing temperature, but the shift was not only northward but also westward or in other directions.

| Elevational differences of two Taxus species
Elevational differences were observed between the two closely related Pinus species Pinus massoniana and Pinus hwangshanensis. The reason for such elevational differences lie in that species specificity and climatic divergence selection of the candidate genes play a key role in the ecological divergence of these two species (Li et al., 2010;Zhou et al., 2014). Theoretically, closely related species are expected to show more similarity as a consequence of shared climate selection, habitat, and evolutionary history (Miller et al., 2019;Nürk et al., 2015). Recent climatic selection may be species-specific since forest trees typically have the highest adaptation in their own environment (local adaptation), and different species typically occupy different climatic niches (ecological niche differences) (Hua & Wiens, 2013;Savolainen et al., 2007). Our species distribution models (SDMs) results indicated that T. chinensis is mainly located at the elevation above 1500 m, and T. mairei tends to occur at the elevation of around 1000 m since LIG. Our results are in accordance with the description of Flora of China (Fu et al., 1999). Species distribution models results shown that there is no large-scale population migration for the two closely related species. Temperature and precipitation are the main factors determining species distribution regions.
Therefore, this study shown that climatic selection and long-term adaption to a given environment might the main factors that influencing the two species divergence along the elevation. Furthermore, we have observed a hybrid zone (T. chinensis and T. mairei) in the intermediate transition zone between high altitude and low altitude region (unpublished data) in the nature reserve, and this hybrid zone provides good materials for us to explore the dynamic history of the two closely related species at the molecular level.

| Implications for conservation and management
From current to future periods (RCP4.5_2070), the suitable distribution area of T. chinensis and T. mairei will shrink. First, upward migration of species along elevational gradients will lead to range contraction for many species since the total area available at a given altitude usually decreases with increasing elevation in mountains (Parmesan, 2006;Wilson et al., 2005). Second, from the view of biological characteristics, Taxus prefers a shady and humid environment (Wu & Wen, 2017). However, with the rising temperature, the climatic conditions such as strong radiation, drought, wind, and other adverse climatic conditions will be more severe in high-altitude areas (Ouyang et al., 2022;Solomon et al., 2007;Yin et al., 2020). Thus, these adverse conditions will pose relatively more physiological constraints on T. chinensis, thus resulting in range reduction.
However, the reason for the range contraction of T. mairei may be different from that of T. chinensis. Rooting conditions is the third most important factor that affects the distribution of T.
mairei. Besides, the increasing studies have shown that temperature and precipitation are the main factors that affect rooting (Fang et al., 2017;Reich et al., 2018). In general, global warming is expected to cause changes in distribution, intensity, and frequency of precipitation (Myhre et al., 2018). Inappropriate hydrothermal conditions are not conducive to the rooting of seeds, eventually resulting in the population number decline. Moreover, the anthropogenic disturbance is stronger in the low-elevation areas than in high-elevation areas, and hence it may also lead to the contraction in species' range.
Taken together, we can establish germplasm resource nurseries to cultivate the seeds from different provenances, especially for the T. chinensis distributed in high-elevation areas. Furthermore, considering the influence of humans, in-situ protection should be enforced for the samples that are easily accessible. Theoretically, evolution can also drive species distribution range shifts in the absence of environmental change (Holt, 2003;Parmesan, 2006), such as the inter-species interactions, hybridization, and introgression, and thus common garden experiments should be undertaken to investigate potential local adaptations and facilitate the development of future genetic studies.
It should be noted that although predictions based on species distribution models (SDMs) effectively uncover the dynamic population history of the two species, there are also some limitations in current study. First, soil variables and human interference index are assumed to be constant in the ancient climate and future climate. However, species could be threatened, or even possibly become extinct in the case of a dramatic increase in human population and land use, thus resulting in habitat loss and fragmentation (Giam et al., 2010). Thus, appropriate protection measures should be taken when anthropogenic disturbances drastically increase. Additionally, the effects of solar radiation on the SDMs results are not taken into consideration in the current study (Ouyang et al., 2022). Because the two species are shade-preferring species and were always the associated tree species, solar radiation may not have a direct impact on them (Liu et al., 2013;Su et al., 2018). However, we suggest that solar radiation should be taken into account in the future research, especially for the heliophilous species and those species distributed in the regions with high solar radiation such as Qinghai-Tibet Plateau. Finally, for woody plants with long generations, the change of climate from suitable to unsuitable does not mean the disappearance of species distribution in a specified area, instead, unsuitable climate may involve more environmental stress that species may need to suffer. Thus, we should consider the effect of climate on species distribution pattern when evaluating the endangered category of the species.

| CON CLUS ION
Our findings enhance our understanding of the past and present

ACK N OWLED G M ENTS
We would like to thank the editors and two anonymous reviewers for their insightful comments and efforts in improving the clarity of this manuscript. This study was supported by the National Natural

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that supports the findings of this study are available in the supplementary material of this article.