Species Distribution Modeling of Sassafras Tzumu and Implications for Forest Management

: Sassafras tzumu (Chinese sassafras) is an economically and ecologically important deciduous tree species. Over the past few decades, increasing market demands and unprecedented human activity in its natural habitat have created new threats to this species. Nonetheless, the distribution of its habitat and the crucial environmental parameters that determine the habitat suitability remain largely unclear. The present study modeled the current and future geographical distribution of S. tzumu by maximum entropy (MAXENT) and genetic algorithm for rule set prediction (GARP). The value of area under the receiver operating characteristic curve (AUC), Kappa, and true skill statistic (TSS) of MAXENT was significantly higher than that of GARP, indicating that MAXENT performed better. Temperate and subtropical regions of eastern China where the species had been recorded was suitable for growth of S. tzumu . Relative humidity (26.2% of permutation importance), average temperature during the driest quarter (16.6%), annual precipitation (12.6%), and mean diurnal temperature range (10.3%) were identified as the primary factors that accounted for the present distribution of S. tzumu in China. Under the climate change scenario, both algorithms predicted that range of suitable habitat will expand geographically to northwest. Our results may be adopted for guiding the preservation of S. tzumu through identifying the habitats susceptible to climate change.


Introduction
Understanding the relationship between the geographical distribution of species and that of environmental characteristics is a central issue in ecology [1,2]. Changes in environmental characteristics, especially climate, have significant implications for life-history traits, flowering and fruiting phenology, as well as habitat requirements that have consequently altered the forest ecosystem functions and structures [2]. Indeed, the changes in global climate that have been detected in the last several decades have generated numerous variations in species distributions, which may serve to be the primary cause leading to species extinction in the short run, based on the direct or synergistic effect with additional extinction enablers [3,4]. Therefore, understanding the distribution of habitat (with suitable climate) can enable forest managers to accurately assess ecosystem and species vulnerabilities in the face of climate change [5]. However, a lack of appropriate biogeographical data accounts for a leading drawback that hinders the implementation of appropriate conservation measures [6].
When an insufficient amount of distribution data is available, species distribution models (SDMs) are known to be the efficient approach for evaluating the possible geographical distribution of species; these models have been used in various fields, such as biogeography, ecology, conservation, and evolutionary biology [7][8][9][10][11][12]. There are various SDMs, such as the generalized additive model, generalized linear model [9,13], general rule set prediction (GARP) [14], DOMAIN [9], bioclimatic envelopes [9], and MAXENT (maximum entropy) [15]. These models use different species records (presence/pseudo-absence, presence/absence, or presence-only data) to predict the distribution of various species.
Recent integrated comparisons of presence-only modelling approaches have shown that GARP and MAXENT perform well in prediction when compared with other approaches [16][17][18]. MAXENT has been used to estimate species distributions through examining the possible maximum entropy distribution [15]. GARP produces mathematical rules in the random set after an iterative rule selection process, by testing, incorporating, and rejecting possible outcomes [10]. Such rule sets have been integrated randomly for generating the possible species niche restrictions based on environmental conditions [7]. Evaluating the performance of GARP and MAXENT will help to reveal the variations in their ability to accurately predict the potential distribution of species.
Sassafras tzumu (Chinese sassafras) is an endemic deciduous tree of China belonging to the Lauraceae family. This species grows well in moist, well-drained, or sandy loam soils and tolerates a variety of soil types [19]. It occurs in Anhui, Fujian, Guangdong, Guangxi, Guizhou, Hubei, Hunan, Jiangsu, Sichuan, Yunnan, and Zhejiang Provinces of China [19], and grows in both broad-leaved evergreen and in mixed deciduous forests with an elevation of 100-1900 m [19]. This species plays a significant role in increasing the ability of soil to retain water and reducing the loss of soil to surface runoff [20]. In addition, the roots and bark of S. tzumu are often harvested from the wild and the species is used as a medicinal plant locally for treating inflammation and traumatic injuries [21]. The essential oils obtained from its bark, leaf, and root vary significantly in chemical composition, which suggests that they might vary in their pharmacological effects as well [22]. Moreover, this speciesimportant timber tree is valued for its decorative and often fragrant wood. This wood is strong, extremely durable, handsome, and maintains its shape well [19]. The profit of one hectare of S. tzumu timber is ca. 250,000 RMB [20].
Over the past few decades, increasing market demands and unprecedented man-made activity in its natural habit have created new threats to S. tzumu [22]. Nonetheless, the potential geographic distribution of S. tzumu has rarely been investigated in existing studies. Determining the effects of climate changes on habitat suitability is also a critical problem linked to its economic value and ecological significance. In this study, the MAXENT and GARP models were adopted for predicting the possible distribution of S. tzumu. Specifically, this study aimed to (1) model this species' potential geographical distribution; (2) identify the critical environmental factors shaping the distribution of S. tzumu; and (3) discuss the variations of appropriate habitat distributions in the face of climate changes. Our results will allow researchers to identify the future suitable habitat and help in the use, management, and cultivation of S. tzumu.

Species Occurrence Data
Occurrence data of S. tzumu was acquired from various sources, such as the Chinese virtual herbarium [23], the Global Biodiversity Information Facility [24], and our field survey, as well as literature reports. The longitude and latitude of few specimens (ca. 20) that give an exact village location were found using Google Earth. To reducing the possible errors of specific points, the locations that were duplicated or lack detailed information were deleted. Eventually, a total of 202 records were screened and used to model the distribution of S. tzumu ( Figure 1).

Environmental Parameters
In total, 28 environmental parameters possibly affecting S. tzumu distribution were selected. These include 19 bioclimatic variables having a 30-second spatial resolution obtained from the World Climate Database [25], three topographic variables, i.e. slope degree, aspect, and altitude, were extracted from the digital elevation model (DEM, 20 × 20 m resolution) obtained from Geospatial Data Cloud [26], and three soil variables, i.e. the type, pH, and organic carbon content of soil which were acquired from the Center for Sustainability and the Global Environment [27]. In addition, data related to the normalized difference of vegetation index (NDVI), relative humidity, and sunshine duration in growing season were acquired based on the China Meteorological Data Sharing Service System [28]. The annual average NDVI was calculated based on a 10-year time-series of SPOT vegetation imagery with a 1 km 2 resolution; then it was used to measure vegetation cover [28].
The Beijing Climate Center Climate System Model version 1.1 (BCC-CSM1.1) were obtained based on representative concentration pathways (RCPs) 2.6 and 8.5 for 2050 and 2070, issued via the Intergovernmental Panel on Climate Change (IPCC) Assessment Report 5 (AR5), which were used as the simulation data for future distribution. The BCC-CSM1.1 has been recommended for use in studies on climate change and on the operational short-term forecast of climate across China [29]. The RCP 2.6 reflects potential radiative forcing by 2100 compared with the pre-industrial values of +2.6 W/m 2 which is optimistic, while RCP 8.5, a more pessimistic scenario, reflects high emission levels of greenhouse gases, and leads to 8.5 W/m 2 of radiative forcing in the year of 2100 [25]. The other variables were not changed to analyze the SDM in the projected future climate scenarios.
Each bioclimatic layer had a resolution of 1 km, and each environmental layer was processed using a system with identical projection, cell size, and extent (WGS84 Longitude-Latitude projection) in ArcGIS 10.1 (ESRI, Redlands, CA, USA). To eliminate the multicollinearity among variables, the Pearson correlation coefficient were used. Only one parameter was selected for those with high crosscorrelation (r 2 > 0.90) [30]. The decision of which variable was kept were based on principal component analyses (For more detail see Yi et al. [30]). Eventually, the number of predicting factors was decreased to 16 (Table 1).

GARP Analysis
The GARP 1.1 was adopted for predicting species distribution [31]. The occurrence records were classified into two groups, among which 75% data were randomly screened in training the model, while the remaining 25% were used to test the model. Afterwards, the environmental variables shown below were adopted to fit the GARP model, including runs of 100 times, a convergence limit of 0.01, and a maximum of 1000 iterations. Additionally, four types of rules (range, atomic, logistic regression, and negated range) were used, while the remaining variables were set as defaults. Eventually, the 'best subset' models were chosen for result mapping. Then, the ArcGIS 10.1 spatial analysis tools were adopted to average the top ten models for generating one singular distribution map.

MAXENT Analysis
MAXENT 3.3.3k [18] was used in analysis. Similar to GARP, we also used 25% of occurrence records to test the model and 75% of occurrence records for training. A bias file layer was used to avoid bias of sampling within those species' occurrence records [32]. This file was generated using occurrence point data by deriving a Gaussian kernel density map that was rescaled from 1 to 20 based on Elith et al. [33]. Altogether, four possibilities, including bioclimatic rules, logistic regression, negated range rules, and range rules, were applied. The algorithm was run until convergence or for 1,000 iterations. The internal jackknife of MAXENT was also adopted for testing and assessing the significance for all environmental parameters in the prediction of the distribution of S. tzumu. Other optimization parameters are selected by default.

Environmental Parameters That Determined S. tzumu Distribution
Maxent can be used in several ways to quantify how each variable contributes to the model. Here, 'permutation importance' was used to identify the most important bioclimatic variables that predicted the geographical distribution of a specific group. Permutation importance measures the decrease in the training area under the receiver operating characteristic curve (AUC) that results from randomly permuting the value of a specific variable when training the model. The decrease in the training AUC is proportional to the importance of each variable to the model [18].

Shift of The Distribution Center of Suitable Habitat
The distribution center of S. tzumu under current and future climate change scenarios were calculated using the SDM toolbox [34] and the primary shifts of the suitable habitat of S. tzumu were analyzed. The movement of the distribution center was tracked to examine the distributional shift under climate change.

Model Assessments
Three approaches were used for model evaluation, including AUC [35], Cohen's Kappa, as well as the true skill statistic (TSS) [36]. To compare the model performance of GARP and MAXENT, a Mann-Whitney U test were used. The AUC has been an approach used to measure model performance, which is independent of a threshold and has been universally accepted [35]. The AUC value can range from 0-1, and AUC values close to 0.9 are recognized as favorable performance. The value of Cohen's Kappa was between −1 and + 1, in which + 1 suggests excellent performance, while values of ≤ 0 indicate that a performance that was not superior to a random result [15]. The value of TSS also falls within the interval of −1 to + 1, among which ≤ 0 suggests poor model performance, while 1 predicts the best possible performance. The R software ''ROCR'' package was adopted for calculating the TSS, kappa and AUC values of those ten models in each algorithm.

Model Accuracy and Comparison
The test AUC, Kappa, and TSS value of both MAXENT and GARP were > 0.85, indicating the good performance of both models. According to the results from a Mann-Whitney U test, MAXENT had significantly higher AUC, Kappa and TSS scores than those of GARP ( Table 2), indicating that MAXENT models performed better than GARP models.

Regions for Possible S. Tzumu Distribution
Based on the prediction results of MAXENT and GARP models (Figure 2), the climate in most regions of southeastern China, such as Anhui, Chongqing, Guangdong, Guangxi, Guizhou, Hubei, Hunan, Fujian, Jiangxi, and Zhejiang, is suitable for the growth of S. tzumu. However, the current modeled distributions via both algorithms are inconsistent: GARP forecast that large areas in Jiangsu, Anhui, Shaanxi, and Hainan were suitable, which were forecasted to be marginally suitable via MAXENT. Moreover, GARP predicted that the areas of potential distribution with high suitability was continuous and covers a large area, whereas those predicted by MAXENT were scattered and small.

Environmental Parameters that Determined S. Tzumu Distribution
Environmental parameters with the greatest contributions to constructing the possible S. tzumu distribution model based on the MAXENT algorithm included relative humidity (Hum, 26.2% contribution), mean temperature at the driest quarter (Bio 9, 16.6%), annual precipitation (Bio 12, 12.6%), and the mean diurnal temperature range (Bio 2, average of (max temp-min temp) every month) (10.3%). For those four environmental parameters, the cumulative contribution reached up to 65.7%, proving that these four parameters would make vital contributions to the possible future distribution of S. tzumu. According to those standard curves regarding the S. tzumu distribution model (Figure 3), we acquired the thresholds for leading parameters (with the probability of occurrence of > 0.3), including Hum of 75%-83%, Bio 9 of 4-20 °C, Bio 12 of > 1000 mm, and Bio 2 of 70%-90%.

Variations of Future Habitat Under Climate Change
In the RCP 2.6-2050 climate scenario, MAXENT forecasted that the gains in suitable habitat areas of S. tzumu would be in the central Sichuan, southern Anhui, southern Shaanxi, northern Hubei, and south Jiangsu (Figure 4 A1), which amounted to ca. 1.14 × 10 5 Km 2 (Table 3). In addition, losses of suitable habitat area of ca. 1.12 × 10 5 km 2 were predicted for Fujian and northern Guangdong. In the RCP 2.6-2070 climate scenario, MAXENT forecast that area of expanded distribution increase to 1.18 × 10 5 (Table 3), while the decreased habitat continued to decrease and accounted for ca. 0.98 × 10 5 km 2 ( Table 3). In the RCP 2.6-2050 and RCP 2.6-2070 climate scenarios, GARP predicted losses of appropriate habitat area accounting for ca. 1.38 × 10 5 km 2 and 1.76 × 10 5 km 2 , respectively, in Guangdong and Guangxi (Figure 4 A1, Figure 4 B2). Meanwhile, the areas of increased suitable habitat accounting for ca. 0.89 × 10 5 km 2 and 0.27 × 10 5 km 2 (Table 3) under RCP 2.6-2050 and RCP 2.6-2070 climate scenarios, respectively, were predicted for areas in southern Shanxi, and southern Anhui (Figure 4 A1, Figure 4 B2). Overall, under RCP 2.6 scenario, MAXENT predict an increase in distribution while GARP predict a decrease in distribution. Nonetheless, both MAXENT and GARP predicted the distribution shift to northwest ( Figure 5).  Under the RCP 8.5 climate scenario, both MAXENT and GARP predicted that the area of change was located in the same provinces as were predicted with RCP 2.6 climate scenario ( Figure 6) and showed a similar distribution shift ( Figure 5). MAXENT predicted that the increased area accounted for ca. 1.19 × 10 5 km 2 and 2.36 × 10 5 km 2 under the RCP8.5-2050 and RCP 8.5-2070 climate scenarios, respectively, while GARP predicted that the increased area accounted for 0.95 × 10 5 km 2 and 0.33 × 10 5 km 2 under the RCP 8.5-2050 and RCP 8.5-2070 climate scenario, respectively. MAXENT predicted that the decreased area accounted for ca. 1.09 × 10 5 km 2 and 0.68 × 10 5 km 2 under the RCP8.5-2050 and RCP 8.5-2070 climate scenarios, respectively, while GARP predicted that the decreased area accounted for 1.49 × 10 5 km 2 and 3.58 × 10 5 km 2 under the RCP 8.5-2050 and RCP 8.5-2070 climate scenarios, respectively.

Discussion
Understanding the relationship between the geographical distribution of a species and that of environmental characteristics needed by a species is a central issue in ecology. S. tzumu is an ecologically and economically important tree species. However, because S. tzumu has a weak natural ability to regenerate and because unprecedented damage has occurred in its natural habitat, wild germplasm resources are being exhausted. In this study, we used MAXENT and GARP to model the potential distribution of this species and the effects of climate change on its habitat suitability. Our results indicate that MAXENT performed better than GARP. Relative humidity, average temperature during the driest quarter, annual precipitation, and mean diurnal temperature range were identified as the primary factors that accounted for the present distribution of S. tzumu in China. Under RCP 2.6 and 8.5 climate change scenarios, both algorithms predicted that the range of suitable habitat will expand geographically to the northwest. SDM has been recognized as an efficient and universal approach to guide the management of forest in the face of global climate change [8]. However, model performance differs among different SDMs based on the sampling size, the study area, and the species modelled [9,16,[37][38][39]. Using the performance measures employed in this study, the results show that the MAXENT has superior performance over GARP and thus can be used to predict the geographical distributions of species, or of other species of similar ecological traits. Such findings conformed to those obtained from Elith et al. [9], Hernandez et al. [16], and Elith and Graham [40], who rated MAXENT higher than GARP.
MAXENT is advantageous when compared with additional ecological niche models for the prediction of possible species distributions. Specifically, MAXENT allows for model to be run based on a small number of sample locations [9,41]. As indicated in the present study, GARP predicted similar appropriate areas of distribution to those predicted by MAXENT, and both algorithms drew maps consistent with the known distribution of the species. However, GARP predicted extensive areas compared with those by MAXENT. These findings may have occurred because GARP and MAXENT have basic differences; GARP tends to result in models with a greater number of errors of commission than MAXENT, i.e., it would predict broader areas of suitable habitat [42]. Apparent errors of commission were derived from regions with potentially suitable habitat where the species was predicted to be present, but where presence cannot be demonstrated because the species has not been verified to exist there. This lack of verification may have several causes such as the existence of map pixels that indicate presence that lack documentation of the species simply because biologists have not adequately sampled an area [40,42].
Determining which environment factor is shaping and maintaining a species geographical distribution is a critical issue in ecology and evolution. According to our findings, the most important ones that explained the species' environmental requirements best were relative humidity, mean temperature during the driest quarter, the average diurnal range of temperature, and annual precipitation. In fact, each area possessing a relatively large probability for the occurrence S. tzumu was associated with rainy seasons having high humidity during the growing season; this might offer an appropriate growth environment for S. tzumu. In the meantime, those areas that had a reduced S. tzumu occurrence probability might experience a dry growing season. The above findings indicated the distribution of S. tzumu would occur within habitats with a high relative humidity.
The tolerance of a particular range of temperatures is one of the most important features used to explain the latitudinal distribution of various species [43]. S. tzumu generally grows in warm and humid regions with an average winter (December-February) temperature > 4°C. This finding agrees with the known climatic preference of S. tzumu [44]. Variations in temperature affect the distribution S. tzumu through affecting germination, water absorption, photosynthesis, transpiration, respiration, reproduction, and growth. Low winter temperature has been suggested to affect the dormancy breaking of S. tzumu seeds [45]. Moreover, based on the mean air temperature of winter in field records (occurrence points), S. tzumu occurrence was not detected within areas with means of <−1°C. The response curve of annual precipitation showed that the optimal range for S. tzumu was from >1000 mm. Patterns and annual amounts of precipitation serve as important factors in plant regeneration and survival for this species as well as in other ecosystem functions [46]. Previous studies have shown that plants of S. tzumu have a plastic response of seedling growth [47], biomass allocation [48,49], physiological characteristics [22,48], and phenology [50] to variations in annual precipitation.
The above parameters might affect the ultimate S. tzumu distribution and its ability to adapt ecologically. However, the species distribution models usually depict the fundamental niches of species, but not the realized niches [51]. Numerous parameters that affect those realized niche dimensions are usually not considered when predicting the possible species geographical distribution [52]. In our results, some parameters that may affect the distribution of S. tzumu were not considered because of a lack of robust data, such as interspecific competition, the effects of natural enemy predation, anthropogenic influence, and geographical barriers. Therefore, future studies need to incorporate these factors in their analysis [53]. Moreover, each bioclimatic layer used in our studied had a resolution of 1 km, predictions based on low resolution data are likely to overestimate the spatial extent of a species' habitat [14,17], and therefore, a higher resolution is needed in further studies.
At present, climate change has received unprecedented attention [54,55]. Although the future potential distributions of S. tzumu predicted by the two algorithms were inconsistent, both models predicted that the range of habitat suitable will expand geographically to the northwest. Typically, southern Shaanxi, and Anhui Province has been recognized to be inappropriate for this species under the present situations based on the literature [19] as well as by our model constructed based on present climatic conditions.
Climate changes can lead to the local or regional disappearance of S. tzumu and loss of entire ecosystems; alternatively, some ecosystems may be replaced by additional types of ecosystems [56]. The variations of temperature and precipitation regimes, such as the duration of the dry season, can lead to S. tzumu trees experiencing phenological shifts, which indirectly affects the faunal and floral species that depend on S. tzumu. In addition, climate changes may have adverse effects on a number of insects, mammals, and terrestrial birds that are indirectly and directly dependent on S. tzumu seeds, fruits, and flowers [57].
Our model results may be tailored to suit the preservation guidelines of S. tzumu through recognizing those critically susceptible and climatically appropriate habitats in which artificial regeneration must be considered during reforestation. Using our models, the appropriate climatic habitats for S. tzumu were detected to shift northwards outside the current area of natural distribution. However, the speed of natural migration is low; therefore, assisted migration of S. tzumu can serve as a potential preservation strategy if extensive future climate change creates new areas of climatically appropriate habitat [58]. Therefore, plantations of S. tzumu can be proposed within newly climatically appropriate habitats; besides, more attention should be paid to preserve the natural regeneration for such trees at sites determined to be in high-risk areas under future climate scenarios. Moreover, the non-change and gain within the climate space shown among the various S. tzumu ecoregions may be recognized to be the underlying refugia of climate change. S. tzumu trees currently grow under various climates, providing a way for the species to adapt to the new climates within specific natural habitats [59]. Consequently, the persistence of S. tzumu can be improved by the application of phenotypic plasticity and the selection of genotypes that potentially adapt to future climate scenarios.
To promote the sustainable development of S. tzumu and the related forest ecosystems, we recommend the following strategies: (1) in-situ and ex-situ conservation methods for wild S. tzumu and related information networks should be developed; (2) anthropogenic disturbance should be minimized; (3) a database for long-term and systematic monitoring of the population dynamics of S. tzumu should be established; (4) for high risk areas under climate change, forest managers should introduce other species that were known to be appropriate for specific climatic situations, rather than continue to create new plantations of S. tzumu; and (5) we recommend forest management agencies be given additional funding so they can conduct intensive surveys to analyze and understand the present situation related to wild S. tzumu. Land managers would then need to develop a management strategy and practical measures to conserve valuable wild S. tzumu resources.

Conclusions
Based on our findings as well as on previous biological data, the distribution of S. tzumu are suggested to be primarily driven by relative humidity, mean temperature during the driest quarter, mean diurnal temperature range, and annual precipitation. The temperate and subtropical regions of eastern China where the species has been recorded were also identified as having high environmental suitability for this species. Under RCP 2.6 and 8.5 climate change scenarios, both algorithms predicted that the spatial extent of the potential distribution would move to the northwest. Additional studies using additional environmental factors, such as anthropogenic influence and geographical barriers, as well as a higher resolution of bioclimatic layer will be needed to validate the effects of climate on the distribution of S. tzumu. Nonetheless, the maps produced in our study allows the quantitative assessment of risks related to regional climate change that could have an effect on S. tzumu cultivation; therefore, the present study offers valuable data related to forest management planning efforts in the face of global climate change. In addition, the findings in this work may also be used to optimize future field surveys and the detection of unknown populations.

Conflicts of Interest:
The authors declare no conflict of interest.