Habitat differentiation and conservation gap of Magnolia biondii, M. denudata, and M. sprengeri in China

The flower buds of Magnolia biondii, M. denudata, and M. sprengeri are the materials of Xinyi, a traditional Chinese medicine. The harvest of flower buds and habitat fragmentation caused by human disturbance heavily threatens the natural regeneration and survival of these three Magnolia species. With the aim to support the conservation and improve the effectiveness of conservation, we performed an assessment on habitat suitability, influences of environmental variables on habitat suitability, and the conservation gap of these three Magnolia species, based on the Maxent modeling method. The results indicated that: (1) altitude, annual mean temperature, extreme temperature, temperature fluctuation, annual precipitation, and extreme precipitation are the most influential environmental variables for the distribution of M. sprengeri, M. biondii, and M. denudata; (2) obvious habitat differentiations were observed among M. biondii, M. denudata, and M. sprengeri. M. sprengeri tends to be located in further northern areas with higher altitudes, lower temperatures, and lower precipitation compared to M. biondii and M. denudata; and (3) a large proportion of suitable habitats have been left without protection. Woodland and forest shared the largest area out of the suitable habitats. However, grassland, agricultural land, residential land, and mining and industry areas also occupied large areas of suitable habitats.


INTRODUCTION
Xinyi is a famous traditional Chinese medicine. It is made from the dried flower buds of Magnolia biondii Pampan., M. denudata Desr., and M. sprengeri Pampan. (Chinese Pharmacopoeia Commission, 2015). Xinyi was first recorded in the ancient book, "Sheng Nong's Herbal Classic", in which it was regarded as a high-grade medicine. Xinyi was known to be used as a medicinal material in the treatment of cold-headaches, nasal obstruction, and nasosinusitis (Wang, 2006;Huang & Yan, 2009;Ma, Zhang & Chen, 2017). In recent years, Xinyi has become more widely used in medicine, and its demand has greatly increased. However, overharvesting of flower buds exerts negative influences on the natural regeneration of these three Magnolia species (Jiang & Sheng, 1997). In addition disturbance within the suitable habitats. We hope that this research can provide technical support for the conservation and cultivation of M. biondii, M. denudata, and M. sprengeri.

Species distribution data
The location information of M. biondii, M. denudata, and M. sprengeri was derived from specimens preserved in the herbarium (http://www.cvh.ac.cn/) (Fig. 1). In the past, longitudes and latitudes were not recorded because of the lack of global positioning systems (GPS). However, the addresses where the specimens were collected were recorded. Those place names were used to infer longitude and latitude using Google Earth (https://www.google.com/earth/). Duplicate records were deleted from the dataset. In addition to this, records of specimens made from cultivated trees were excluded from this dataset. Finally, 162, 93, and 90 specimen records were used for the modeling of M. biondii, M. denudata, and M. sprengeri, respectively.
Topographic variables (altitude, slope, and aspect) were extracted from the ASTER Global Digital Elevation Model (GDEM). The ground resolution of the ASTER GDEM is 30 m, the vertical accuracy is estimated to be 20 m, and the horizontal accuracy is 30 m, with 95% confidence (Tachikawa et al., 2011). The ASTER GDEM was resampled and interpolated to the same spatial resolution as that of the bioclimatic variable (2.5′). Slope and aspect were calculated using the spatial analyst tools provided by ArcGIS (V9.3). For the aspect, the northern slope was set as 0 and increased gradually to 180 at the southern slope. This served to transform the aspect from a circular variable into a continuous variable.
In order to decrease the influences of high correlations between pairs of environmental variables on model fitting, we set 0.7 as the maximum correlation coefficient allowed between environmental variables (Dormann et al., 2013). Table S1 presents the Figure 1 The location where the specimens of Magnolia biondii, Magnolia denudata, and Magnolia sprengeri were collected. A-C are, respectively, the location map of Magnolia biondii, Magnolia denudata, and Magnolia sprengeri.
Full-size  DOI: 10.7717/peerj.6126/ fig-1 correlation coefficients between pairs of environmental variables. Altitude, slope, aspect, mean diurnal range, temperature seasonality, minimum temperature of coldest month, mean temperature of wettest quarter, annual precipitation, precipitation seasonality, and precipitation of warmest quarter were used to perform the model fitting.
Nature reserves and land use/cover data Vector data for national nature reserves were provided by the Nanjing Institute of Environmental Sciences, Ministry of Environmental Protection, China (http://www.nies.org/). Land use/cover data of China (2015) with one km spatial resolution were provided by the Resources and Environment Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/).

Model fitting
Maxent software was employed to assess the habitat suitability of M. biondii, M. denudata, and M. sprengeri in this research. Among the records of specimens, 25% were randomly selected for the test data and the remainder was used as training data to fit the model. The output format of the model values was cloglog, and the model gives an estimate between 0 and 1 for probability of presence. Other parameters were set as the following: "maximum iterations" was 500, "maximum number of background points" was 10,000, "replicates" was 1, and "replicated run type" was cross-validate. The "convergence threshold" was set as 0.00001, and the regularization multiplier was set as 1.
The Maxent software presented different approaches to assess the contribution of environmental variables: percent contribution, permutation importance, and jackknife test.
Percent contribution and permutation importance are based on the gain of the mode; gain is similar to deviance, in that it can be used to measure the goodness of model fitting (Phillips, 2017). Percent contribution represents the contribution that the environmental variable made to the fitted model, and depends on the path the model adopted to achieve the optimal solution (Songer et al., 2012;Smart et al., 2012). The permutation importance depends on the final model, and is denoted by the change of area under the receiver operating curve (AUC) due to the variation in the values of predictors between presence and background points (Songer et al., 2012). In addition, when high correlation exists between pairs of environmental variables, the percent contributions should be interpreted with caution (Phillips, 2017). The results of percent contribution and permutation contribution are shown in Table S2.
The jackknife test was conducted based on the training gain, test gain, and AUC. In the jackknife test, a model was first developed by using all the environmental variables. Then, each environmental variable was omitted in turn, and a model was created with the other environmental variables. Additionally, a model was developed using only one environmental variable. Through comparing the change in model gain under different model fitting strategies, the importance of each environmental variable was assessed. The result of the jackknife test is shown in Fig. S1.
Based on Table S2 and Fig. S1, we ranked these environmental variables according to the contribution of each variable (Tables S3-S5) for M. biondii, M. denudata, and M. sprengeri.
Response curves were also fitted to show how the habitat suitability changed when environmental variables varied from minimum to maximum. In order to diminish the influences of correlation among environmental variables, each response curve was fitted by using only one corresponding environmental variable, omitting other variables. The response curves of M. biondii, M. denudata, and M. sprengeri are presented in Figs. S2-S4, respectively. The suitable range of environmental variables of these three species (Table 2) were estimated based on the response curves of habitat suitability to environmental variables.

Model validation
In many studies, only presence data was available for species distribution modeling; no negative records (absence) could be used to measure specificity. It was hard to use the receiver operating characteristic curves (ROC) to assess the performance of the fitted models. Phillips, Anderson & Schapire (2006) presented an approach that could avoid this problem by distinguishing presence from random, rather than presence from absence. Therefore, ROCs were used to validate the fitted models in the Maxent software.
In ROCs, the y-axis indicated the sensitivity, and the x-axis indicated the false-positive fraction (1-specificity). The AUC was calculated as a measure of prediction success. The AUC ranged from 0 to 1, with a value higher than 0.5 indicating a better-than-random performance event (Fielding & Bell, 1997). The following system was often used to indicate the performance of the model (Swets, 1988): 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). In this research, the AUC of the training data and test data were both higher than 0.9 (Fig. S5). This indicated that the performance of the fitted models was excellent.

Threshold selection
The output ASC II grid is a continuous probability data that ranges from 0 to 1. Many thresholds were used to transform the continuous probability into binary data; those used include 0.5 (Waltari & Guralnick, 2009), 0.8 (Ramírez-Barahona et al., 2009, the minimum predicted value (Phillips, Anderson & Schapire, 2006), the 10th percentile training presence threshold (Brito et al., 2009), the 20th percentile training presence threshold (Donegan & Avendaño, 2010), the threshold that results in a sensitivity of 95% (Newbold et al., 2009), and maximization of the sum of sensitivity and specificity (Liu, White & Newell, 2013;Liu, Newell & White, 2016). Different methods were adopted (including the above methods) to estimate thresholds by the Maxent software and the thresholds estimated in this research were all lower than 0.5 (Tables S6-S8). Li et al. (2014) deemed that grids with probabilities higher than 0.6 could be regarded as the most suitable habitats (core area). Liu, White & Newell (2013) and Liu, Newell & White (2016) investigated which thresholds can be used confidently with presence-only data, and the thresholds estimated by most methods were lower than 0.8.
The threshold should be chosen based on the objective for generating the species distribution model (Wilson et al., 2005). In this research, the main objective was to determine the location of the most suitable habitats, and provide guidance for the conservation of these Magnolia species. The area of the suitable habitats cannot be too large for the consideration of conservation cost. Therefore, a large threshold should be selected to transform the continuous probability into binary data.
Based on the objective of this research and the thresholds used in previous studies, we selected 0.8 as the threshold to transform the continuous probability into binary data in this research. Grids with a probability higher than 0.8 were regarded as suitable habitats.

Assessment of conservation gap and human disturbance
The suitable habitats of M. biondii, M. denudata, and M. sprengeri were merged into a single suitable habitat layer by using the "merge" tool provided by ArcGIS (v9.3). The overlay analysis of nature reserves and suitable habitats was performed to assess the conservation gap of these species. Additionally, the overlayer analysis of land use and suitable habitats was conducted and the percent of each land use type within the suitable habitats was calculated to indicate the degree of human disturbance.

The importance of environmental variables
The contribution of environmental variables were presented in Table S2 and Fig. S1. Among these 10 selected variables, altitude is the most influential topographic variable on habitat suitability (Table S3), temperature seasonality and min temperature of coldest month were regarded as the most important temperature variables (Table S4), annual precipitation was recognized as the most influential precipitation variables (Table S5).
We could also see that high correlations existed between altitude and annual mean temperature, max temperature of warmest month, mean temperature of warmest quarter, temperature seasonality and temperature annual range, min temperature of coldest month and mean temperature of driest quarter (Table S1). Additionally, high correlations existed between the annual precipitation and precipitation of the wettest month, precipitation of the driest month, precipitation of the wettest quarter, and precipitation of the driest quarter (Table S1).
Based on these, we concluded that altitude, annual mean temperature, extreme temperature, temperature fluctuation, annual precipitation, and extreme precipitation are the most influential environmental variables for the distribution of these Magnolia species.
Habitat characteristics of M. biondii, M. denudata, and M. sprengeri Compared to M. biondii and M. denudata, M. sprengeri tended to grow in higher altitude areas; they could tolerate greater temperature fluctuations (annual range), lower extreme temperatures, less precipitation, and greater precipitation fluctuations (Table 2). Between M. biondii and M. denudata, M. denudata favors lower and warmer areas with smaller temperature fluctuations, more precipitation and smaller precipitation fluctuations (Table 2).

Spatial distribution of suitable habitats
The suitable habitats of M. denudata were mainly located in the middle sub-tropics (Fig. 2). Small areas of suitable habitats could also be found in the southeastern areas of Tibet and Taiwan. The spatial distribution patterns of M. biondii were similar to that of M. denudata. M. sprengeri extended to locations in further northern areas (Fig. 2), and the area of its suitable habitat was smaller than that of M. denudata and M. biondii.

Conservation gap and human disturbance
The overlay analysis of nature reserves and suitable habitats indicated that only 4.75% of the suitable habitats were under the protection of national reserves (Fig. S6). Large areas of suitable habitats were left without any protection.
The overlay analysis of land use and suitable habitats indicated that woodland and forest shared the largest area of the suitable habitats ( Fig. 3; Fig. S7). However, in the land use/cover data, orchards, mulberry fields, tea gardens, and other kinds of plantation forest were not differentiated from natural woodlands and forest. It was difficult to determine the percent of plantation forest in the woodland and forest categories. Agricultural land, residential land, mining land, and industrial areas also occupied a large portion of suitable habitats (Fig. 3). In addition to this, grassland is a very important land cover within the suitable habitats. However, in the southern area of China, most grasslands are secondary vegetation caused by deforestation. Heavy human disturbance still exists in these suitable habitats.  (Table 2; Fig. 2). This information is very important for us to perform in situ or ex situ conservation of these Magnolia species. Moreover, this information can assist in the selection of suitable areas to cultivate these species. However, a deficiency of this study is that only variables of topography, temperature, and precipitation were used to estimate the habitat suitability; other variables, such as soil, hydrology, and geology, were not included in the modeling. This might cause some uncertainties in the assessment of habitat suitability (Gaston, 2003). In addition, the geographical coordinates of the specimens used in the modeling were not measured by GPS. They were inferred from the recorded names of the villages and towns located near the sites where the specimens were collected. Villages and towns are usually located in flat areas with relative low altitudes. This might lead to mismatches between the altitude derived from the Digital Elevation Model (DEM), and the actual altitude of collection. Nature reserves are the most important and effective refuges for wild plants and wild germplasm resources. However, only a small percent of the suitable habitats are under the protection of national nature reserves (Fig. S6). Heavy human disturbance still exists in these suitable habitats (Fig. 3). This remind us that more measures and actions should be taken to conserve the wild resources of M. biondii, M. denudata, and M. sprengeri. However, in this research, only national nature reserves were used to estimate the conservation gap. Other nature reserves, such as provincial nature reserves, municipal nature reserves, county level nature reserves, national forest park, scenic spots, and geological parks were not included in this research. Therefore, the percent of protected suitable habitats could be higher than 4.75% (Fig. S6). The habitat suitabilities estimated by our research provide a scientific base for the selection of priority conservation areas of M. biondii, M. denudata, and M. sprengeri. Besides habitat suitabilities, genetic diversity should also be considered for the conservation of endangered species populations (Cires et al., 2013). Genetic diversity was regarded as the key to the survival of a species or population (Beardmore, 1983;Antonovics, 1984;Barrett & Kohn, 1991). It is essential for us to understand the genetic variation of inter-populations and genetic diversity of intra-populations of endangered species when we develop strategies for in situ and ex situ conservation activities (Hogbin & Peakall, 1999). Previous studies have found that the genetic diversity or genetic variation in a population is not correlated with the population size (Zhang et al., 2010). Therefore, genetic diversity within different sizes of populations and genetic variation among different populations should be fully assessed. More attention should be paid to the conservation of populations with high genetic diversity and large amounts of "private" fragments to prevent further loss of genetic diversity (Zhang et al., 2010).
China is the world center of distribution for Magnoliaceae; more than 40% of the species originated from southern China (Wang & Jiang, 2001;Cicuzza, Newton & Oldfield, 2007). In past decades, Magnoliaceae species in China have suffered heavy human disturbance, such as the harvesting of timber and medicinal material, as well as deforestation (Wang & Jiang, 2001). Many taxa of Magnoliaceae are threatened due to overexploitation and habitat destruction (Cicuzza, Newton & Oldfield, 2007;Cires et al., 2013). According to the "China Biodiversity Red List," 71.7% of Magnoliaceae species are classified as "Threatened" (Ministry of Ecology and Environment of China and Chinese Academy of Sciences, 2013). Therefore, more studies should be performed on more threatened species; field surveys of threatened species and assessments of habitat suitability and genetic diversity are required to support the ex situ and in situ conservation of Magnoliaceae species.

CONCLUSIONS
(1) Altitude, annual mean temperature, extreme temperature, temperature fluctuation, annual precipitation, and extreme precipitation are the most influential environmental variables for the distribution of M. sprengeri, M. biondii, and M. denudata.
(2) Obvious habitat differentiation was found among M. biondii, M. denudata, and M. sprengeri. Compared to M. biondii and M. denudata, M. sprengeri tended to be located in further northern areas with higher altitudes, lower temperatures, and lower precipitation.
(3) Large proportions of suitable habitats were left without any protection. Heavy levels of human interference still exist in the suitable habitats of these three Magnolia species.