Evaluation of Limiting Climatic Factors and Simulation of a Climatically Suitable Habitat for Chinese Sea Buckthorn

Chinese sea buckthorn (Hippophae rhamnoides subsp. sinensis) has considerable economic potential and plays an important role in reclamation and soil and water conservation. For scientific cultivation of this species across China, we identified the key climatic factors and explored climatically suitable habitat in order to maximize survival of Chinese sea buckthorn using MaxEnt and GIS tools, based on 98 occurrence records from herbarium and publications and 13 climatic factors from Bioclim, Holdridge life zone and Kria' index variables. Our simulation showed that the MaxEnt model performance was significantly better than random, with an average test AUC value of 0.93 with 10-fold cross validation. A jackknife test and the regularized gain change, which were applied to the training algorithm, showed that precipitation of the driest month (PDM), annual precipitation (AP), coldness index (CI) and annual range of temperature (ART) were the most influential climatic factors in limiting the distribution of Chinese sea buckthorn, which explained 70.1% of the variation. The predicted map showed that the core of climatically suitable habitat was distributed from the southwest to northwest of Gansu, Ningxia, Shaanxi and Shanxi provinces, where the most influential climate variables were PDM of 1.0–7.0 mm, AP of 344.0–1089.0 mm, CI of -47.7–0.0°C, and ART of 26.1–45.0°C. We conclude that the distribution patterns of Chinese sea buckthorn are related to the northwest winter monsoon, the southwest summer monsoon and the southeast summer monsoon systems in China.


Introduction
Chinese sea buckthorn (Hippophae rhamnoides subsp. Sinensis) is a subspecies that belongs to the family Elaeagnaceae [1,2], and is also commonly called acid vinegar salix or black thorn. It has the largest distribution and planting area of all species belonged to the same genus, and it is the most important economic and ecological shrub species. Thus far, the ecological effects of Chinese sea buckthorn have been well documented. It improves soil properties, reduces pollution and prevents soil erosion. Economically, the plant species has not only been considered as an ornamental shrub but also been used as a source of nutritious food, medicine and firewood [3][4][5]. Therefore, Chinese sea buckthorn has been arousing international attention as a promising new crop, and it is predicted by some as the next major health food fad [3]. Since 1985, Chinese sea buckthorn has been used to control accelerated desertification under an initiative of the government and the cultivation of this species has been widely promoted in northwest China. In the past 20 years, the planting area of Chinese sea buckthorn has reached approximately 14,667.4 km 2 (22 million acres) [6], and increased by more than 800 km 2 annually [7]. Detailed information of ecological characteristics, suitable habitat and potential distribution of Chinese sea buckthorn is needed for improvement of scientific research and general cultivation of this plant species.
Climate has been documented to play a critical role in determining large-scale species distributions [8][9][10][11]. Climatic factors are considered effective surrogates for resource needs or other limiting factors, which may have a stronger forecast ability than biophysical factors in predicting species distribution and growth [12,13]. Gu et al. [14] studied the factors that determined the distribution of Chinese sea buckthorn on the Tibetan Plateau and found that precipitation was the key factor affecting distribution of the plant species. Lian and Chen [15] believed that Chinese sea buckthorn has a stable distribution area, which could be used as an indicator for China's three major vegetation zone boundaries (forest, steppe and alpine zones). These studies greatly improve our understanding of the influence of climatic factors on geographic distribution of Chinese sea buckthorn. However, these studies are based on empirical data, and they failed to provide/quantify the relative importance of climatic factors, response curves and climatic thresholds of this species.
The current dominant method to explore climate-species distributions is using species distribution models (SDMs), which is based on species niche theory [16][17][18][19]. Several methods are available to simulate species distribution patterns [20,21], and these models can be classified into mechanistic modeling approaches and statistical modeling approaches. Mechanistic modeling approaches mainly simulate the fundamental niche of a species, while statistical modelling approaches investigate the realized niche of a species [22,23]. Mechanistic modeling approaches, such as the CLIMEX model [24][25][26][27][28][29], have been used extensively to simulate various species of agricultural and non-agricultural crops as well as pests. But Mechanistic modeling approaches are based on many eco-physiological parameters, which are rarely obtained for a majority of species, including Chinese sea buckthorn. Statistical modeling approaches can be classified into two groups: those requiring presence-absence data and those requiring presence-only data [22,30,31]. The latter group has been developed recently and is gaining popularity because of the availability of presence-only data. Financial constraints, difficulty in conducting field inventories and the lack of reliable or meaningful absence data has also made the models requiring presence-only data more appealing (e.g. genetic algorithms [32], maximum entropy approach [33,34], bioclimatic envelope [35] and distance algorithms [36]). Such modeling techniques are recommended when absence data of a species is not available or reliable [30,37]. Although there is high consistency between models, many studies have shown that maximum entropy modeling (MaxEnt) is widely used and usually produces accurate predictions of species distributions [20,21,38]. Furthermore, Elith et al. [39] expanded the ability of the MaxEnt model by increasing the limiting factor mapping and similar surface mapping for range-shifting species, which is especially suitable for predicting potential species distributions and interpretation of limiting factors.
The purpose of our study was to examine how climate factors limit the distribution of Chinese sea buckthorn and to simulate the climatically suitable habitat for Chinese sea buckthorn by utilizing the available data of species occurrence, features of the MaxEnt model and a geographic information system (GIS). Our study mainly focuses on the following objectives: (1) to identify the important climatic factors and where the limiting factors affect the distribution boundaries of this species; (2) to identify the climatic thresholds and find out the climatically suitable habitats for Chinese sea buckthorn; and (3) to explain the distribution patterns based on the perspective of East Asian monsoon climate. Our study could provide theoretical support for top-level design of the cultivation and planting of Chinese sea buckthorn for policy-makers and planners.

Study area and species records
The study was conducted in China, which is located in the eastern part of Asia, on the western coast of the Pacific Ocean (3°52 0 -53°33 0 N, 73°40 0 -135°2 0 E), and spans are over 5,000 km from east to west and over 5,500 km from north to south. The summer monsoons dominate the climate of the country, and the topography descends in a three-step staircase from west to east (Fig 1). China covers a land area of about 9.6 million km 2 , occupying one-fifteenth of the world's total land area. Because of its large area and extent, as well as its high population, China has a wide variety of climates, topography, soil and vegetation types [40].
Species records of Chinese sea buckthorn in continental China were mainly collected from published resources [14,15,41] and the Chinese Virtual Herbarium [42], which integrates herbarium data of national natural museums from 14 institutes in China. We found 368 specimens, and then removed duplicate specimens and specimens with no location/coordinates information. The longitude and latitude of the records were accurate at the county level. Consequently, 97 records were collected from all sources. Spatial distribution pattern of 97 records relating to topography of China can be seen in Fig 1 and latitude/longitude coordinates of each record was stored in an Excel database (S1 File) for MaxEnt model building.

Climatic variables
Many climatic variables have been used to characterize hydrological-thermal climatic niches in China, such as BIOCLIM variables [43,44], Holdridge life zone variables [8] and Kira's index variables [45,46]. These variables have been widely used in research on the relationship between species/vegetation and climate at a regional or global scale. The BIOCLIM variables have been widely used in SDMs studies, as the data can be easily downloaded from the World-Clim database with no further calculation requirement. Holdridge life zone and Kira's index variables are seldom used in SDMs studies, but they still provide considerable power for explaining species distributions. Here, we integrated these climatic variables and selected 13 climatic variables (Table 1), including 5 variables from Holdridge life zone and Kira's index variables and 8 variables from Bioclim variables, mainly by reducing redundant information in Bioclim variables and considering extreme climatic variables that usually limit species distribution [47]. Baseline climatic layers were downloaded from the WorldClim database at a spatial resolution of 10 arc-min, which were generated by using thin-plate smoothing splines with latitude, longitude, altitude, and monthly temperature and precipitation records from 1950-2000 from climate stations [43,44].

Model selection and evaluation
This study utilized MaxEnt software (version 3.3), a machine-learning algorithm designed by Phillips et al. [33,34,48]. MaxEnt expresses the suitability of a grid cell as a function of the environmental or climate layers at that grid cell in a landscape, together with a set of sample locations where the plant species has been observed. The distribution of suitability is proved to be the same as the Gibbs distribution that maximizes the product of the probabilities of the sample locations, where the Gibbs distribution takes the form (Eq 1): Here c 1 , c 2 , c 3 ,. . . are constants, f 1 , f 2 , f 3 ,. . . are the features, and Z is a scaling constant that ensures that P sums to 1 over all grid cells.
Elith et al. [39] enhanced the ability to explain the MaxEnt model outputs by integrating multivariate environmental similarity surface (MESS) analysis and limiting factors mapping in the new version of MaxEnt (3.3). MESS analysis assists in revealing whether there is possible model-predicted novel habitat (extrapolation), which informs the credibility of model output. If we let min i be the minimum value of variable V i over the reference point set and the similarly for max i , p i be the value of variable v i at point P, f i be the percent of reference points whose value of variable V i is smaller than p i , then the similarity of P with respect to variable V i is calculated by Eq 2: Finally, the multivariate similarity of P is the minimum of its similarity with respect to each variable. The MESS is similar to a BIOCLIM analysis [35] that makes use of percentiles, but it is extended to give negative values so as to differentiate levels of dissimilarity when outside the range of the reference points, where the reliability of model output depends on the accuracy of species response curves. MESS results allow the mapping of locations where limiting factors are important. The assumption is that information on which variable is driving the MESS value at any given point can be extracted and mapped by finding the most dissimilar variable. The most dissimilar variable for a point P is the variable with the smallest value of similarity to P [39]. According to map of limiting factors, we could easily deduce which factors are limiting physiological and ecological processes.
The MaxEnt model is usually evaluated by the threshold-independent receiver operating characteristic (ROC) approach [calculating the area under the ROC curve (AUC) as a measure of prediction success], which is a graphical method that represents the relationship between the false-positive fraction (1-specificity) and the sensitivity for a range of thresholds [49]. It has a range of 0-1, with a value greater than 0.5 indicating a better-than-random performance event. A rough classification guide is the traditional academic point system [50]: poor (0.5-0.6), fair (0.6-0.7), good (0.7-0.8), very good (0.8-0.9), and excellent (0.9-1.0). Simulating process and statistical analysis Coordinates of the Chinese sea buckthorn presence points were uploaded into the MaxEnt software, together with all of the 13 climatic variables. Linear, quadratic, product, threshold, and hinge features were used to generate feature types. A linear feature is simply one of the continuous environmental variables. A quadratic feature is the square of one of the continuous environmental variables. A product feature is the product of two continuous environmental variables; when used with linear and quadratic features, product features constrain the output distributions to have the same covariance for each pair of environmental variables as the samples. A threshold feature is derived from a continuous environmental variable. For a threshold value v, the threshold feature is binary (taking values 0 and 1) and is 1 when the variable has value greater than v. A hinge feature is also derived from a continuous environmental variable. It is like a linear feature, but it is constant below a threshold v. The convergence threshold (10 −5 ), maximum number of iterations (500), and 10,000 global background points were used to run the MaxEnt model. The logistic output (scaled up in a non-linear manner) was used to estimate the probability of presence (ranging from 0-1), which was easier to use and interpret than the other two outputs (raw and cumulative formats). A jackknife test (systematically leaving out each variable) and the regularized gain change [log of the number of grid cells minus the log loss (average of the negative log probabilities of the sample locations)] were then used to evaluate which climatic factors were the most important in determining the potential distribution of the species. Response curves were created by using the MaxEnt model only with the corresponding variable to avoid the potential problem of strongly correlated factors on the explanation of species response curves. MESS analysis and the limiting mapping technique were used to determine the location of novel habitat and where the limiting factors affected the distribution boundary of the species (each reference point indicates the occurrence records for the species).
Considering the large amount of occurrence data, we used a 10-fold cross-validation (sometimes called rotation estimation) method to evaluate the performance of the MaxEnt model. This method randomly partitioned the original sample into 10 subsamples of equal size. Of the 10 subsamples, 1 subsample was used as testing data, while the remaining 9 subsamples were used as training data. The cross-validation was then repeated 10-fold, and each observation was used for validation exactly once. The AUC values produced from the 10-fold cross-validation were then averaged to indicate the performance of the MaxEnt model. Cross-validation is important in testing model performance, especially where further sample collections are hazardous, costly, or impossible to collect.
A suitable habitat map for Chinese sea buckthorn was produced by utilizing the AUC weight averages of the 10 logistic output maps produced by 10-fold cross-validation, in which the relative suitability range was from 0 to 1. Normally, a threshold of 0.5 is set, which means that a grid cell with a value equal to or greater than 0.5 is recorded as presence, while a value smaller than 0.5 is recorded as absence. In order to retain the maximum predictive information, we divided the predicted suitability into 4 habitat class levels: core area (0.75-1.00), moderately suitable area (0.50-0.75), marginal area (0.25-0.50), and unsuitable area (0.00-0.25). The climatic thresholds (conditions) for each habitat class were analyzed using a GIS tool by ArcGIS 9.3 [51].

Current and potential distribution of Chinese sea buckthorn
According to the locations of Chinese sea buckthorn occurrence records, a map showing current distribution of this plant species was produced (Fig 1). Chinese sea buckthorn mainly occurs in seven provinces: Xizang, Sichuan, Qinghai, Gansu, Ningxia, Shaanxi and Shanxi. When the maps obtained from the 10-fold cross-validation were AUC weight averaged, a climatically suitable habitat map for Chinese sea buckthorn was created in Fig 2. The potential distribution of Chinese sea buckthorn is mainly located in semi-humid and semi-arid regions (Fig 2D), and begins from southwest of the Hengduanshan Mountains and extends to the southern boundary of the Daxinganling Mountains (extending southwest to northeast at 45°). The altitude decreases from the southwest (about 3500-4100 m) in the Hengduanshan Mountain region to northeast (1000-1500 m) in the Loess Plateau region. The width of the range is about 560 km with an area of approximately 1.12 million km 2 , occupying 11.7% of the land area of China, spanning 12 provinces including Xizang, Yunnan, Sichuan, Gansu, Qinghai, Ningxia, Shaanxi, Shanxi, Hebei, Inner Mongolia, Henan, and Liaoning with the core areas mainly located in Xizang, Sichuan, Gansu, Qinghai, Ningxia and Shaanxi provinces (Fig 2A). The distribution is located in the transition zone of China's three major vegetation zones (forest, steppe and alpine zones), primarily in the transitions between subtropical evergreen broadleaf forest region to alpine meadow and steppe region, and temperate deciduous broad-leaved forest region to the temperate steppe region (Fig 2C) [40,52].
Based on the suitability level map of Chinese sea buckthorn simulated by MaxEnt (Fig 2), Four habitat categories were defined: the core area (0.75 to 1.00), the moderately suitable area (0.50 to 0.75), the marginal area (0.25 to 0.50), and the unsuitable area (0.00 to 0.25). The climatic thresholds for the habitat categories are shown in Table

Model performance and importance of climatic factors
The 10-fold cross-validation method and areas under the receiver operating characteristic curve method were used to assess the accuracy of the MaxEnt model. The AUC values of test data and training data based on 10-fold cross validation in ascending order by training AUC value were shown in Fig 3, which indicated that the MaxEnt model performance was highly accurate (AUC value greater than 0.9), with an average test AUC value of 0.93 (0.9077 to 0.9745) and an average training AUC value of 0.95 (0.9471 to 0.9523). The coefficient of variation in test AUC values was only 2.2% among 10 model simulations, which showed that the 10-fold cross-validation method didn't influence the predicted performance of the MaxEnt model.   The relative contribution of climatic variables (in descending order) in determining the potential distribution of Chinese sea buckthorn was shown in Fig 4. The results indicated that PDM, AP, CI and ART were the most influential climatic factors in limiting the potential distribution of Chinese sea buckthorn and these 4 factors explained 70.1% of the variation (12.4-21.1% for each factor), followed by PWM, WI, PSD and MTCM, which explained another 24.9% of the variation (4.3-7.9% for each factor). The remaining 5 climatic factors were unimportant in limiting the distribution of the species, and they accounted for 5% of the variation (0.1-1.8% for each factor). These significant factors can be divided into hydrological-related factors (PDM, AP, PWM and PSD, accounting for 55.6% of the variance) and thermal-related factors (CI, ART, WI and MTCM, accounting for 39.4% of the variance). The results indicated   Table 1. The results clearly illustrated that PDM was best described by exponential decay and CI showed an exponential increase relationship, whereas AP and ART showed a bell curve relationship. Response peaks to habitat suitability for PDM, AP, CI and ART were 0 mm, 540 mm, 0°C and 38°C, respectively.

Discussion
SDMs have become the subject of an active field of research in large-scale ecology and biogeography, and they have been used to solve many ecological issues in recent decades [19,22,23]. These models have been used for biodiversity assessment, biological reserve design, habitat management and restoration, invasive species and pest threat management, species introduction and cultivation, and predicting the effects of climate change on species distributions [53][54][55][56][57][58][59]. The assumption underlying the SDMs is that species distribution is in equilibrium with climatic conditions and that the sample records used for the training model are sufficient and they represent the niche of the species. In this study, the potential distribution of Chinese sea buckthorn was found to be similar to its occurrence records (Fig 2B), which demonstrates that the distribution pattern of Chinese sea buckthorn is in equilibrium with climatic conditions. The results of MaxEnt model showed that little novel habitat outside the blue polygon as shown in the MESS map (Fig 6), which indicates that the MaxEnt model output had high reliability, because the habitat was interpolated, but not extrapolated, by the MaxEnt algorithm. Therefore, the results demonstrated that the occurrence records were sufficient to simulate the MaxEnt model, which was further shown by the excellent AUC performance with 0.93 (significantly better than random). In addition, the MaxEnt model is an open source and free software, and the workflow integrating MaxEnt and GIS tools used in this study are ideal for other species around the world.
The MaxEnt model shows that the potential distribution pattern of Chinese sea buckthorn (southwest to northeast orientation) is perpendicular to the East Asian monsoon sweep Climatic Factors Limiting the Distribution of Chinese Sea Buckthorn direction (southeast to northwest orientation) (Fig 7) [60]. We speculated that the distribution patterns were related to the northwest winter monsoon, the southwest summer monsoon, and the southeast summer monsoon systems in China. In winter, the Siberian cold winter monsoon blows from the northwest to southeast, which brings dry and cold air. In summer, the southeast summer monsoon blows from the Pacific Ocean and the southwest summer monsoon blows from the Indian Ocean, which brings moist and warm air [61]. In the northwest of potential distribution area of Chinese sea buckthorn, the plant species mainly locates in the center of the Loess Plateau, which belongs to inland area of China and is far away from East China Sea. Therefore, air mass of Eastern Asian monsoon arrives in Loess Plateau later and continues for shorter period, while that of Siberian cold winter monsoon reaches this area earlier and lasts for longer time within a year [62,63]. In the southwest of potential distribution area of Chinese sea buckthorn, the plant species mainly locates in the Hengduanshan Mountains, where the climate is mainly influenced by the southwest monsoon and Tibetan Plateau (the roof of the world) [60]. From Yun-Gui Plateau to Tibetan Plateau, there are many dry-hot valleys of north-south orientation across complex terrain conditions (Fig 1), which are the suitable habitats for Chinese sea buckthorn (Figs 1 and 7). When the southwest monsoon blows from Yun-Gui Plateau to Tibetan Plateau in summer, the warm and moist air mass can only sweeps northward through the north-south valleys. The distribution pattern of Chinese sea buckthorn is very similar to that of liaotung oak (Quercus wutaishannica), a species endemic and native to China [64]. However, it is worth mentioning that the south boundary of liaotung oak is not able to extend to the northern slope of the Qinling Mountains, which represent a dividing line between the warm temperate and northern subtropical zone in China. The exotic species black locust (Robinia pseudoacacia) could survive in the same region as Chinese sea buckthorn, but the distribution of black locust does not show the southwest-northeast direction [47] because black locust is mainly limited by low temperatures unlike those of the Eastern Asia monsoon climate. Accordingly, we speculated that other native species in this region are influenced by monsoon climate, but further verification is necessary.
Evaluation of climatic factors revealed that the distribution of Chinese sea buckthorn is mainly affected by PDM, AP, CI and ART. These factors are associated with the strong, cold Siberian winds and high altitudes of the Tibetan Plateau, which demonstrate that winter climate conditions play an important role in determining the potential distribution of Chinese sea buckthorn. Chinese sea buckthorn is a drought-tolerant plant that requires little water in winter ( Table 2 and Fig 5). The limiting factors mapping (Fig 8) shows that high rainfall is the limiting factor determining the southern boundary of the species, as shown in the response curves of PDM and AP. It is possible that excessive rainfall causes Chinese sea buckthorn to be out-competed by forest species, leading to the restriction of southern boundary. Low temperature is the limiting factor determining the northern boundary of this species, which is shown by the response curves of CI and ART. Low temperatures may limit plant physiological tolerance, leading to low survivability of Chinese sea buckthorn in the far northern areas.
A serious soil erosion problem exists within the range of potential distribution region of Chinese sea buckthorn. Chinese sea buckthorn is an effective nitrogen fixing plant with a welldeveloped root system, it has been widely planted for preventing soil erosion. At the same time, Chinese sea buckthorn is a plant species of high economic value (nutritious food, high medicine and firewood), and its potential distribution region is suitable for the construction of Chinese sea buckthorn farms, especially the core distribution area, which can improve the current cropping system and planting structure. We found that the climatic thresholds of the core area of Chinese sea buckthorn are 1.0-7.0 mm for PDM, 344.0-1089.0 mm for AP, -47.7-0.0°C for CI, and 26.1-45.0°C for ART, which were estimated based on the map of the climatically suitable habitats (Fig 2). Many studies have revealed that shifts in the climatic niche of plant species seldom occurred (niche conservative) [65][66][67]. Therefore, we can calculate the fitness of Chinese sea buckthorn under any given climatic conditions by using the climatic thresholds reported in Table 2 and the response curves shown in Fig 5 or S1 Fig. We could also use local weather station data to indicate when this species could be planted and cultivated.
It is worth mentioning that the predicted distribution map is based on statistical models MaxEnt, which does not currently consider soil factors and land use types. Thus, the predicted suitable area might partly be unsuitable for cultivation of the species (e.g. urban land, water body). This study mainly focused on the most suitable climatic habitat and the climatic limiting factors for Chinese sea buckthorn at a national scale with a 10 arc-min resolution. It is the first step of macro planning, and it is necessary to consider local site conditions when this predictive map is applied in practice. Future research in the predicted core suitable area should focus on experimental introduction and cultivation, breeding, selecting prior planting area, and integrating macro-planning of cropping systems and planting structure. Considering the potential economic and ecological value of Chinese sea buckthorn, the plant species would contribute greatly to controlling soil erosion and increasing the income of local farmers in the future.

Conclusions
Chinese sea buckthorn is a plant species of high economic value (nutritious food, high medicine and firewood) and a plant species of effective nitrogen fixing with a well-developed root system. It has been widely planted for preventing soil erosion. For scientific cultivation of this species, there is a need to identify suitable habitat where the species is most optimal for survival. To realize this target, this study collected 97 occurrence records of Chinese sea buckthorn from herbarium and publications, and used 13 climatic variables to simulate potential suitable habitat for Chinese sea buckthorn using MaxEnt and GIS. The simulated potential distribution area of Chinese sea buckthorn is mainly located in semi-humid and semi-arid regions and begins from southwest of the Hengduanshan Mountains and extends to the southern boundary of the Daxinganling Mountains, which occupies 11.7% of the land area of China and spans 12 provinces including Xizang, Yunnan, Sichuan, Gansu, Qinghai, Ningxia, Shaanxi, Shanxi, Hebei, Inner Mongolia, Henan, and Liaoning with the core areas mainly in Xizang, Sichuan, Gansu, Qinghai, Ningxia and Shaanxi provinces. The northern boundary of this species is mainly limited by extreme cold temperate (CI and ART), and the southern boundary of this species is mainly affected by abundant rainfall (AP), especially in winter (PDM). We can calculate the fitness of Chinese sea buckthorn under any given climatic conditions by using the climatic thresholds reported in Table 2 and the response curves shown in Fig 5 or S1 Fig. We can also use local weather station data to indicate when this species can be planted and cultivated. This study mainly focuses on the most suitable climatic habitat and the climatic limiting factors for Chinese sea buckthorn at a national scale with a 10 arc-min resolution. Thus, the predicted suitable habitat might partly be unsuitable for cultivation of the species. It is the first step of macro planning, and it is necessary to consider local site conditions when this predictive map is applied in practice. We believe that identification of priority areas for planting of Chinese sea buckthorn can be realized by integrating remote sensing data and boundary of suitable habitat of the species.