Potential distribution of Notopterygium incisum Ting ex H. T. Chang and its predicted responses to climate change based on a comprehensive habitat suitability model

Abstract Notopterygium incisum Ting ex H. T. Chang is a rare and endangered traditional Chinese medicinal plant. In this research, we built a comprehensive habitat suitability (CHS) model to analyze the potential suitable habitat distribution of this species in the present and future in China. First, using nine different algorithms, we built an ensemble model to explore the possible impacts of climate change on the habitat distribution of this species. Then, based on this model, we built a CHS model to further identify the distribution characteristics of N. incisum‐suitable habitats in three time periods (current, 2050s, and 2070s) while considering the effects of soil and vegetation conditions. The results indicated that the current suitable habitat for N. incisum covers approximately 83.76 × 103 km2, and these locations were concentrated in the Tibet Autonomous Region, Gansu Province, Qinghai Province, and Sichuan Province. In the future, the areas of suitable habitat for N. incisum would significantly decrease and would be 69.53 × 103 km2 and 60.21 × 103 km2 in the 2050s and 2070s, respectively. However, the area of marginally suitable habitat would remain relatively stable. This study provides a more reliable and comprehensive method for modelling the current and future distributions of N. incisum, and it provides valuable insights for highlighting priority areas for medicinal plant conservation and resource utilization.


| INTRODUC TI ON
The growth, development, and reproduction of plants are limited by climatic and other environmental factors. Global warming is likely to result in changes in biological habitats, losses of regional species diversity, and increases in the risk of species extinction (Anderson, 2013;Marzloff et al., 2018;Pacifici et al., 2015;Urban, 2015).
Topography affects the redistribution of moisture and heat in the natural environment; therefore, topography is an important driver of plant species distributions, especially for alpine plants, and in mountainous regions, macrotopographies are usually large enough to provide refuge for plant species under changing climates (Myan, Walker, & Paramor, 2013). The Qinghai-Tibet Plateau uplifted more than 3,000 m in the Quaternary period, which dramatically changed the topography and climate in this region and formed various climate types (Zhang, Fengquan, & Jianmin, 2000). Meanwhile, microtopography will hinder the ability of a species to shift poleward and upslope because it can cause a relatively stable and closed climate environment within a short distance (Patsiou, Conti, Zimmermann, Theodoridis, & Randin, 2015). Hence, microtopography created a biological refuge for many rare, relict, and endemic plants during the last glacial maximum (Elsen & Tingley, 2015;Yang, Zhou, Li, Song, & Chen, 2011). Global warming will inevitably affect the living environment and hydrothermal conditions of plants in the Qinghai-Tibet Plateau region, and the complex topography will exacerbate this impact .
Models are useful tools for simulating the impact of future climate change on plant species distribution, especially in this study, which considered a large spatiotemporal scale. In recent years, species distribution models (SDMs) have been popular tools for assessing the spatial-temporal variations in species distributions under different climate scenarios (Anderson, 2013;Mcquillan & Rice, 2015;Zhang et al., 2014). In the past two decades, advancements in statistical methods have promoted the development of SDMs, and numerous statistical methods and software programs have been developed to describe the niche characteristics of species and predict species distribution patterns. The popular algorithms are as follows: surface range envelope (SRE, i.e., BIOCLIM) (Booth, Nix, Busby, & Hutchinson, 2014), flexible discriminant analysis (FDA) (Basile et al., 2016), generalized linear model (GLM) (Lopatin, Dolos, Hernández, Galleguillos, & Fassnacht, 2016), generalized additive model (GAM) (Muñoz-Mas, Papadaki, Martinez-Capel, Zogaris, & Ntoanidis, 2016), multiple adaptive regression splines (MARS) (Friedman, 1991), generalized boosting model (GBM) (Moisen et al., 2006), classification tree analysis (CTA) (Thuiller & Lavorel, 2010), artificial neural network (ANN) (Segurado & Araujo, 2004), random forest (RF) (Mi, Huettmann, Guo, Han, & Wen, 2017), and maximum entropy (MaxEnt) (Phillips, Anderson, & Schapire, 2006). However, differential niche requirements of species shape the geographic distribution of species within an environment. Hence, the following factors should be considered when selecting an appropriate SDM algorithm for species distribution research on specific spatiotemporal scales, that is, species niche characteristics, especially those related to specific traits of species; environmental characteristics, especially the limiting factors in species' habitats; and data quality, including data availability and the spatial and temporal resolution of data (Bell & Schlaepfer, 2016;Guisan, Thuiller, & Zimmermann, 2017;. Simultaneously, for the same species with the same input data, different response curves and variable weights can be adopted according to the corresponding statistical properties of different algorithms, which will lead to different simulation results (Guo, Li, Zhao, & Wei, 2018) and will notably increase the uncertainties in the predicted species distribution. Therefore, the ensemble model (EM) strategy is proposed to solve this problem.
This strategy applies several SDM methods within a consensus modelling framework and reduces the predictive uncertainty of individual models by combining the results from multiple models (Araújo & New, 2007;Guisan et al., 2017). Previous studies have indicated that EMs can substantially improve model accuracy and applicability (Grenouillet, Buisson, Casajus, & Lek, 2011).

Notopterygium incisum
Ting ex H. T. Chang, also known as the "Emissary of Qiang Nationality" and "messenger of King Hu," is a Chinese herbal medicine associated with many beautiful and moving folk stories, and it is also an endangered traditional Tibetan medicinal plant (Committee of Flora of China, 1986). According to field investigations and the literature, N. incisum is mainly distributed in Shaanxi Province, Sichuan Province, Qinghai Province, and the Tibet Autonomous Region (Committee of Flora of China, 1986). Modern science has confirmed that N. incisum can be used to relieve inflammation, cure arrhythmia and myocardial ischemia, promote cerebral circulation, remove thrombosis and bacteria, and carry out other pharmacological actions (Committee of National Pharmacopoeia, 2015;Li, Zhang, Wang, & Lei, 2003).
Due to the increasing demands in the domestic and international medical markets and increasingly intense human activities, N. incisum is becoming increasingly threatened, primarily due to habitat disturbance and destruction (Liu, 2006;Sun et al., 2015;Zhou et al., 2003). In addition, this plant has been included on the Red List of endangered species in China since 2005. Thus, information on the potential geographic distribution and its response to climate change is vital to the protection and resource utilization of N. incisum.
Existing studies on N. incisum mainly focus on the identification of medicinal components, pharmacological analyses, cultivation, industrialization, and resource investigations (Gao & Fang, 2007;Li et al., 2003;Peng, Dong, Zhu, & Yan, 2006), while the emphasis of our research is on the dis-

| Occurrence data
In this research, the occurrence data for this species were from a variety of sources, including published scientific literature, field survey reports (Jiang et al., 2006;Liu, 2006;Sun et al., 2015), and online plant databases, such as the Global Biodiversity Information Facility (GBIF, 2019) and Chinese Virtual Herbarium (CVH, 2019). To improve the geographic accuracy and perform scientific modelling, we selected only data with precise longitude and latitude information and removed duplicate coordinates and incomplete information.
Finally, we obtained 99 occurrence data points to build the model ( Figure 1, Table S1).

| Pseudo-absence data
In practice, most SDM models require species absence data; however, true absence data are usually not available, so we use pseudoabsence data as a substitute (Zhang et al., 2016). The numbers and qualities of pseudo-absence data affect the accuracy of the SDMs (Zhang et al., 2016). In this research, using the "Sample by Buffered Local Adaptive Convex-Hull" tool in the SDMtoolbox 2.0 (Brown, Bennett, & French, 2017), we generated pseudo-absence data (Barbet-Massin, Jiguet, Albert, & Thuiller, 2012) within 200 km buffers around the occurrence points. We selected three random subsets of the background to generate three groups of pseudo-absence data, and each group included 500 points ( Figure 2).

| Environmental variables
In this study, we chose four categories of environmental datasets with a total of 24 environmental variables to characterize the environmental demands of N. incisum (Table 1) The topographic variables, including elevation, slope, and aspect, have the same resolution as the bioclimatic variables. The elevation variable was also acquired from the WorldClim dataset, and the slope and aspect variables were generated by the ArcGIS spatial analysis function based on the elevation variable.
In this research, the 1:1 million soil database of China was used for the soil type variables, and the 1:1 million China vegetation data were used for the vegetation type variables. Both of these data types were acquired from the National Tibetan Plateau Data Centre TA B L E 1 Environmental element index used for predicting the potential geographic distribution of Notopterygium incisum the same resolution as the other environmental variables. Because changes in vegetation and soil types lag behind climate change (Wu et al., 2015), we used the same soil type and vegetation type data in the future climate change scenarios.

| Methods
The "biomod2" package is an object-oriented, expandable, and repro- The accuracy of each model algorithm was evaluated by the area under the receiver operating characteristics (ROC) curve (AUC) and true skill statistic (TSS) (Thuiller et al., 2009). Then, we selected the algorithm that met the precision requirements to construct an EM to simulate the migration trend of N. incisum. Finally, we used biotic variables (vegetation and soil) to further refine the EM results.

| Input data processing
In this research, the processing of input data included the removal of collinearity between environmental variables and the generation of pseudosampling points. First, principal component analysis (PCA) (Supplementary materials 1) and Pearson correlations were used to select a subset of environmental variables (Guisan et al., 2017;Guo, Li, Zhao, & Nawaz, 2019). Finally, we selected seven bioclimatic variables and two topographic variables (i.e., annual mean air temperature, bio1; temperature seasonality, bio4; mean temperature of the warmest quarter, bio10; mean temperature of the coldest quarter, bio11; annual precipitation, bio12; precipitation seasonality, bio15; precipitation of the coldest quarter, bio19; slope; and aspect) (Table 1). Second, three sets of pseudosampling points were randomly generated, and there were 500 pseudosampling points in each group. We obtained three sets of model inputs using the three sets of pseudosampling points and the N. incisum occurrence data.

| Single modelling technique
We separately input the three sets of data into the model. For each modelling process, the sampling data of the N. incisum points involved in the modelling were divided into two parts. Seventy percent of the sampling data were used as a training set, and the remaining data were used as testing data. In addition, we set an equal weight for the total presence of sampling points and pseudosampling points.
We used TSS and AUC to evaluate the model, and Equation (1) for TSS is as follows: where a refers to the number of true positives, b refers to the number of false positives, c refers to the number of false negatives, and d refers to the number of true negatives. TSS values > 0.6 are considered good (Allouche, Tsoar, & Kadmon, 2006;Jones, Acker, & Halpern, 2010). The AUC indicates the area under the ROC curve, and the ROC plots display the relationship between sensitivity (Equation (2)) and 1-specificity (Equation (3)) over a range of threshold values (0-100). For AUC, which does not require selection of a habitat suitability threshold, values > 0.9 are considered good, and values from 0.7 to 0.9 are considered moderate (Allouche et al., 2006;Jones et al., 2010).
In addition, for each modelling algorithm, we repeated this modelling process 10 times with a bootstrapping sampling strategy; thus, in total, we built 300 single models (three sets of sampling data × 10 single modelling techniques × 10 repeats). Then, we constructed the EM using Equation (5):

| Ensemble model
where BT i refers to the potential habitat suitability index of the evaluation unit (grid) i, w j refers to the weight of the results of model j, and

| Comprehensive habitat suitability model
In its natural habitats, N. incisum is not the dominant species and is usually an associated species that grows in alpine shrublands, (1) TSS = Sensitivity + Specificity − 1 (2) Sensitivity = a∕(a + c) alpine meadows, and woodlands (Mcquillan & Rice, 2015). Hence, the vegetation types are the important limiting factors, reflecting the limitations of the species' migration ability, and here, according to the literature, we conducted binary conversion of the vegetation type data. A suitability value of 1 was defined as a suitable vegetation type for N. incisum, namely, alpine forest land, shrubland, and alpine meadows, while all other vegetation types received a value of 0. Because there are no clearly defined soil types that are suitable for the growth of N. incisum, we used the soil type variables and occurrence data of N. incisum to create a MaxEnt model to assess the soil suitability requirements. The model settings were as follows: 75% of the occurrence data were used as a training set, and the remaining data were used as testing data. Ten replications were performed, and the AUC was used to evaluate model performance.
Then, based on the EM and considering the effect of soil and vegetation conditions, we built a CHS model to further identify the distribution characteristics of N. incisum-suitable habitats, and the index of CHS for N. incisum was defined using the following equation (Equation (6)

| Model performance
The statistical accuracy results for the 10 models showed that RF was the best model, and the TSS and AUC values were 0.929 and 0.984, respectively. This model was followed by FDA and GBM, and the TSS values for these two models were above 0.90. The accuracy of the SRE was lowest, and the TSS and AUC values were 0.50 and 0.75, respectively (Figure 3). In addition, in the nine models involved in the EM, the average TSS value of each single model was >0.74, and the average AUC value of each single model was >0.91 ( Figure 3). Therefore, the EM provided satisfactory results, and the TSS and AUC values were 0.80 and 0.96, respectively. The AUC value of the MaxEnt model for the soil environment requirements was 0.87, which means that the model results were scientific and reasonable.

| Distribution of suitable habitats in the current climate environment
In this study, different statistical algorithms led to different simulation results, but the spatial distribution patterns of N. incisum were consistent. In addition, all model results showed that the suitable distribution area for N. incisum was concentrated in the transitional zone from the first step to the second step in the "ladder topog- which had a sporadic distribution. In addition, the habitats of N. incisum were fragmented with large local patches of suitable habitat that were not spatially contiguous. We also calculated the areas of suitable habitats and marginally suitable habitats for N. incisum, and the results showed that the suitable habitats had an area of approximately 83.76 × 10 3 km 2 , and the marginally suitable habitats had an area of approximately 102.72 × 10 3 km 2 (Table 2). In summary, the proportions of N. incisum-suitable habitats were small and narrow, and the distribution region was mainly around the eastern Qinghai-Tibet Plateau.

| Model rationality
Currently, several studies have been devoted to identifying the current distribution range of N. incisum (Shang et al., 2015;Sun et al., 2015), but there are many shortcomings in these studies, and one of them is the use of a simple algorithm, which can delimit only the general distribution scope and cannot provide many landscape details (Sun et al., 2015). Other studies have focused on only the local distribution of this species (Shang et al., 2015). In contrast to previous studies, we used an EM based on nine different models with different mathematical algorithms to simulate the potential distribution of N. incisum throughout China. This model strategy avoids the selection of a single best model, thus eliminating (or at least limiting) model selection bias and, even more importantly, reducing the uncertainty of the modelling results caused by different modelling algorithms (Guisan et al., 2017;  thus, at the regional scale, the vegetation and soil conditions dictate the distribution of the species, even in a favorable climate.
Thus, during the modelling process, the soil and vegetation types were used as overlying variables to further test the prediction results, and including these variables made the prediction results more realistic.

| Habitat area estimates for Notopterygium incisum
According to previous studies, the optimal potential habitats for N. incisum in some parts of Sichuan Province, Tibet Autonomous Region, Qinghai Province, and Gansu Province have areas >1.42 × 10 5 km 2 (Sun et al., 2009). The most suitable habitat is in Sichuan Province, and more than 60% of the most suitable area for N. incisum is located in western Sichuan, such as in Aba and Ganzi prefectures (Sun et al., 2009(Sun et al., , 2015. The EM results showed that the suitable habitats for N. incisum were mainly distributed in Sichuan Province, Gansu Province, Qinghai Province, Tibet Autonomous Region, and Shaanxi Province (Figure 4). In addition, the current suitable habitat area for N. incisum in mainland China covers 83.76 × 10 3 km 2 , and the suitable habitats are mainly distributed in Sichuan Province, with an area of approximately 38.71 × 10 3 km 2 , which is approximately 8.02% of the total area of Sichuan Province (Table 2). Our model result is consistent with that of previous studies in terms of the spatial range (Sun et al., 2009(Sun et al., , 2015 but provides more details about the distribution of N. incisum habitats. The model results indicated that by the 2070s, the areas of suitable habitats will have persistently decreased (Figures 4-6). Moreover, the suitable habitats for N. incisum are narrowly distributed, and in the centralized distribution area, most suitable habitats will remain roughly unchanged in the future. This

| The dominant environmental index response to suitability
Previous studies have shown that cold tolerance, growth-season temperatures, and available water supply are the main environmental indexes that influence the distribution of alpine vegetation (Woodwand, 1987). In this study, according to the model results, the major variables that affected the distribution of N. incisum habitats were annual mean air temperature, temperature seasonality, mean temperature of warmest quarter, mean temperature of coldest quarter, annual precipitation, and precipitation seasonality. To further clarify the suitable ranges of environmental factors and the critical values, we drew the response curves for each variable mentioned above ( Figure 7) and defined the suitable ranges of the variables (logistic probability of presence > 0.3). The results showed that the suitable ranges were −2.5 to 16.9°C for annual mean air temperature, 6.0-8.6°C for temperature seasonality, 6.8 -16.3°C for mean temperature of warmest quarter, −12.8 to 7.7°C for mean temperature of coldest quarter, 419-1181 mm for annual precipitation, and 71.4-98.5 mm for precipitation seasonality. The habitat characteristics of N. incisum, such as cold humid climates, were summarized through these environmental indexes. The results were more detailed than those of previous studies (Shang et al., 2015;Sun et al., 2015).
Climate change will influence N. incisum physiology through several interrelated processes, which will lead to changes in N. incisum distributions. As a perennial herb, freezing winter temperatures (i.e., bio11) at middle to high latitudes will harm the underground rhizome of N. incisum, which will be a main limiting factor affecting the N. incisum distribution; however, this situation

F I G U R E 7
The dominant environmental index response curve. Annual mean temperature (bio1), temperature seasonality (bio4), mean temperature of warmest quarter (bio10), mean temperature of coldest quarter (bio11), annual precipitation (bio12), and precipitation seasonality (bio15) may change in some areas with a warming climate. In addition, climate change will contribute to the increase in accumulated temperatures and lengthen the duration of plant growth (Zhang, Li, Ding, & Zheng, 2017). As a medicinal plant, the root and rhizome have been used in traditional Chinese medicine; hence, with increasing accumulated temperatures, the production of medicinal material will increase. Furthermore, in high and cold regions, increasing temperatures will lengthen the duration of soil moisture necessary for plant growth and reduce frost damage. However, if temperatures exceed their physiological optimum range, the growth and reproduction of N. incisum will be affected, and suitable habitats will gradually disappear.
Similarly, changes in precipitation patterns can influence the N. incisum distribution. Global warming will lead to an increase in precipitation in western China, especially the Qinghai-Tibet Plateau; hence, the water supply available for plants will increase.
Nevertheless, prior work has suggested a relatively low impact of precipitation on the distribution of N. incisum; rather, the impact of temperature change will be the strongest (Shang et al., 2015;Sun et al., 2015). Indeed, in this study, nearly all 10 single models showed that among the nine environmental variables, all four temperature variables had key impacts on the distribution of N. incisum, and for the EM, the contribution rates of the temperature variables were higher than those of the precipitation variables, which indicated that the N. incisum distributions were less sensitive to precipitation variation than temperature variation.

| CON CLUS IONS
Numerous studies have indicated that species will shift their geographic distributions poleward and upslope due to the dramatic topographic changes in the Qinghai-Tibet Plateau; this shift will be complex and heterogeneous. Here, we conducted an analysis of the

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

AUTH O R CO NTR I B UTI O N S
ZFZ, YLG, HYW, and WG conceived and designed the study. ZFZ and YLG collected the samples and established the model. ZFZ, QR, JL, and QZZ analyzed the data. ZFZ, YLG, HYW, and WG wrote the paper. All authors read and approved the final product.

DATA AVA I L A B I L I T Y S TAT E M E N T
All environmental variables used in the manuscript are already publicly accessible, and we have provided the download address in the manuscript; relevant sampling site information can be found in Table   S1 in the online version.