Integration of multiple climate models to predict range shifts and identify management priorities of the endangered Taxus wallichiana in the Himalaya–Hengduan Mountain region

Taxus wallichiana Zucc. (Himalayan yew) is subject to international and national conservation measures because of its over-exploitation and decline over the last 30 years. Predicting the impact of climate change on T. wallichiana’s distribution might help protect the wild populations and plan effective ex situ measures or cultivate successfully. Considering the complexity of climates and the uncertainty inherent in climate modeling for mountainous regions, we integrated three Representative Concentration Pathways (RCPs) (i.e., RCP2.6, RCP4.5, RCP8.5) based on datasets from 14 Global Climate Models of Coupled Model Intercomparison Project, Phase 5 to: (1) predict the potential distribution of T. wallichiana under recent past (1960–1990, hereafter “current”) and future (2050s and 2070s) scenarios with the species distribution model MaxEnt.; and (2) quantify the climatic factors influencing the distribution. In respond to the future warming climate scenarios, (1) highly suitable areas for T. wallichiana would decrease by 31–55% at a rate of 3–7%/10a; (2) moderately suitable areas would decrease by 20–30% at a rate of 2–4%/10a; (3) the average elevation of potential suitable sites for T. wallichiana would shift up-slope by 390 m (15%) to 948 m (36%) at a rate of 42–100 m/10a. Average annual temperature (contribution rate ca. 61%), isothermality and temperature seasonality (20%), and annual precipitation (17%) were the main climatic variables affecting T. wallichiana habitats. Prior protected areas and suitable planting areas must be delimited from the future potential distributions, especially the intersection areas at different suitability levels. It is helpful to promote the sustainable utilization of this precious resource by prohibiting exploitation and ex situ restoring wild resources, as well as artificially planting considering climate suitability.


Introduction
The Himalaya-Hengduan Mountain (HHM) region is located in the world's Third Pole. It is a biodiversity hotspot and home to Taxus wallichiana Zucc. (Himalayan yew), a relict tree species of the Quaternary glaciation. It is well-known because of taxol, one of the most successful anticancer drugs derived from natural sources. Recent research on the taxonomy of Taxus species  laid the foundation for our study of species distribution and regional conservation. T. wallichiana is the most widely distributed yew in the HHM region. As a relict species from the Quaternary glaciation, it has a history of 2.5 million years on earth and it is known as a living fossil of the plant kingdom. The genus, Taxus, is known for its medicinal value. Taxol was first extracted from T. brevifolia in 1971 (Wani et al. 1971). Early experiments demonstrated that taxol could completely inhibit the exponential growth of cancer cells at low concentrations (Schiff et al. 1979). Subsequently, the effectiveness of taxol has been documented for the treatment of various cancers, inflammatory conditions, and Acquired Immune Deficiency Syndrome (Yan et al. 2013;Yang et al. 2017). As one of the most successful anticancer drugs derived from natural sources (Yang et al. 2017), taxol has a huge and expanding market demand (Miao et al. 2015). The bark and leaves of T. wallichiana are now used to produce taxol, and this is the reason for its exploitation (Thomas and Farjon 2011;Uniyal 2013). The population in the HHM region has declined significantly ([ 50% in China,B 90% in Nepal and India) over the last 25-30 years (a single generation) (Thomas and Farjon 2011).
In 1995, T. wallichiana. was assessed as an endangered species by the International Union for Conservation of Nature (IUCN 2018) and trade in the species was regulated by the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES 2007). It has been listed as a national first-class protected plant in China (State Forestry Administration of China 1999), and it is on the Negative List of Exports of India (Sajwan and Prakash 2007). Some protected areas are established, e.g., Three Parallel Rivers of Yunnan Protected Areas in China, Sagarmatha National Park in Nepal, and BiDoup-NuiBa National Park in Vietnam. However, in order to restore wild populations or cultivate T. wallichiana successfully, the impact of global environmental change (especially climate warming) on T. wallichiana must be considered due to its relatively restricted and scattered geographical distribution (Su et al. 2005; Thomas and Farjon 2011;Poudel et al. 2014), and its weak ecological adaptability (e.g., poor seed regeneration, long pre-reproductive phase in nature) (Paul et al. 2013;Uniyal 2013). Global climate warming has accelerated (IPCC 2013), especially in areas supporting habitats suitable for T wallichiana in HHM (Fig. S1). Globally, this is affecting the distribution and abundance of species (Hughes 2000;Root et al. 2003;Pacifici et al. 2015;Asner et al. 2016), and is leading to increasing challenges for the conservation of biodiversity (Myers et al. 2000;Pereira et al. 2010;Garcia et al. 2014). Thus, it is important to accurately predict the impact of climate change on species distribution to as a basis for the planning of nature reserves and to protect priority areas for endangered species (Pyke et al. 2005).
Species distribution models (SDMs) which can be drived by Global Climate Model (GCMs) data are powerful tools for predicting species distribution, genetic and evolutionary research, zoning management and protection in the face of climate change. For example, Liu et al. (2013) used molecular biology theories and the MaxEnt model to analyze cryptic speciation and predict changes in the potential distribution of two T. wallichiana lineages [i.e., EH lineage (East Himalaya to the Yunnan Plateau region) and HM lineage (South Hengduan Mountains region)] from the Last Interglacial (LIG, ca. 120 ka) and Last Glacial Maximum (LGM, ca. 21 ka) to the present (ca. 1950-2000) in the HHM region. In a warming climate, Poudel et al. (2014) used similar methods to examine genetic diversity and population differentiation, and to predict current and future distributions of three Taxus species (T. contorta, T. mairei and T. wallichiana) in the Central Himalayas. Gajurel et al. (2014) input climatic variables generated by the HadGEM2-ES climate model into the MaxEnt model to predict the current and future potential distribution of T. wallichiana in the Nepal Himalayas. In general, previous predictions of the potential distribution of Taxus species were mainly based on climate scenario datasets generated by single or several climate models. However, due to the complexity of earth's climate systems, as well as the high uncertainty surrounding climate modeling for mountainous regions (due to the paucity of weather stations) (Hijmans et al. 2005), we set out to integrate multiple scenario datasets from multiple climate models to predict the potential distribution of Taxus species.
Our objectives were to: (1) investigate changes in the potential distribution of T. wallichiana between current  and future (2050s and 2070s) climate scenarios by integrating three Representative Concentration Pathway (RCP) (i.e., RCP2.6, RCP4.5, RCP8.5) datasets from 14 GCMs of CMIP5 (Coupled Model Intercomparison Project, Phase 5); (2) to quantify the major climatic factors influencing the distribution; and (3) to propose management interventions based on the resulting predictions.

Study area
We set our study area boundaries as 17°50 0 41 00 N to 33°43 0 40 00 N, and 75°55 0 52 00 E to 107°38 0 15 00 E based on the current distribution region of T. wallichiana (Liu et al. 2012(Liu et al. , 2013Poudel et al. 2014), the species distribution database of the Chinese Virtual Herbarium (CVH), and the Global Biodiversity Information Facility (GBIF). This covered an area of 5,478,120 km 2 , ranging in elevation from 0 to 8000 m above sea level (Fig. 1). Globally, T. wallichiana is mainly distributed in the HHM region. The general circulation pattern in this region is characterized by the influence of Indian and East Asian monsoons in summer, and Westerlies in winter Yao et al. 2012). The habitats currently suitable for T. wallichiana in the HHM region are expected to experience intense warming and slight humidification in the future scenarios (Fig. S1).

Sample data of T. wallichiana
Records for the geographical distribution of T. wallichiana were collected from the literature based on the most recent knowledge of the taxonomy of Taxus species in the study area (Liu et al. 2012(Liu et al. , 2013Poudel et al. 2014). Sixty-nine records were found with species name, longitude, and latitude (Table S1). These data were used to train and validate the current and future distribution simulations of T. wallichiana.

Climate data
Global climatic data were downloaded from WorldClim (www.worldclim.org) (Hijmans et al. 2005). Among them, the current climate data  were produced based on the observed data in the weather stations using the thinplate smoothing spline algorithm implemented in ANUS-PLIN (Hutchinson 2004). The cross-validation error is between 0 and 1°C for temperature in most area and less than 10 mm/month for precipitation in the vast majority of places, respectively, despite the high uncertainty in mountainous areas (see details in Hijmans et al. 2005). The future climate data (outputs of GCMs) were downscaled (interpolated) and calibrated from its original resolution (about two to three geographical degree) to 30 arc-second (0.93 km 9 0.93 km = 0.86 km 2 at the equator) using current climate data as the baseline. We downloaded one geographical variable and 19 climatic variables (including the current conditions and predictions for future scenarios (outputs from 14 GCMs)) at a spatial resolution of 30 arcseconds (Table S2). All global spatial data adopted a WGS84 coordinate system and a cylindrical orthographic projection. The climatic data from 1960 to 1990 were used to generate the current potential distribution of T. Integration of multiple climate models to predict range shifts and identify management… wallichiana, and the climatic data for three scenarios (RCP2.6, RCP4.5, RCP8.5) in the 2050s (2041-2060) and 2070s (2061-2080) outputted from 14 GCMs (Table S3) were used to generate future potential distribution. A total of 20 Climatic and Geographical Variables (20 CGVs) for the HHM region were clipped and extracted, and their spatial resolution, coordinate system and projection were the same as the original spatial data.
To remove multicollinearity among variables, a Pearson's correlation (r) matrix (Table S4) was computed with SDMtoolbox 2.0 software (Brown et al. 2017) for all CGVs. Then, the input variables for distribution simulation were selected according to the rule |r| \ 0.80 Bosso et al. 2017) combined with the following ecological features of T. wallichiana. In Table S4, CGVs BIO2, BIO3, and BIO15 were reserved according to the rule |r| \ 0.80. BIO1 was selected among BIO5, BIO6, BIO8, BIO9, BIO10, BIO11 and elevation (topographic variable) because the annual average temperature has a more general and direct impact on plant survival and growth than extreme temperature and elevation (Theurillat and Guisan 2001). Because the study area has a monsoon climate Yao et al. 2012;Schickhoff et al. 2015) (Fig. 1), BIO4 (temperature seasonality) is a more accurate indicator for plant growth and development than BIO7 (annual temperature range). Therefore, BIO4 was reserved. BIO12 was also selected because precipitation is a useful ecological index of T. wallichiana's survival. BIO14 (the driest monthly precipitation) and BIO19 (the coldest seasonal precipitation) were also selected as they are the stress variables of T. wallichiana under low temperature and dry conditions. As a result, eight CGVs (Table 1) were selected to simulate the current and future potential distribution of T. wallichiana. Figure 2 shows the overall technical flowchart. First, CGVs were selected for inputting to the MaxEnt model after eliminating multicollinearity among all variables. Then, the MaxEnt model was parameterized to simulate the potential distribution of T. wallichiana under current and future climate scenarios. The Receiver Operator Characteristic (ROC) analysis and uncertainty analysis (Spatial Inconsistency Rate and Coefficient of Variation for area) were used to validate the simulation accuracy and assess the simulation differences among the 14 climate models. Finally, the spatial and temporal changes in the potential distribution of T. wallichiana were investigated, and the percent contribution rate and the response range of influencing factors were detected.

Methods
Predicting the potential suitable distribution and the factor contribution rate of T. wallichiana Based on ecological theories and statistical analyses, SDMs simulate species distribution by relating species observations with environmental variables (Guisan and Zimmermann 2000). The MaxEnt model (Elith et al. 2006;Phillips et al. 2006Phillips et al. , 2017 is one SDM which has been widely and successfully applied to the prediction of plant habitats worldwide (Costion et al. 2015;Bocksberger et al. 2016;Booth 2018;Zhang et al. 2018;Dyderski et al. 2018). It is particularly suitable for the prediction of scarce distribution points compared to other presence-only SDMs (Elith et al. 2006;Hernandez et al. 2006;Wisz et al. 2008). Therefore, we used the MaxEnt model (version 3.4.1) to simulate the potential distribution of T. wallichiana.
The MaxEnt model was parameterized according to Bosso et al. (2017). First, the sample points of T. wallichiana and the selected eight CGVs for the current and for one of the future scenarios were input to the MaxEnt model. The cross-validation method was adopted for each simulation (Bosso et al. 2017) and 30 replicates were set.
There was a total of 84 simulations (14 GCMs 9 3 RCPs 9 2 periods). Each simulation simultaneously generated both current and future potential distributions. Because there were 30 replicates for each simulation, the final results could be statistically represented as the mean of 30 replicates to reduce random error. The final average map had a complementary log-log (cloglog) output format with continuous suitability values from 0 (completely unsuitable habitat) to 1 (most suitable habitat) (Phillips et al. 2017). The contribution rate for each factor is automaticall computed with the MaxEnt model. It can derive many features for each predictive variable, each of which is a simple mathematical transformation of the predictive variable (Merow et al. 2013). The model assigns the increase in the gain to the predictive variable(s) on which the feature depends. Finally they are converted into percentages to represent the contributions rates (Phillips 2005).

Verifying the simulation performances and differences
The performance of simulation and the significance of each environmental variable were evaluated by the Area Under Curve (AUC) of ROC with the Jackknife test. Different AUC values indicate different suitable levels: poor (AUC \ 0.80), fair (0.80 B AUC \ 0.90), good (0.90 B AUC \ 0.95), or excellent (0.95 B AUC B 1.00) (Thuiller et al. 2005).
The Spatial Inconsistent Rate (SIR) was used to evaluate the spatial distribution differences among simulations. First, the common suitable area among the 14 GCMs outputs under the future scenarios was computed (the common suitable area indicates the suitable area intersection of 14 GCMs outputs under the RCPs scenarios.). Then, the SIR was calculated according to formula (1): where i indicates a single simulation, SIR i the spatial inconsistent rate for the i single simulation, SA i the predicted suitable area for the i single simulation, and CA the common suitable area.
The Coefficient of Variation (CV) for area was used to evaluate the differences of suitable areas at different levels among future scenarios of 14 GCMs outputs. The CV of area was calculated according to formula (2): where i indicates a future scenario, j a suitable level, u ij the average area among 14 simulations for the i future scenario at the j suitable level, and , ij the standard error.

Classifying the climate suitability level
The climate suitability level was classified according to the statement of ''possibility'' in IPCC AR5, which matched the probability results of the MaxEnt model output. The study area was first classified as either a non-suitable area (0 B P \ 0.05) or a suitable area (0.05 B P B 1.00), and the suitable area was further classified as a marginally suitable (0.05 B P \ 0.33), moderately suitable (0.33 B P \ 0.66), and highly suitable (0.66 B P B 1.00).

Computing the areas for different levels of suitability
The predicted results from Maxent had the same spatial resolution, geographic coordinate system and projection as the original data. However, because all the spatial data adopted a cylindrical orthographic projection, the areas for each pixel at different latitudes were not equal. So we first set the area of a pixel at the equator (0.93 km 9 0.93 km = 0.86 km 2 ) as the benchmark, and the area of a pixel at other latitudes was then weighted by the square root of the cosine of its corresponding latitude (Wallace et al. 2000;Marshall et al. 2012).  Integration of multiple climate models to predict range shifts and identify management…

Testing statistical significance
We used paired t test to compare mean simulated areas by suitability level and by elevation between current and future scenarios (Bosso et al. 2017;Liu et al. 2018). The threshold p \ 0.05 was selected to indicate significance of differences between means, and p \ 0.01 indicated extremely significant difference.

Simulation performances
Among the 84 simulations, the highest and lowest AUC training mean ± standard values were 0.989 ± 0.000, and 0.986 ± 0.000, respectively, and the highest and lowest AUC test values were 0.987 ± 0.012, and 0.984 ± 0.015, respectively (Table S5). All the simulations showed excellent performance. Fifty-seven of the 69 actual distribution points were located within the moderately and highly suitable areas of the current predictions (Table S6).
Among the 84 simulations, 57 showed a reduction of the suitable area under future scenarios compared with the current mean condition (Table S7, Table S9). Significant reduction in suitable area was predicted under the RCP2.5 and RCP4.5 scenarios for the 2050s. The relative reduction ratios (area of reduction to current area, or ratio of reduction) were 11% (p \ 0.05) and 9% (p \ 0.05) ( Table 2), respectively. The relative rates of reduction (the linear regression slope of the relative reduction ratio and year) were 2%/10a and 1%/10a, respectively (Table 3). No significant reduction in the area of suitable habitat was predicted in other future scenarios. In contrast, a significant increase was predicted under the RCP8.5 scenario for the 2070s, and the relative increase ratio was 21% (p \ 0.05) ( Table 2). The relative rate of increase was 2%/10a ( Table 3).
The upslope shift of the moderately suitable and marginally suitable areas were predicted to be greater than the predicted average elevation change of the highly suitable area (see details in Tables 4, 5, and Fig. 5).

Main climatic factors affecting distribution of potentially suitable habitats
The potential distribution of T. wallichiana was mainly influenced by the annual mean temperature (contribution rate ca. 61 ± 3%), annual precipitation (ca. 17 ± 2%), isothermality (ca. 10 ± 1%), and temperature seasonality (ca. 10 ± 2%), which indicated that annual mean temperature and its seasonal variation, as well as annual precipitation, were the main climate variables affecting the habitat of T. wallichiana (Table 6, Table S11). Seasonality of (ca. 1 ± 0%), the precipitation of the coldest quarter (ca. 1 ± 0%), mean diurnal range of temperature (ca. 0 ± 0%), Integration of multiple climate models to predict range shifts and identify management…   Integration of multiple climate models to predict range shifts and identify management… and precipitation of the driest month (ca. 0 ± 0%) were other important factors (Table 6, Table S11).

Discussion
Uncertainties in the simulated potential distribution of T. wallichiana Our models predicted substantial differences in the potential area and distribution range of T. wallichiana under future climate scenarios (Table 7). Among the differences, the CV for the simulated future suitable areas among 14 GCMs was large (12-29%), especially for the moderately suitable areas (13-34%) and the highly suitable areas (38-64%). Moreover, the SIR for the simulated future suitable areas among 14 GCMs reached 91 ± 26% to 259 ± 104%. These results showed the considerable uncertainty when a single GCM dataset was used, and highlighted the need for integrating multiple scenario datasets from multiple climate models. Some previous research reported considerable uncertainty when using  GCMs datasets to predict the potential future distribution of species. Bosso et al. (2017) used six GCMs to predict the current and future potential suitable area of the Italian tree pathogen, Diplodia sapinea, and the simulations with different climate models showed large differences. Ikegami and Jenkins (2018) used 14 GCMs to predict the current and future risks of pine wilt disease for 21 pine species based on the MaxEnt model, and the predicted risk area change rate (max-min) was greater than 20% for nearly half of the species among 14 GCMs. Cuyckens et al. (2016) used three GCMs (MIROC-ESM, HAdGEM2-ES, IPSL-CM5A-LR) to predict future changes in suitable habitat for Polylepis tarapacana forest in the semi-arid mountainous region of South America. Under the RCP6.0 scenario in 2080, the predicted area with IPSL-CM5A-LR decreased by 42%, while the predicted area with MIROC-ESM decreased by 83%. Wenger et al. (2013) reported that the greatest contribution to uncertainty in predictions of area of suitable habitat was climate uncertainty, followed by the parameter uncertainty and model uncertainty. The variance of future climate (GCMs and SRES) is the largest contributor to the overall uncertainty of the impact of climate change (Bagchi et al. 2013). Integrated forecasting is a promising approach to deal with this uncertainty (Araújo and New 2007). It uses the multi-model inference or multimodel averaging value. In general, the predictions that use more environmental factors (assuming that multicollinearity is eliminated) can be closer to an actual species ecological niche when considering the uncertainties of environmental factors in simulating potential distribution. Our simulations were mainly based on climatic variables because they are the dominant factors affecting species distribution and change on a regional scale ([ 200 km) (Pearson and Dawson 2003). However, on a local scale (\ 10 km), the growth of a plant species is affected by other environmental factors, such as soil characteristics (Pearson and Dawson 2003), biological traits , interspecific competition, and disturbance (Liang et al. 2018). For example, Liang et al. (2016) pointed out that the upslope shift rate of a tree line was largely determined by interspecific relationships (e.g., mutual benefits, interspecific competition and their intensities) with future climate warming. Based on 10-year micro-meteorological observations and a 4-year seedling transplantation experiment of Smith-fir (Abies georgei var. smithii) in the Sergyemla Mountains, southeast Tibet, it was found that climate warming advanced the growth season, which caused early-season freezing, low temperature photoinhibition, and high sunlight significantly reduced seedling survival above the tree line . Thus, the upslope shift of the tree line will be influenced by multiple factors and become lower than the predicted value based on climate variables. The main distribution area of T. wallichiana is located in the high vegetation cover area of southwest China which is an important biodiversity hotspot (Myers et al. 2000). Climate oscillations here have a profound effect on population size (i.e., change, extinction, or promoting local adaptation) (Davis and Shaw 2001), as do the ecological relationships between plant species and other trophic levels (Parmesan and Hanley 2015). Human over-exploitation is an important factor affecting the population of wild T. wallichiana, so the current distribution of the species might not coincide with the most suitable habitat. However, loggers are unlikely to deforest all T. wallichiana in the short term (there are no records of deforestation), and the sampling points from the literature were mainly collected in 2010, when the actual distribution of T. wallichiana should have been more widespread but was constricted by deforestation. So the simulated current distribution of the species still represented the most suitable and available habitat, and simulated current distribution tended to reflect the potential Integration of multiple climate models to predict range shifts and identify management… niche of the species rather than the real one, since the biological factors and human activities were not considered.
Comparison and integration of multi-SDMs can help reduce uncertainty in the prediction of species distribution and can optimize prediction results compared with a single SDM. Guo et al. (2018) used 9 SDMs (MaxEnt, SRE, FDA, MARS, GBM, CTA, GLM, ANN and RF) to simulate suitable habitats for Populus euphratica in the Heihe River Basin of China. They found that the spatial distribution of suitable habitats for P. euphratica output from 9 SDMs varied widely, while only the simulation of MaxEnt model accurately depicted the distribution characteristics along the river for P. euphratica (verified by remote sensing images). The simulated results in our study also achieved a high accuracy, 83% observed sample points lay within the moderately and highly suitable areas under the current condition. Moreover, the probability response range (presence probability range [ 0.50) for annual mean temperature was 8.2-12.6°C, similar to the result of Su et al. (2005) (7.8-11.2°C).

Changes in the horizontal and vertical potential distribution of T. wallichiana
Compared with the current condition, there were significant and rapid decreases in the highly and moderately suitable habitats for T. wallichiana in the future climate scenarios. Poudel et al. (2014) also reported that the potential suitable areas of three yew species in the Central Himalaya would decrease substantially in the future (year 2080) under the A2A and B2A scenarios. Liu et al. (2013) predicted the potential suitable areas of two T. wallichiana lineages in the HHM region with the MaxEnt model and showed that the potential suitable areas in the Last Glacial Maximum (LGM, ca. 21 ka) were significantly greater than those in the Last Interglacial Period (LIG, ca.120 ka) and in the current climate conditions. Moreover, the potential suitable areas in the Last Interglacial were significantly less than the current climate condition because the temperature in the Last Interglacial was at least 5°C higher than the current temperature, which was consistent with the trend of reducing suitable area of T. wallichiana due to climate warming (Table 2). There was, however, a significant increase in suitable habitat area for T. wallichiana under the RCP8.5 scenario during the 2070s, and this can be related to the change of structure and composition of terrestrial vegetation under a high GHG (global greenhouse gas) emission scenario. IPCC (2013) reported that the rate of warming under the RCP8.5 scenario will be 65 times as high as the average warming during the last deglaciation, and the probability of large compositional changes and large structural changes in terrestrial vegetation are both greater than 60% under the RCP8.5 scenario (Nolan et al. 2018).
The study also showed that the mean elevations of suitable habitats for T. wallichiana would shift upslope in the future scenarios based on the outputs from multiple climate models. During the Quaternary climatic oscillations, the unique topography in the HHM region provided a climatic gathering place for plant migration (Liu et al. 2013). Global climate warming has promoted plants to migrate to higher elevations and latitudes and this has been confirmed by worldwide observations (Parmesan and Yohe 2003;Lenoir et al. 2008Lenoir et al. , 2009. Some studies have also predicted that global climate warming will continue to change the distribution of plants by forcing them to move to higher elevations (Huang et al. 2018;Du et al. 2018;Broadhurst et al. 2018). These shifts would lag behind those of climatic changes (Bertrand et al. 2011;Alexander et al. 2018) and significant shifts in elevation would occur only over the long-term (Schickhoff et al. 2015). With global climate warming, Chen et al. (2011) reported that the average upslope shift rate was ca. 11 m/10a based on a meta-analysis of 1367 species. Lenoir et al. (2008) observed that the average upslope shift rate of tree species was ca. 29 m/10a in Western Europe. Huang et al. (2018) predicted that future warming could cause Chinese sea buckthorn to shift upslope by 43-128 m/10a. Our results suggested that the potential distribution of T. wallichiana could shift upslope by 42-100 m/10a. The larger upslope shifts of suitable habitats in high elevation regions such as the HHM region were related to the strong sensitivity of habitat factors to global warming (Kohler et al. 2010) as well as monsoon interactions (Schickhoff et al. 2015;Huang et al. 2018).

Main climatic factors affecting the potential distribution of T. wallichiana
In general, terrestrial ecosystems are highly sensitive to temperature change (Nolan et al. 2018). As mean annual temperature increases, other ecologically important variables, such as seasonal temperatures, seasonal precipitation, climatic extremes, and degree of variability will change, often in multivariate, complex (Jackson and Overpeck 2000;Jackson et al. 2009), episodic, and nonlinear ways (Millar and Stephenson 2015;Johnstone et al. 2016).
In our study, annual mean temperature consistently contributed 61% towards the formation of potential suitable habitats for T. wallichiana. Moreover, variation in isothermality and temperature seasonality contributed 20%. This means that with warming climate, the shift in potential suitable habitats for T. wallichiana in the HHM region was mainly controlled by the annual average temperature and by temperature variation. On the other hand, the BIO1 (annual mean temperature) was strongly correlated with these variables including BIO5 (maximum temperature of warmest month), BIO6 (minimum temperature of coldest month), BIO8 (mean temperature of wettest quarter), BIO9 (mean temperature of driest quarter), BIO10 (mean temperature of warmest quarter), and BIO11 (mean temperature of coldest quarter), whose Pearson's correlation coefficients were 0.96, 0.99, 0.98, 0.98, 0.99, and 0.99, respectively (Table S4). This indirectly showed that extreme temperatures also played an important role in the formation of potential suitable habitats for T. wallichiana. In addition, this study showed that annual precipitation contributed 17% to the formation of suitable habitats for T. wallichiana. Su et al. (2005) also concluded that the three main climatic factors affecting the geographical distribution of T. wallichiana were low temperatures, humidity, and high temperature, in that order.

Management priorities of T. wallichiana under climate change
Areas of highly and moderately suitable habitats for T. wallichiana would sharply decline and would shift upslope in future climate scenarios. Therefore, the effective ex situ measures for the scattered mature individuals and seedlings in marginally suitable areas may need to plan based on field investigation and planting experiments.
Macroscopically, a map of potentially suitable areas for T. wallichiana was created by integrating multi-climate model results (Fig. S3-S6). The intersection area of the suitable area at different levels, for example, under the RCP4.5 scenario in 2050s, was more credible than the differential area (union-intersection) and should have priority for protection to ensure the survival of habitats and trees occupying areas that are currently over-exploited due to unprotected. The suitable areas at higher elevations should be prioritized for management. The differential area (union-intersection) has uncertainty, but it retains potential for secondary protection and development in terms of climate suitability. These areas take the intersection as the core area, expand to surrounding areas of mountains and foothills, extend to upper reaches of rivers, or continue to stretch over the high elevation mountainous areas. This accommodates the ecological requirements of T. wallichiana under climate warming. However, planting tests are needed to verify this.
The potential suitable habitats for T. wallichiana were mainly located in China, Nepal, India, and Bhutan in the HHM region, so the statutory protection and planting of T. wallichiana requires coordinated action between these countries. It is reported that the unlawful logging of yew in China has been severe for nearly a decade, and the number of criminal cases exceeded 1000 (Wang 2018). The exploitation of wild resources must be strictly prohibited through the implementation of laws and policies. Meanwhile, although the cultivation of yew and the industrial extraction of taxol was launched in the late 1990 s, the management of conservation, cultivation, and commercialization still need multi-sectoral efforts. At present, cultivation of yew can be one of the most effective means of utilizing sustainably and protecting the yews what remain in the wild (Thomas and Farjon 2011;TRAFFIC East Asia Report 2007). Among these, the delineation of planting areas will be an important link based on survival rate of yew with respect to climate change. For example, Duque-Lazo et al. (2018) found there was a significant linear positive correlation (r 2 = 0.414, p \ 0.05) between the survival rate of cork oak afforestation and the potential probability of presence by predicting the suitable area distribution of cork oak based on SDMs.

Conclusions
In the future (2050s, 2070s) scenarios (RCP2.6, RCP4.5, RCP8.5), the potentially highly and moderately suitable habitats for T. wallichiana would decline in area, and would shift upslope. Annual mean temperature and seasonal variability in temperature, combined with annual precipitation were the main climate variables affecting the habitat of T. wallichiana. Prior protected areas and suitable planting areas must be delimited from the future potential distributions, especially the intersection areas at different suitability levels. It is helpful to promote the sustainable utilization of this precious resource by prohibiting exploitation and ex situ restoring wild resources, as well as artificially planting considering climate suitability.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.