The possible impact of human activity and climate change on the potential suitable habitats of Taxus wallichianavar.mairei (Taxaceae) evaluated by ensemble modeling techniques

In the current work, we examined the possible changes in potential suitable habitats of Taxus wallichiana var. mairei (Taxaceae) under the influence of human activity and climate change using ensemble modeling techniques. T. wallichiana var. mairei is an endangered subtropical coniferous evergreen tree and is often used in medical applications. It is mainly distributed in southern China and has experienced several degenerations due to human activity in recent decades. Its growth is also very sensitive to climate change. The area of its highly suitable habitat under current climate conditions was approximately 2.31 × 106 km2, as estimated by our ensemble models, accounting for 81.85% of its total suitable habitat. The ensemble model results showed that the mean diurnal range (bio2) was the most critical environment variable affecting the performance of the ensemble models. Human activity had an overall negative influence on the suitable habitats of T. wallichiana var. mairei under current climate conditions. Under the influence of future climate change, the ensemble model predicted that the highly suitable habitat would convert to low- to medium-suitability or be lost, with the predicted loss of highly suitable habitat at the end of this century ranging from 30.40% to 96.55%. The predicted changes were more intense with increases in the severity of global warming. In addition, the ensemble model also predicted a shift in the distribution of total suitable habitats toward higher elevations. The results of this study should provide information for conservation purposes.


Introduction
Taxus wallichiana var.mairei (synonym: Taxus chinensis var.mairei; Taxaceae) is a subtropical coniferous evergreen tree mainly endemic to southern China (Zhang et al 2018, Rathore et al 2019, Wang et al 2019).T. wallichiana var.mairei is of great value in China because it can be used in many ways.Since it has a highly compact and solid wood texture, it is well known as an excellent timber tree (Cao et al 2016, Zhang et al 2018).T. wallichiana var.mairei is also a popular ornamental tree (also called 'beauteous yew' in China) because it can provide beautiful red seeds (arils) and evergreen leaves in autumn (Zhang et al 1994, Liao et al 2002, Zhang et al 2018).Moreover, the bark of T. wallichiana var.mairei contains Paclitaxel, which has been known to be an anticancer compound since the 1970s (Wang et al 2004, Gao et al 2007, Sun et al 2015, Shao et al 2021).Paclitaxel wallichiana var.mairei.Some previous studies have noted that due to the lack of reliable data, it is difficult to evaluate the impact of human activity (Ahmad et al 2020).
In this study, we first applied ensemble modeling approaches to SDMs to analyze the possible impact of human activity on the potential suitable habitats of T. wallichiana var.mairei.Then we used ensemble modeling techniques proposed in our previous work (Wang et al 2022) to both SDMs and climate models to investigate the possible influence of climate change on the potential suitable habitats of T. wallichiana var.mairei using four shared socioeconomic pathway (SSP) scenarios in four future periods.To our knowledge, this study is the first systematic analysis of the influence of both human activity and climate change on potential suitable habitats of T. wallichiana var.mairei.The results of this study may provide information for suitable habitat assessments and species conservation.

Occurrence records
The occurrence records of T. wallichiana var.mairei with location information in latitude and longitude were collected from the Global Biodiversity Information Facility (GBIF, www.gbif.org;Lane and Edwards 2007) and the Chinese Virtual Herbarium databases (CVH, www.cvh.ac.cn;Zhang et al 2004).There were 254 records obtained from GBIF and 76 records obtained from CVH.All duplicated records were removed and then spatially rarefied using SDM toolbox version 2.5 (Brown 2014) to ensure that the distance between any two records is larger than 10 km following Penado et al (2016).The rarefication process was conducted because SDM usually requires that input data to be spatially independent (Guisan et al 2017, Ahmed et al 2021).During this process, the records from CVH had a higher priority of remaining than those from GBIF.Finally, 177 records remained and were used to construct the SDMs.A site map of all the records used in the current work was produced using the National Center for Atmospheric Research (NCAR) Command Language (NCL, www.ncl.ucar.edu)(figure 1).

Human activity data
In this study, the 1 km spatial resolution Human Footprint (HF) data 2018 Release (Venter et al 2016), which provides the degree of human activity in 1993 (https://sedac.ciesin.columbia.edu/data/set/wildareas-v3-1993-human-footprint) and 2009 (https://sedac.ciesin.columbia.edu/data/set/wildareas-v3-2009-humanfootprint)(hf1993 and hf2009 hereafter), are used as human activity data in this study.This dataset integrates eight human-induced ecological pressures (i.e., built environments, population density, night-time lights, croplands, pasture, roads, railways and navigable waterways) into a single map.Although climate change is also affected by human activity (e.g., burning of fossil fuels; IPCC 2013, IPCC 2021), in this study, we do not include the climate change induced by humans as human activity.A large (small) value in this dataset represents high (low) ecological pressure from human activity.

Current environment variables
The 19 bioclimatic variables and three topographic variables (STable. 1) with a 2.5-minute spatial resolution used in this study are obtained from the World Climate Database version 2.1 (www.worldclim.org)(Hijmans et al 2005).These variables are regarded to be more biological relevant than the others and are widely used to build SDMs (Guisan et al 2017).Following Wang et al (2022), we used a variance inflation factor (VIF) analysis to select environment variables.The VIF can identify the hidden correlation structures between variables (e.g., Graham 2003, Dormann et al 2013, Guisan et al 2017, Qin et al 2017, Yan et al 2020, Zhao et al 2020, Qin et al 2020a, 2020b) and help to avoid the multicollinearity problem in SDM (e.g., Harrell 2001, Guisan et al 2002, Guisan et al 2006).The variables with a VIF value greater than 10 were removed (Guisan et al 2017).After VIF analysis, seven environment variables were selected and used in this study to construct SDM: bio 2 (mean diurnal range), bio 3 (isothermality), bio 8 (mean temperature of wettest quarter), bio 9 (mean temperature of driest quarter), bio 14 (precipitation of driest month), bio 15 (precipitation seasonality) and bio 18 (precipitation of warmest quarter), which is the same as Wang et al (2022).
To ensure that using human activity data and environment variable data together will not introduce hidden correlation structures, the hf1993 and hf2009 data were interpolated onto environment variable grids by bilinear interpolation.VIF was then applied to these environment variables.It can be noticed that using the human activity data and the environment variable data together did not introduce hidden correlation structures (SFig. 2 in the supplement material).
In this study, the multivariate environmental similarity surface (MESS) was adopted to analyze climate anomalies (e.g., Elith et al 2010, Feng et al 2019).A grid cell with climate anomalies suggests that the future climate in this grid cell is different from that at present, and climate change is thus significant in this grid cell.The formulae and calculation process of the MESS can be found in Elith et al (2010).In brief, when a grid cell has a MESS value of 100, it means that the future climate is identical to the current climate at this grid cell.When a grid cell has a positive MESS value, the larger the value is, the smaller the climate differences are; when a grid cell has a negative MESS value (climate anomaly), it indicates that there is one or more environment variable that are out of the current climate range.

Constructing ensemble SDMs
Ten techniques are used to construct SDMs (STable.3 in the supplement material).These techniques are provided by the R package 'Ensemble Platform for Species Distribution Modeling (biomod2)' (Thuiller et al 2016).The SDMs were trained and calibrated using 70% of the occurrence records with a 5-fold cross-validation and evaluated using the remaining 30% occurrence records (test dataset).In addition, 5 pseudo-absence datasets with sizes of 500, 1000, 2000, 5000 and 10000 were randomly generated and used to train the SDM, and 250 models (5 pseudo-absence datasets × 5 folds × 10 techniques) were produced.The performances of SDMs are evaluated with the receiver operator characteristic (ROC) and the true skill statistic (TSS).Typically, ROC > 0.9 or TSS > 0.8 are criteria used to select SDMs with excellent performance (e.g., Ben Rais Lasram et al 2010, Lin and Chiu (2019)).In this study, only those SDMs with both ROC > 0.9 and TSS > 0.8 were selected to create ensemble SDMs.
Three ensemble SDMs are produced in this study.Two ensemble SDMs are constructed with the seven environment variables and with hf1993 and hf2009, respectively (denoted as HF1993 and HF2009 hereafter).
The third ensemble SDM is constructed with the seven environment variables only, which is denoted as BIOCLIM in the following text.All the training and validation settings of the three ensemble SDMs were the same, and the only difference among them was whether hf1993 or hf2009 was used.There were 184 SDMs selected to form BIOCLIM and 166 SDMs selected to form HF1993 and HF2009 with the criteria ROC > 0.9 and TSS > 0.8.Comparing the predicted potential suitable habitats under current climate conditions between HF1993, HF2009 and BIOCLIM, we can quantify the possible impact of human activity on the potential suitable habitats of T. wallichiana var.mairei.Such a comparison can reveal how suitable habitat changes are associated with ecological pressure from human activity.The differences between hf1993 and hf2009 can also provide information on the extent of human activity's impact on potentially suitable habitats.
Usually, the SDM output must be converted into presence-absence maps when SDMs are used in conservation planning (Allouche et al 2006).Thus, habitat suitability was divided into four categories according to SDM prediction results: unsuitable (<0.2), low suitability (0.2 to 0.4), medium suitability (0.4 to 0.6), and high suitability (0.6 to 1.0) (e.g., Convertino et al 2014, Ansari and Ghoddousi 2018, Zhang et al 2019, Qin et al 2020b).In this study, regions that had high suitability (0.6 to 1.0) are referred to as highly suitable habitats, and regions that had both low suitability and higher suitability (>0.2) are referred to as total suitable habitats.
To evaluate the possible impact of future climate change on the potential suitable habitats of T. wallichiana var.mairei, BIOCLIM was projected onto 128 future climate model datasets (4 time periods × 4 future scenarios × 8 climate models).By analyzing the ensemble projection of BIOCLIM at different time periods under different future scenarios, we are able to quantify the possible impact of climate change on the potential suitable habitats of T. wallichiana var.mairei.The centroids of total suitable habitats were calculated based on the BIOCLIM predictions using SDM toolbox version 2.5 (Brown 2014).The changes in the centroids revealed how changes in climate conditions affected the total suitable habitat of T. wallichiana var.mairei over time.c)), it can be noticed that the ecological pressures from human activity were not uniformly changed spatially.For example, the human induced ecological pressures in Guangdong Province (the location of various provinces in China can be found in SFig.1), one of the most developed regions in China, decreased in 2009 compared with 1993.However, the ecological pressures from human activity increased from 1993 to 2009 in most regions of China.

Ensemble model performances
The ROC and TSS scores of the ten techniques (STable.3)used in BIOCLIM, as well as the scores of BIOCLIM, HF1993 and HF2009 are presented in SFig. 3 in the supplement material.Among the ten modeling techniques, Random Forest (RF) exhibited the best performance, while Surface Range Envelope (SRE) had the worst performance (SFig.3).The ensemble model BIOLCIM had a higher performance (represented by a red star in SFig.3) than any other individual technique, which had an ROC value of 0.97 and a TSS value of 0.85.The ensemble model HF1993 (represented by a yellow triangle in SFig.3) had an ROC score of 0.99 and TSS score of 0.87, which suggested that HF1993 performed slightly better than BIOCLIM.Similarly, the ensemble model HF2009 (represented by a solid blue circle in SFig.3) had an ROC value of 0.99 and a TSS value of 0.89, which outperformed BIOCLIM and HF1993.The better performances of HF1993 and HF2009 than BIOCLIM may be due to the additional variable of human activity, which was independent of the seven environment variables (SFig.2).Overall, the predictive skills of HF1993, HF2009 and BIOCLIM were high, and therefore were used in the following analysis.

Variable importance
To compare the relative importance of human activity and environment variables under the current climate conditions in the ensemble model performance, we calculated the mean contribution of each variable in BIOCLIM, HF1993 and HF2009.We then converted their contributions to percentages for better comparison (figure 3(a)).In BIOCLIM, the mean diurnal range (bio 2), the mean temperature of the driest quarter (bio 9), the precipitation seasonality (bio 15) and the precipitation of the warmest quarter (bio 18) made the most significant contributions, suggesting that they were the most important environment variables that significantly affected the BIOCLIM performance.In HF1993 and HF2009, the relative contributions to the performance from the environment variables are generally consistent with those of BIOCLIM.However, obvious changes in the contributions can also be noticed when human activity is considered.By adding the variable hf1993, the contributions of bio 3, bio 8, bio 15 and bio 18 were apparently increased, bio 14 was slightly increased, and the contributions of bio 2 and bio 9 were apparently decreased when compared with BIOCLIM (figure 3(b)).Compared to HF1993, when the variable hf2009 was added, the contributions of bio 3, bio 8, bio 15 and bio 18 were further increased.In contrast, bio 2, bio 9 and bio 14 were apparently decreased (figure 3(b)).As a result, for HF2009, the mean diurnal range (bio 2) and the precipitation of the warmest quarter (bio 18) made the largest contributions.
In sum, the relative contributions of the variables used in the three ensemble models are generally consistent with each other.The mean diurnal range (bio 2) made the largest contribution to all three ensemble models.Adding the human activity variable leads to apparent changes in the contributions of the environment variables, especially for the mean temperature of the wettest quarter (bio8) and the precipitation of the warmest quarter (bio18).

The impact of human activity on the potential suitable habitats
In this section, the possible influence of human activity on the potential suitable habitats of Taxus wallichiana var.mairei (Taxaceae) under current climate conditions are examined.BIOCLIM predicted the potential suitable habitat of T. wallichiana var.mairei under current climate conditions is shown in figure 4(a).The highly suitable habitat regions (0.6 to 1.0) were mainly located in southern China (south of the Yangtze River), covering Zhejiang, Anhui, Fujian, Jiangxi, Hubei, Hunan, Henan, Guangxi, Guangdong, Shanxi, Chongqing, Sichuan, Guizhou, Yunnan, Gansu and Taiwan.The highly suitable habitats of T. wallichiana var.mairei encompassed an area of approximately 2.31 × 10 6 km 2 , accounting for 81.85% of the total suitable habitats (>0.2).
The habitats of species can contract, expand or shift under the influence of human activity.The HF1993 (figure 4(b)) and HF2009 (figure 4(c)) predicted potential suitable habitats of T. wallichiana var.mairei generally matched that given by BIOCLIM (figure 4(a)).To see the possible influence of human activity, figure 5 depicts the differences in total and highly suitable habitats between HF1993 and BIOCLIM, as well as between HF2009 and BIOCLIM.Under the pressure of human activity in hf1993, the total suitable habitats mainly contracted north and south of those predicted by BIOCLIM (figure 5(a)).The contracted highly suitable habitats due to human activity in hf1993 were mainly located to the north of those given by BIOCLIM (figure 5(b)), and the expanded highly suitable habitats were mainly located to the south of those given by BIOCLIM (figure 5(b)).Some highly suitable habitats in the middle of the highly suitable habitats of BIOCLIM were contracted under the influence of human activity in hf1993 (figure 5(b)).The locations of the contracted and expanded total and highly suitable habitats under the impact of human activity in hf2009 were similar to those in hf1993, but the area of contracted total (highly) suitable habitats given by HF2009 was slightly larger (smaller) than that given by HF1993.
To obtain a better view of the changes in the total and highly suitable habitats associated with human activity in hf1993 and hf2009, the expanded and contracted areas were converted to percentages (figure 6).In general, under pressure from human activity, the contracted areas of total and highly suitable habitats were much larger than the expanded areas for both hf1993 and hf2009.For the total suitable habitat area, the contracted areas using hf2009 were more prominent than those using hf1993.However, for the highly suitable habitats, the contracted areas using hf1993 were larger than those using hf2009.In sum, human activity obviously has a negative influence on the suitable habitats of T. wallichiana var.mairei under current climate conditions.

The impact of future climate change on potential suitable habitats
To better understand the evolution of the total suitable habitats of T. wallichiana var.mairei associated with changes in the climate in the future, the MESS analysis results were first presented in figure 7.They showed that under the SSP126 and SSP245 scenarios, climate anomaly areas (MESS < 0) were quite small and limited in Xinjiang and Hainan Provinces.In the SSP370 and SSP585 scenarios, climate anomalies expanded from Hainan to its north in southern China with time.Under the most extreme warming scenario (SSP585), climate anomalies were distributed in Hainan, Guangxi, Guangdong, Hunan, Hubei, Jiangxi and Xinjiang, which indicated that climate change imposed large pressures on the suitable habitats of T. wallichiana var.mairei (figure 7(p)).
The possible impact of future climate change on potential suitable habitats of Taxus wallichiana var.mairei was examined (figure 8).A comparison to figure 4(a) shows that the suitable habitats under all four future scenarios tended to contract with time.In addition, comparing the suitable habitats at the same time periods but   under different scenarios, it is also clear that the more significant the increase in global temperature was, the more severely the suitable habitats were contracted.For example, in 2081-2100, there were much fewer highly suitable habitats under the SSP585 scenario than under the SSP126 scenario (figures 8(m)-(p)).
Further analysis is performed to investigate the differences in highly suitable habitats between the four future climate scenarios and the current climate (figure 9).Under the mildest warming scenario (SSP126), the highly suitable habitats mainly contracted in the south and west and slightly expanded to the north of the current highly suitable habitats (figures 9(a), (e), (i), (m)).The area of the expansion regions increased from 2021 to 2060, but it almost stopped increasing after 2061.In contrast, the area of the contracted regions continued to increase with time since 2021.In other warming scenarios, the more severely the temperature increased, the greater the area of highly suitable habitats contracted, and the contracted highly suitable habitats tended to increase with time.Meanwhile, the area of expanded regions in SSP245, SSP370 and SSP585 increased during the period from 2021 to 2060 and then decreased during the period from 2061 to 2100.However, the areas of expanded regions in SSP245, SSP370 and SSP585 were limited and almost the same.After 2061, the more severely the global temperature increased, the more highly suitable habitat regions decreased.In the most extreme warming scenario (SSP585), T. wallichiana var.mairei was predicted to lose most of its highly suitable habitats, with limited highly suitable habitats constrained in western Hubei and central Taiwan Province (figure 9

(p)).
To more clearly see the areas of the contracted and expanded highly suitable habitats, we converted the prediction results to percentages (figure 10).In general, with the increase in the severity of climate change, the areas of the contracted highly suitable habitats increase.Under the most extreme warming scenario (SSP585), T. wallichiana var.mairei was predicted to lose up to 96.55% of its current highly suitable habitats during the period 2081 to 2100.The expanded highly suitable habitats were limited and less than 4% for all four climate scenarios.These results suggested that T. wallichiana var.mairei would probably face severe pressure from future climate change due to the loss of highly suitable habitats.
The above analysis revealed that T. wallichiana var.mairei may lose most of its highly suitable habitats under the future severe climate scenario; however, it was not clear how the total suitable habitats (>0.2) would change.
To clarify this point, the difference in total suitable habitats between the future climate scenarios and the current climate is presented in figure 11.In the mildest scenario (SSP126), the area of expanded and contracted total suitable habitats tended to increase during the period from 2021 to 2060, but they were almost stable after 2061.In the worst scenario (SSP585), by the end of the 21st century, most contracted total suitable habitats were located in Guangdong, Guangxi, Fujian, Yunnan, Sichuan and Henan Provinces (e.g., figure 11(p)).The total suitable habitats tended to expand in the northwest (e.g., figure 11(p)).In other climate scenarios, the more severe the temperature increases were, the more total suitable habitats expanded and contracted.Additionally, under all scenarios except SSP126, the area of expanded and contracted total suitable habitats increased with time.In the most warming scenario (SSP585), Yunnan, Henan, Guangdong, Guangxi and Fujian were predicted to lose most of their total suitable habitats, while most expanded total suitable habitats were found in western Sichuan and Xizang Provinces (figure 11(p)).Comparing figure 11(p) and figure 1, the total suitable habitats expanded toward higher elevations.Comparing figures 9 and 11, vast regions with highly suitable habitats (e.g., Zhejiang, Jiangsu, Anhui, Jiangxi, Hubei, Hunan, Guizhou and Chongqing) will convert to medium-and lowsuitability habitats at the end of this century under the influence of climate change under the SSP585 scenario.
As shown in figure 10, we counted the expanded and contracted total suitable habitat areas and converted them to percentages (figure 12).Under SSP126, the area of contracted and expanded total suitable habitats was predicted to increase until 2080 and then begin to decrease.In other scenarios, the area of both contracted and expanded total suitable habitats was predicted to increase with time.Again, the more severe the temperature increases were, the greater the predicted expansion and contraction.In the most extreme warming scenario (SSP585), T. wallichiana var.mairei was predicted to lose 38.47% of its current total suitable habitats and expand 15.94% of its current total suitable habitats.
The centroids of the total suitable habitats under current climate conditions and future scenarios at different time periods are presented in figure 13.Under the most extreme warming scenario (SSP585), T. wallichiana var.mairei was predicted to lose most of its current highly suitable habitats; therefore, we only displayed the centroid shift of total suitable habitats.Figure 13 shows that the predicted centroids in the period from 2021 to 2040 under all climate scenarios were located to the north of the current centroid, suggesting a tendency for shifting to higher latitudes.In addition, with the increase in the severity of global warming, the centroids of total suitable habitats shifted toward the northwest of the current centroid, and the more intense the global temperature increases were, the greater the centroid moved toward northwestern China (higher elevation, figure 1).Overall, future climate changes were predicted to shift the total suitable habitats toward a higher latitude and higher elevation.

Discussion
Our ensemble model predicted that under the current climate conditions, the highly suitable habitats (0.6 to 1.0) of T. wallichiana var.mairei.were mainly aggregated in southern China (south of the Yangtze River), consistent with the documentation in the literature (e.g., Li et al 2007).In addition, our model predictions showed that the low-and medium-suitability habitats (0.2 to 0.6) were mainly distributed in Jiangsu, western Yunnan and southeastern Xizang.There were some plant bases in Jiangsu and western Yunnan, as documented in the literature (e.g., He et al 2017).Previous studies reported that T. wallichiana var.mairei was discontinuously distributed over the Hindu Kush-Himalaya and adjacent regions (e.g., Poudel et al 2012), and some previous studies also predicted suitable habitats in similar regions in Xizang (e.g., Li et al 2021).These results verified the reliability of our model predictions.In addition, compared to previous SDMs constructed using a single modeling technique, our ensemble SDMs showed an obviously better performance in terms of ROC and TSS skills, indicating that the ensemble SDMs had a higher prediction reliability than individual SDM, consistent with previous research (e.g., Convertino et al 2014, Ansari and Ghoddousi 2018, Zhang et al 2019, Qin et al 2020b).
By introducing human activity as an additional variable to BIOCLIM, our ensemble HF1993 and HF2009 had higher prediction skill than BIOCLIM.This was because hf1993 and hf2009 were independent to the other variables, which would provide additional information for the ensemble SDMs.Results showed that 12.45% (7.88%) and 7.49% (9.82%) highly (total) suitable habitats contracted, while 1.94% (1.55%) and 4.06% (2.24%) highly (total) suitable habitats expanded due to human activity in 1993 and 2009, respectively.Overall, human activity imposes a negative effect on the potential suitable habitats of T. wallichiana var.mairei.Figure 9.The predicted changes in the current highly suitable habitats in the future by BIOCLIM.Blue denotes that the highly suitable habitats can remain at that time period, pink denotes that the current highly suitable habitats are lost at that time period, and green denotes that the regions become highly suitable at that time period.2) in the future by BIOCLIM.Blue denotes that the highly suitable habitats can remain to that time period, pink denotes that the current highly suitable habitats are lost at that time period, and green denotes that the regions become highly suitable at that time period.Wang et al (2019) predicted that the suitable habitats of T. wallichiana var.mairei would substantially shrink in the 2070 s, while Li et al (2021) predicted that the current suitable habitats of T. wallichiana var.mairei would expand in the 2050 s under the mildest warming scenario but would contract in the same period under other severe warming scenarios.Guo et al (2022) predicted that T. wallichiana var.mairei would lose 28.9% of its current suitable habitats at the end of this century under the mildest warming scenario using a climate suitability model and a distribution limitation model.In this study, we applied ensemble techniques to ten modeling techniques and eight climate models to systematically predict the potential suitable habitats of T. wallichiana var.mairei in four time periods under four future climate scenarios.Our ensemble model results showed that even under the mildest warming scenario (SSP126), T. wallichiana var.mairei was predicted to lose up to 26.87% of its current highly suitable habitats at the end of this century.The greater the severity of global warming increases was, the more current highly suitable habitats were predicted to be lost.In the most extreme warming scenario (SSP585), T. wallichiana var.mairei was predicted to lose 95.75% of its current suitable habitats at the end of this century.In addition, our ensemble model predictions suggested that vast highly suitable habitats in China (e.g., Zhejiang, Jiangsu, Anhui, Jiangxi, Hubei, Hunan, Guizhou and Chongqing) would convert to low-and medium-suitability habitats at the end of this century.Under SSP585, T. wallichiana var.mairei was predicted to lose 22.53% of its current total suitable habitats in the 2100 s.The low-and medium-suitability habitats were predicted to expand toward higher elevation regions with time.Overall, with the increasing temperature in the future, the highly suitable habitats were predicted to contract or convert to low-and medium-suitability habitats and the low-and medium-suitability habitats were predicted to expand toward a high elevation.The total suitable habitats were predicted to shift toward higher latitudes and higher elevations.
Some previous studies predicted that the tropics might expand under the influence of global warming (Quan et al 2014, Tselioudis et al 2016).As a result, the subtropical evergreen broad-leave forests (as mentioned in the introduction, T. wallichiana var.mairei mainly grows in this kind of forest) in southern China may move northward toward central China (Weng and Zhou 2006).However, T. wallichiana var.mairei mainly grows on mountains and requires taller trees to serve as nurse plants, and regions north of the Yangtze River are mainly plains and lack mountains, T. wallichiana var.mairei is unable to move northward with broad-leave evergreen forests.Therefore, south of the current highly suitable habitats of T. wallichiana var.mairei will be lost under future climate scenarios.Some studies reported that broadly distributed species usually have higher ecological adaptability (the capability of species to address environment variations) than narrowly distributed species (e.g., Zhang et al 2015, Ma andSun 2018).The conversion of highly suitable habitats to low-and medium-suitability habitats (rather than being lost) and the expansion of low-and medium-suitability habitats toward higher elevation regions under the severe global warming scenario were probably due to the high ecological adaptability of T. wallichiana var.mairei.Overall, in the future, continuously increasing temperature was predicted to have a negative effect on T. wallichiana var.mairei.
Our ensemble models revealed that, among all the environment variables used in this work, the mean diurnal range (bio2) has the largest contribution to SDMs.This result suggested that temperature was a critical factor that affected the distribution of T. wallichiana var.mairei.Xie (2014) identified that the annual average temperature was the most important factor affecting the distribution of T. wallichiana var.mairei using principal component analysis.Both high and low temperatures have a negative effect on the distribution of T. wallichiana var.mairei, and the best temperature for its growth is approximately 16 °C.Cao et al (2016) also noted that temperature, heat, humidity, and sunshine were dominant factors affecting the distribution of T. wallichiana var.mairei.Their experimental study revealed that temperature had a larger effect on the accumulation of primary metabolites than other factors, such as radiation Li (2017).These works support the influence of temperature variation on T. wallichiana var.mairei.Our results further revealed that the mean diurnal range was important for T. wallichiana var.mairei distribution, indicating that T. wallichiana var.mairei needed a certain range of temperatures for growth, which agreed with the findings of previous studies.In addition, the predicted contracted suitable habitats under the pressure of human activity were smaller than those under the pressure of future climate change.These results suggest that the negative effect imposed by climate change on suitable habitats may be larger than that imposed by human activity.

Conservation recommendations
Based on our analysis, we have the following conservation recommendations: (1) Human activity generally imposes a negative influence on the suitable habitats of T. wallichiana var.mairei.
Therefore, human activity should be regulated and controlled, which may help to conserve suitable habitats of T. wallichiana var.mairei; (2) In the most severe warming scenario (SSP585), T. wallichiana var.mairei may lose most of its highly suitable habitats.Although there are numerous nature preserves located in the current highly suitable habitats, these nature preserves may convert to low-or medium-suitability habitats or be unsuitable for T. wallichiana var.mairei growth.Model prediction showed that some highly suitable habitats remained near the boundary of Hubei and Chongqing (red ellipse in SFig.4) and Taiwan (yellow ellipse in SFig.4).Therefore, we suggest that nature reserves should be established in the above two regions.

Limitations
This study has several modeling limitations: (1) Our ensemble SDMs do not consider the biological characteristics of T. wallichiana var.mairei.As mentioned in the introduction, T. wallichiana var.mairei usually grows under taller trees, which serve as its nurse plants (Zhang et al 2018).The seed germination and seedling survival rates of T. wallichiana var.mairei is low (Li and Fu 1997, Fu et al 1999, Zhang et al 2018).Lei et al (2011) revealed that the larvae of Scarabaeoidea can damage the roots of T. wallichiana var.mairei.Fusarium oxysporum can cause root rot of T. wallichiana var.mairei (Qiao et al 2015).These biological factors can also affect the distribution of suitable habitats; however, how these factors affect the suitable habitats under current and future climate conditions remains unclear.
(2) There is a lack of human activity data under future scenarios.Therefore, it is unable to compare the influence of human activity and climate change on the suitable habitats of T. wallichiana var.mairei in future climate.The discussion of how suitable habitats are affected by human activity is limited under current climate conditions.

Conclusion
In the current work, ensemble modeling techniques were applied to SDMs and climate models to investigate the influence of human activity and climate change on the potential suitable habitats of T. wallichiana var.mairei in China.The mean diurnal range (bio2) was identified as the most important factor affecting the distribution of T. wallichiana var.mairei.Overall, human activity imposed a negative effect on the total and highly suitable habitats, but such a negative effect may be relatively smaller than that of climate change.The highly suitable habitats were predicted to contract in the future under the pressure of climate change, regardless of which climate scenario was applied in the ensemble model.Meanwhile, the more severe the global warming was, the more highly suitable habitats were predicted to be lost or converted to low-and medium-suitability habitats.
Under the most extreme warming scenario, T. wallichiana var.mairei was predicted to lose most of its highly suitable habitats.The present work also provided some suggestions for conservation strategy makers to alleviate the negative influence of human activity and climate change on the distribution of T. wallichiana var.mairei.

Figure 1 .
Figure 1.Topography (elevation above sea level in meters, color shading) and occurrence records of T. wallichiana var.mairei used in this study (red dots).

Figure 2
shows the map of hf1993 (figure 2(a)) and hf2009 (figure 2(b)) and their differences (figure 2(c)), which were created by NCL.Generally, central China has a greater ecological pressure from human activity than other regions (figures 2(a) and (b)).Comparing the human activity in 1993 and 2009 (figure 2(

Figure 2 .
Figure 2. Map of human activity in (a) 1993 and (b) 2009.(c) Changes in human activity between 2009 and 1993.

Figure 3 .
Figure 3. (a) The contribution of each variable in BIOCLIM, HF1993 and HF2009.(b) The variable contribution differences between HF1993 and BIOCLIM, and between HF2009 and BIOCLIM.The y-axis represents the (a) variable contribution and the (b) variable contribution differences in percentage.

Figure 5 .
Figure 5.The predicted changes in the (a), (c) total suitable habitats (suitability > 0.2) and (b), (d) highly suitable habitats (suitability > 0.6) between (a), (b) HF1993 and BIOCLOM and (c), (d) HF2009 and BIOCLIM.Blue denotes that suitable habitats can remain, pink denotes that suitable habitats are lost, and green denotes that the regions become suitable under the influence of human activity.

Figure 6 .
Figure6.Percentage changes in the total suitable habitats (suitability > 0.2) area and highly suitable habitats (suitability > 0.6) area between HF1993, HF2009 and BIOCLIM.The y-axis represents the changes of corresponding area in percentage.

Figure 7 .
Figure 7. MESS value maps of different future scenarios at different time periods.

Figure 10 .
Figure10.Percentage changes in the current highly suitable habitats area in the future.'-' represents the habitat contracted case, and '+' represents the habitat expanded case.The y-axis denotes the changes of corresponding habitats area in percentage.

Figure 11 .
Figure11.The predicted changes in the current total suitable habitats (suitability > 0.2) in the future by BIOCLIM.Blue denotes that the highly suitable habitats can remain to that time period, pink denotes that the current highly suitable habitats are lost at that time period, and green denotes that the regions become highly suitable at that time period.

Figure 12 .
Figure12.Percentage changes in the current total suitable habitats (suitability > 0.2) in the future.'-' represents the habitat contracted case, and '+' represents the habitat expanded case.The y-axis denotes the changes of corresponding habitats area in percentage.

Figure 13 .
Figure 13.Changes in total suitable habitat (suitability > 0.2) centroids in the future.