Incorporating Local Adaptation Into the Forecasts of Paeonia Veitchii Distribution

Local is and the new approach of separate SDMs for separate populations is suggested. The implementation of a conservation strategy for PV_Daodi based on a population model could include in-situ conservation and assistance in northward migration. The predicted spatial and temporal pattern of range shifts will be a useful reference for the cultivation and conservation of high-quality P. veitchii. Our study provided a cornerstone for the distribution prediction of high-quality medicinal plants.


Background
Species distribution models (SDMs) typically link the presence (or absence) of species in multiple locations to associated environmental covariables to estimate habitat preferences or predict their distribution [1,2]. SDMs, as a useful tool for mapping the potentially suitable habitat of species, have become a crucial tool in conservation, ecology, biogeography, and wildlife management research.
However, the effect of local adaptation is largely ignored in most prediction studies, in which all populations of a species are assumed to respond consistently to the environmental conditions experienced by the entire species. Local adaptation is widely observed [3,4] , and several attempts have recently been made to incorporate local adaptation into ecological analyses under climate change [5][6][7].
These studies indicated that populations can be discriminated based on their response to climate change. However, there is little knowledge about local adaptation in SDM concerning traditional Chinese medicine.
Daodi herb is a concept that has been widely recognized in Chinese medicinal history for centuries. The quality and clinical effects of Daodi herbs surpass those of the same botanical origin produced from other regions (non-Daodi). Thus, this herb is widely recognized and has long enjoyed a good reputation [8,9]. The ecological quality of traditional Chinese medicine is a new frontier in the eld of ecological research. This eld is a comprehensive and cross-applied discipline composed of the quality and ecology of traditional Chinese medicine. This emerging discipline also explores the interaction and mechanism between the quality of Chinese medicine and the ecosystem [10][11][12][13]. Paeonia veitchii is one of the major and essential Daodi herbs in Sichuan Province, China. The main chemical compounds found in P. veitchii are glycoside compounds, such as paeoni orin, hydroxy-paeoni orin, paeonin, benzoylpaeoni orin, oxypaeoni orin, and benzoyloxypaeoni orin [14], which have many pharmacological effects, including antioxidative, anti-thrombotic, anti-atherosclerotic, immunomodulatory, and anti-in ammatory effects [15][16][17]. Paeoni orin is the main component of P. veitchii, and the paeoni orin content of P. veitchii from Daodi (PV_Daodi) is relatively higher than that of P. veitchii from non-Daodi (PV_non-Daodi) [18]. Previous studies primarily paid attention to chemical components, pharmacological effects, pollen fertility, and photosynthetic physiology [19,20]. Little is known about its habitat distribution and the driven eco-factors that form the suitable habitat of PV_Daodi and PV_non-Daodi at the population level.
In this study, principal component analysis (PCA) and hierarchical cluster analysis (HCA) were employed to delineate the PV_Daodi population with presence sites in Sichuan Province and the PV_non-Daodi population with presence sites in non-Sichuan places. Then, we used maximum entropy algorithm (MaxEnt) to simulate the migration trend of the potential distribution of each population under several climate change scenarios. The aims of the present study were to 1) identify probable speci c populations with different local adaptations, 2) identify the essential climate variables that shape the distribution range of each population, and 3) predict the change in the potential habitat distribution of each population under global climate change. This study is the rst to combine the prediction of the geographic distribution of a species with the Daodi herb concept. These results would provide a reasonable basis for the prediction of the habitats of medicinal plants with high quality and clinical effects.

Target species and occurrence data.
We obtained the event records of P. veitchii through an extensive literature search in three databases, namely, the Global Biodiversity Information Facility (https://www.gbif.org/), the National Specimen Information Infrastructure of China (http://www.nsii.org.cn), and the Chinese Virtual Herbarium (http://v5.cvh.org.cn/). We used Google Earth (http://ditu.google.cn/) to determine the latitude and longitude of an area when the exact geographic coordinates are missing. Reasonable sampling site selection was performed based on the following principles: only one record was reserved when the site data have replicates, very close records were deleted according to the resolution of environment variables, and space was ltered so that only one point occurred in each grid cell (10 km×10 km). These operations minimized the spatial autocorrelation of the sample points and reduced their impact on the model's results. A total of 226 speciation records (Table s1) were obtained for analysis using the given process ( Fig. 1). The distribution records of the target species were sorted as ".CSV" format les according to the requirements of MaxEnt software.

Environmental parameters
The relevance and completeness of predictors are the key to building SDMs [21]. In this study, we preliminarily selected 31 environmental factors that may have an impact on the current species distribution patterns. These factors included 19 bioclimate variables downloaded from the World Climate database (www.worldclim) [2], and 12 global UV-B radiation variables from the gIUV database (http://www.ufz.de/gluv/) [22]. Highly correlated variables would mislead the importance of variables and response curves; therefore, Spearman's rank correlation was used to test the cross-correlation (ρ<0.80), and highly correlated environmental factors (Supplementary Fig. 1 and Fig. 2) were removed. Eight of the 31 variables were selected as evaluator variables (Table 1). All environment variables were preprocessed in American Standard Code for Information Interchange format with the aid of ArcGIS.

Population Grouping
All occurrence sites were grouped into PV_Daodi and PV_non-Daodi populations according to the theory of Daodi herbs in traditional Chinese medicine. PCA [23] and HCA [24] were employed to check whether different climate conditions exist between the two populations, and PCA was conducted in all eight variables using Simca-P. The rst two principal components were employed to de ne the populations. The SDM of both populations (PV_whole) was also constructed using all presence sites together in addition to the separate SDMs of the two populations. This practice is common when conducting SDMbased research. In other words, we constructed SDMs for two approaches: one that produced separate SDMs for separate populations (PV_Daodi and PV_nonDaodi) and another that produced a single SDM for the entire species that encompassed both populations (PV_whole).

Past and future climate scenarios
The Last Glacial Maximum community climate (LGM, ca. 22 ka) and the Mid-Holocene (MH, ca. 6 ka) data sets in the Community Climate System Model version 4.0 (CCSM4; http://www.worldclim.org) were employed to forecast the possible distribution areas in the past. CCSM4 is one of the most effective climate models that simulate the effects of past and future climate changes on plant distribution [25]. This model has been successfully applied to similar studies [26].
Future bioclimatic variable data in CCSM4 under three representative concentration pathways (RCPs), namely, RCPs 2.6, 6.0, and 8.5, were used to project the future climate situations according to human economic activity, land use patterns, climate policy, and other factors. These RCPs include a stringent mitigation scenario (RCP 2.6), one intermediate scenario (RCP 6.0), and one scenario with extremely high greenhouse gas emissions (RCP 8.5). We used the climate data of two periods in the 21st century (the 2050s and 2070s, which are the averaged data for 2041-2060 and 2061-2080, respectively) to project habitat distribution changes under each RCP.

Habitat suitability modeling.
The MaxEnt model used was based on the free software version 3.3.3k (http://www.cs.princeton.edu/_schapire/maxent/). The algorithm calculates the most likely potential geographical distribution of a species according to the relationship between geographical data and the known distribution of the target species [27].
We utilized 75% of the occurrence data to calibrate the model and the remaining 25% of the data to test the predictive ability of the model during the modeling procedure. Threshold independent receiver operating characteristic (ROC) analyses were carried out within MaxEnt software, and the areas under the ROC curve (AUC) were calculated to check the performance of the model [28]. The AUC ranges from 0.5 (random) to 1.0 (perfect discrimination), and values > 0.8 indicate satisfactory performance [29,30]. We used the permutation importance and jackknife measures [31] in MaxEnt to understand which climate variables may be important and differ among populations, evaluate the relative contribution of each environmental variable to identify the most important variables, and determine the predicted distribution of the modeled entity. The response curves of all variables in the models were identi ed and examined to check the response of each population to different variable values.

Ecological factors of different populations
The key thresholds of the ecological factors of species were extracted from ArcGIS following the sampling sites of P. veitchii. The violin plots of the ranges of the ecological factor values of PV_Daodi and PV_non-Daodi are shown in Fig. 2, which shows that the ecological factor values varied considerably. The isothermality (Bio3), annual temperature range (Bio7), precipitation of the wettest month (Bio13), precipitation of the driest month (Bio14), mean global UV radiation in January (UV_JAN), and mean global UV radiation in July (UV_JUL) of PV_Daodi and PV_non-Daodi showed signi cantly different values (p<0.05). This result indicated that the populations of Daodi-and non-Daodi-producing areas require different environmental and climatic conditions.

Population Grouping
The distribution of the 226 presence sites was drawn based on the rst two principal components, which explained 78.6% of all variance (Fig. 3). Despite some overlap between the two populations, PV_Daodi and PV_non-Daodi were located separately in PCA maps: PV_non-Daodi occupied a niche featured by annual mean air temperature (Bio1), Bio3, high Bio13, Bio14, UV_JAN, and UV_JUL, whereas PV_non-Daodi was featured by mean diurnal range (Bio2) and Bio7. The HCA result showed that all sites of occurrence can be grouped into three groups: two groups for PV_non-Daodi and one group for PV_Daodi.
The combined results from HCA and PCA indicated that divergence in climate adaptation exists between the two populations despite the little overlap between these populations. In other words, the population from Daodi-processing area experienced a different climate compared with those in other areas.   according to the MaxEnt model under current climatic conditions. The logical output generated by the MaxEnt software is expressed in probability and ranges from 0 to 1. The model results were divided into four levels using the reclassi cation tool of ArcMap 10.2: 0-0.2 is not suitable, 0.2-0.4 is generally suitable, 0.4-0.6 is moderately suitable, and 0.6-1 is highly suitable [32].

3.4Predicting distribution under current conditions
The suitable areas for PV_Daodi were China, India, Bhutan, the United States, Afghanistan, and Syria with smaller distributions in Lebanon, Israel, Nepal, and Burma. The highly suitable areas were concentrated in Eastern China (Fig. 5a). The current areas of the generally, moderately, and highly suitable habitats for PV_Daodi were 3.8×10 5 , 2.6×10 5 , and 1.1×10 5 km 2 , respectively. The suitable areas for PV_non-Daodi were distributed in China with sporadic distributions in Afghanistan, Bhutan, India, Iran, Iraq, Lebanon, Morocco, Spain, Syria, Tajikistan, and the United States. The highly suitable areas were greatly distributed in China (Fig. 5b). The current areas of the generally, moderately, and highly suitable habitats for PV_non-Daodi were 6.4×10 5 , 4.0×10 5 , and 2.0×10 5 km 2 , respectcively. The suitable areas for PV_whole were distributed in China, Afghanistan, India, Bhutan, United States, Lebanon, Syria, Morocco, Burma, and Vietnam, and the highly suitable area is in China (Fig. 5c). The current areas of the generally, moderately, and highly suitable habitats for PV_whole were 8.6×10 5 , 5.4×10 5 , and 2.7×10 5 km 2 , respectively.
3.6 suitable habitat centroid shift from past to future climate.
The centroid was calculated to fully understand the suitable distribution shifts, and the centroid migration path was drawn to illustrate the direction and distance shift of the population under different climate scenarios (Fig. 7). The centers of the suitable distribution area are from Libya and Gaza in LGM to Iran in MH and to China and India at present. The centers will move northward toward Central Asia under the three future climates. The magnitude of migration varies slightly: the future centers of the PV_Daodi population are Iran and Afghanistan, the future centers of the PV_non-Daodi population are Turkmenistan and Afghanistan, and the center of the PV_whole population is Afghanistan. Notably, PV_Daodi will move a smaller distance to the north than PV_non-Daodi under the future climate environments. This difference may be caused by the high temperature of the carbon dioxide emission model, which PV_non-Daodi will not be able to adapt.

Differences Between Populations
P. veitchii is a Daodi herb primary distributed in Sichuan Province, China for thousands of years. Therefore, we hypothesized that the population growth of P. veitchii in Sichuan and non-Sichuan provinces would experience different degrees of adaptation to their local environments. The PCA and HCA results partly supported our hypothesis. PCA results suggested that P. veitchii could be diverging because of different climate variables, such as UV, rainfall, and temperature. A former work evaluated the impact of climate change at the species level. Here, we individually evaluated the impact represented by the theory of traditional Chinese medicine, the Daodi herb, at the population level. This impact is directly related to the quality of the herb. We found that P. veitchii were climatically diverged and responded differently to past and future climates. The pattern and magnitude of the impact of climate change on suitable habitat at the population level uncovered in this study would not be detectable at the species level. We believe that this study will provide a scienti c basis for developing management strategies for the adaptation of Chinese medicinal plants to climate change.

Difference Between three SDMs
We constructed three SDMs for P. veitchii. The rst SDM followed common practice and treated all presence sites as a whole (PV_whole). The other two SDMs assumed that Daodi herb would adapt locally to their local environments. The PCA and HCA results veri ed this hypothesis; thus, separate SDMs were built for separate populations, namely, Daodi-(PV_Daodi) and non-Daodi-processing area populations (PV_non-Daodi). The past, current, and future distributions of P. veitchii Lynch forecasted by these three SDMs were consistent with each other with slight difference. Figure S7 shows that PV_whole failed to predict two countries that both independent SDMs predicted in LGM under RCPs 2.6 and 8.5 in 2070. PV_whole also forecasted some countries that the two separate SDMs did not detect. Therefore, we suggest that the separate SDM approach can correctly evaluate the suitability of a species in regions where local adaptation exist and can provide more details about species distribution than using a single SDM for all populations. The results indicated that climate change led to the differences among the SDMs and special consideration should be given to the selection of suitable methods for predicting the suitable distribution areas of Daodi herbs.

Conservation Strategy.
Different changes occurred to the suitable distribution areas of the two populations under the in uence of global warming. Most of the suitable habitats of PV_Daodi will decrease in 2070 compared with 2050 under RCPs 2.6 and 6.0 but will increase in 2070 compared with 2050 under RCP 8.5. The suitable area of PV_non-Daodi will increase in 2070 compared with 2050 under RCPs 6.0 and 8.5 and decrease in 2070 compared with 2050 under RCP 2.6. However, the suitable area of PV_whole will continue to increase regardless of the scenario; thus, this SDM indicates that global warming might bene t P. veitchii. This result is consistent with other studied peony species, namely, P. delavayi, P. rockii, and P. veitchii, which will have expanded highly suitable habitats in the 2070s under RCP 8.5 according to their models [33,34].
The ranges of the centroid shifts of PV_Daodi, PV_non-Daodi, and PV_whole are different, and the distance of the shifts between PV_non-Daodi and PV_whole is larger than that between PV_Daodi and PV_whole.
The implementation of a conservation strategy for PV_Daodi according to a population model could include in-situ conservation. Most of the areas currently suitable for planting will remain suitable for growing PV_Daodi in the future. Therefore, efforts should be made to conserve these areas from threats posed by human activities. Northward migration should be facilitated by 2070 based on the predicted future climate distributions.

Conclusion
Estimating how climate change will affect the distribution of Daodi herbs is crucial for their speci c cultivation. Our results indicated that local climate adaption exists within the distribution range of P. veitchii. Bio13, UV_JAN, and UV_JUL were de ned as the key elements shaping the distribution of PV_Daodi, and UV_JAN, Bio13, and Bio3 are important when shaping the distribution of PV_non-Daodi.
Local adaptation is worth considering in SDM research, and the new approach of constructing separate SDMs for separate populations is suggested. The implementation of a conservation strategy for PV_Daodi based on a population model could include in-situ conservation and assistance in northward migration. The predicted spatial and temporal pattern of range shifts will be a useful reference for the cultivation and conservation of high-quality P. veitchii. Our study provided a cornerstone for the distribution prediction of high-quality medicinal plants.