According to the Flora of China (Volume 30(1)) (Liu 1996), M. officinalis subsp. officinalis is mainly distributed in the northern and central subtropical regions, which includes the southern area of Shaanxi province, the south-eastern areas of Gansu province and Henan province, the western area of Hubei province, the south-western area of Hunan province, the central and eastern areas of Sichuan province and the north-eastern area of Guizhou province (Fig. 1). Magnolia o. subsp. biloba is mainly distributed in the central subtropical regions, which includes the western areas of Anhui province and Zhejiang province, the southern areas of Jiangxi province, Fujian province and Hunan province, the northern area of Guangdong province and the northern and north-eastern areas of Guangxi province (Fig. 1).

A total of 241 specimens were reviewed from the Chinese Herbarium (http://www.cvh.ac.cn/) to retrieve the distribution information of M. o. subsp. officinalis and M. o. subsp. biloba. As most specimens were identified as M. o. in the original background information of these specimens, we invited plant taxonomists to re-identify these specimens to the level of subspecies.

We also collected geographical distribution information for both subspecies from the published literature (Zhang et al. 2009; Yu et al. 2010; Wang et al. 2016; Yan et al. 2016). Specifically, specimen records of cultivated trees were excluded from the dataset. Finally, 163 records of M. o. subsp. officinalis and 106 records of M. o. subsp. biloba were used to estimate the habitat suitability, respectively (Fig. 1).

Figure 1.

The distributions of Magnolia officinalis subsp. officinalis and M. officinalis subsp. biloba in China.

Information on nineteen bioclimatic variables was retrieved from WorldClim (http://www.worldclim.org/) (Table 2). We used the averages for each of the years from 1970 to 2000 and were of spatial resolution 2.5' (Hijmans et al. 2005).

The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) was used to derive topographic variables in this research. The vertical accuracy of the ASTER GDEM is 20 m and its horizontal accuracy is 30 m with 95% confidence (ASTER GDEM Validation Team 2011). We re-interpolated the ASTER GDEM to the same resolution as the bioclimatic variables (2.5') using a spline function.

Topographic variables used in this research include altitude, slope and aspect. Slope and aspect were calculated using the Spatial Analyst Tools in ArcGIS (v9.3). Slope helps to identify the rate of maximum change in z-value from each cell and the range of slope values in degrees is 0 to 90 (ESRI 2008). Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbours and it can be thought of as the slope direction (ESRI 2008). The output raster value of aspect is the compass direction, in which aspect is expressed in positive degrees from 0 to 360, measured clockwise from north, east, south and again to the north. In this research, we defined the north-facing slope as 0° and the south-facing slope as 180°. As the aspect changes from north to south, it increases gradually from 0 to 180°. We used this technique to transform aspect values from a circular variable (0–360°) into a continuous variable (0–180°).

To reduce the deleterious effects of collinearity on model fit, the maximum coefficient allowed between pairs of variables was set to 0.7 (Dormann et al. 2013). Correlation coefficients amongst environmental variables are presented in Suppl. material 1: Table S1. We then selected the following ten environmental variables for use in fitting the model: altitude, slope, aspect, mean diurnal range, temperature seasonality, min. temperature of coldest month, mean temperature of wettest quarter, annual precipitation, precipitation seasonality and precipitation of warmest quarter.

We retrieved data on the spatial distribution of national nature reserves from the Nanjing Institute of Environmental Sciences, Ministry of Ecology and Environment of China (http://www.nies.org/). We also retrieved land cover data (2015) of 1-km spatial resolution from the Resources and Environment Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/).

Many software algorithms can be used to calculate habitat suitability using only species presence data, such as Bioclim (Busby 1986), Domain (Carpenter et al. 1993), Garp (Stockwell 1999), NFA (Hirzel et al. 2002) and Maxent (Phillips et al. 2006). The Maxent is one of the machine learning algorithm models and it is one of the most widely used and well-predicted species distribution models (Phillips et al. 2006; Elith et al. 2006; Harte et al. 2008; Wisz et al. 2008; Phillips et al. 2009; Yackulic et al. 2013).

In this study, we used Maxent to evaluate habitat suitability in M. o. subsp. officinalis and M. o. subsp. biloba. We randomly selected 25% of the data entries for each subspecies for use as test data, with the remaining 75% being used to train the model (training data). The Maxent model parameters were set as follows: ‘maximum iterations’ was 500, ‘maximum number of background points’ was 10000, ‘replicates’ was 1 and ‘replicated run type’ was ‘cross-validate’. The ‘convergence threshold’ was 0.00001 and the ‘regularization multiplier’ was 1.

The output ASCII grid produced by Maxent is continuous probability data ranging from 0 to 1, which represents the habitat suitability of a species in a specified region. Based on previous research (Liu et al. 2013; Li et al. 2014; Liu et al. 2016), we adopted 0.7 as the threshold to transform the probability data into binary data (0/1). Grids, with probability higher than 0.7, were defined as suitable habitats.

Percent contribution and permutation importance are approaches available in the Maxent software which evaluate the contribution of variables to model predictions. Percent contribution and permutation importance are estimated based on the model gain, which is closely related to the deviance and is used to measure the goodness of fit of the model (Phillips 2017). Percent contribution depends on the path that the model used to reach the optimal solution (Songer et al. 2012; Smarter et al. 2012), whereas permutation importance depends on the final iteration of the fitted model. It is expressed as the differences between the AUC (area under the receiver operating curve) values of presence and background points caused by variation in the predictor variables (Songer et al. 2012). It is worth noting that when pairs of variables are strongly correlated with one another, the percent contributions of each should be interpreted with great care (Phillips 2017). Therefore, in this study, we used permutation importance to assess the contribution of each environmental variable to the fit of the produced models.

We also fitted response curves to evaluate how habitat suitability responded to variation in the environmental variables. In order to reduce the effects of correlation between pairs of environmental variables on model fit, each response curve was fitted using only one environmental variable. The fitted response curves are shown in Suppl. material 2: Figs S1–S6. Based on these response curves, we identified the characteristics of suitable habitats based on their environmental variables (Table 1).

In numerous studies, only presence data (such as herbarium data) was available when attempting to model species distributions and habitat suitability. Absence data (negative records) are rare, despite their usefulness in assessing model specificity. Thus, it was difficult for us to use ROC (receiver operating characteristic curves) to evaluate the performance of the fitted models. Phillips et al. (2006) suggested an alternative approach which could avoid this issue through the distinguishing of presence from random, rather than presence from absence. In this way, ROC could be used to evaluate the performance of the fitted model.

In ROC, the ordinate axis represents the sensitivity, while the abscissa axis represents the false-positive fraction (1-specificity). The AUC (area under curve) is then used to measure the prediction success of the fitted model. AUC values range from 0 to 1; if AUC values are higher than 0.5, the imitative effect of the fitted model is deemed different from random (Fielding and Bell 1997). AUC values were then used to categorise model fits as follows (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 study, the AUC calculated, based on the training data and test data, were both higher than 0.9 (Suppl. material 2: Fig. S7). This indicated that the fitted models performed excellently.

We used overlay analysis and maps of national nature reserves and the identified suitable habitats to identify conservation areas where M. o. subsp. officinalis and M. o. subsp. biloba populations were likely to occur. We also used overlay analysis of land cover and the identified suitable habitats to calculate the percentage for each land cover type in the suitable habitats and evaluated the extent to which these habitats have been disturbed by human activity.

Amongst the ten selected variables, the environmental variables which most affected the distribution of M. o. subsp. officinalis were min. temperature of coldest month, followed by temperature seasonality and altitude (Table 1). Given that many of the environmental variables were highly correlated with one another (Suppl. material 1: Table S1), we can reasonably conclude that extremely low temperatures (min. temperature of coldest month, mean temperature of driest quarter, mean temperature of coldest quarter) and temperature fluctuation (temperature seasonality and temperature annual range) are the most important variables which influence the distribution of M. o. subsp. officinalis. Both altitude and high temperature (max. temperature of warmest month, mean temperature of warmest quarter) also significantly influence the distribution of M. o. subsp. officinalis (Table 1), but to a lesser extent than extreme low temperatures and temperature fluctuations.

M. o. subsp. biloba seems to be more influenced by annual precipitation and extreme precipitation (precipitation of wettest month, precipitation of driest month, precipitation of wettest quarter, precipitation of driest quarter, precipitation of coldest quarter). Temperature fluctuation (temperature seasonality, temperature annual range, isothermality and mean diurnal range) also had notable influences on the distribution of M. o. subsp. biloba.

Habitats suitable for M. officinalis subsp. biloba are mainly located in central subtropical regions (Fig. 2). Small areas of suitable habitats could also be found in the northern areas of Taiwan. Magnolia o. subsp. officinalis habitats extend to locations in areas further north (Fig. 2), i.e. north subtropical areas. Habitats suitable for M. o. subsp. officinalis were also found in the south-eastern area of Tibet.

Compared to M. officinalis subsp. biloba, M. officinalis subsp. officinalis tends to grow at higher altitudes (Fig. 2) and is able to survive in areas experiencing greater fluctuations in ambient temperature, lower extreme temperatures, less precipitation and greater fluctuations in precipitation (Table 2).

Figure 2.

The spatial distributions of habitats suitable for Magnolia officinalis subsp. officinalis and M. officinalis subsp. biloba.

Overlay analysis of the identified suitable habitats and national nature reserves found that only 8.4% of the habitats suitable for M. officinalis subsp. officinalis were located in national nature reserves. Similarly, for M. o. subsp. biloba, only 3.4% of suitable habitats were located in national nature reserves. Thus, large areas of habitats, suitable for both subspecies, are not protected and are thus at risk (Fig. 2).

Overlay analysis of the identified suitable habitats and land cover found that the majority of land in habitats, suitable for both M. o. subsp. officinalis and M. o. subsp. biloba, is woodland and forest (based on area) (Figs 3, 4). Suitable habitat areas were also constituted largely by agricultural land, residential land and mining and industry land (Figs 3, 4). In addition, suitable habitats were found to contain significant areas of grassland.

Figure 3.

Land cover status (2015) within the suitable habitats of Magnolia officinalis subsp. officinalis.

Figure 4.

Land cover status (2015) within the suitable habitats of Magnolia officinalis subsp. biloba.

The large-scale geographical distribution of vegetation is heavily influenced by the climate (Manthey and Box 2007; Punyasena et al. 2008). Amongst the myriad climatic variables, rainfall and temperature are the major factors which influence the distribution of vegetation; in particular, maximum and minimum rainfall and average temperature greatly limit the distribution and dispersal of vegetation (Thuiller et al. 2003; Crawford 2008). In our study, we found that extreme temperatures and precipitation play important roles in determining the distributions of M. o. subsp. officinalis and M. o. subsp. biloba (Table 1). However, the most influential variables varied between the two subspecies. For M. o. subsp. officinalis, temperature (min. temperature of coldest month, mean temperature of driest quarter, mean temperature of coldest quarter) was the most significant determinant of distribution. Conversely, in M. o. subsp. biloba, precipitation (annual precipitation, precipitation of wettest month, precipitation of driest month, precipitation of wettest quarter, precipitation of driest quarter, precipitation of coldest quarter) had a greater influence on distribution than temperature.

The reasons for this difference may be due to the different habitat preferences of the two subspecies. Habitats suitable for M. o. subsp. officinalis are mostly located in northern and central subtropical regions at high altitudes (Fig. 2). These areas are vulnerable to the Siberian cold current, resulting in extremely low temperatures in winter (Table 2). Therefore, extremely low temperature is the most important factor which limits the survival and spread of M. o. subsp. officinalis. Conversely, M. o. subsp. biloba prefers warmer, more southern regions (central subtropical areas) at low altitudes (Table 2). In southern areas, where low temperatures become less limiting, precipitation becomes the most important factor which prevents M. o. subsp. biloba from spreading to more northerly regions (Fig. 2). We also found that temperature fluctuations greatly influence the distributions of both M. o. subsp. officinalis and M. o. subsp. biloba (Table 1). This finding was consistent with the hypothesis that temperature stability, in combination with sufficient rainfall, can promote colonisation by plant species (Slik et al. 2009; Zhang et al. 2013).

The output ASCII grid produced by Maxent is continuous probability data ranging from 0 to 1. A threshold is needed to transform the probability data into binary data (0/1) and to acquire the information about the spatial distribution of suitable habitats. In previous studies, 0.5 (Waltari and Guralnick 2009) and 0.8 (Ramírez-Barahona et al. 2009) were often used as arbitrary thresholds to transform the probability into binary data (0/1). Other thresholds, such as the minimum predicted value (Phillips et al. 2006), the 10^th percentile training presence threshold (Brito et al. 2009), the 20^th percentile training presence threshold (Donegan and Avendaño 2010), thresholds which result in a sensitivity of 95% (Newbold et al. 2009) and the maximisation of the sum of sensitivity and specificity (maxSSS) (Liu et al. 2013; Liu et al. 2016), have been utilised.

Maxent uses a variety of methods to determine thresholds, including the minimum training presence, 10^th percentile training presence, equal training sensitivity and specificity, maximum training sensitivity plus specificity and so on (Suppl. material 1: Table S2). In this study, all thresholds estimated by Maxent were lower than 0.5 (Suppl. material 1: Table S2 and Table S3). Li et al. (2014) deemed that grids assigned probabilities higher than 0.6 could be reasonably defined as suitable habitats. Additionally, Liu et al. (2013, 2016) investigated which thresholds could be used confidently to transform probability data into binary data (0/1), reporting that the thresholds calculated by most algorithms are lower than 0.8.

Wilson et al. (2005) contend that the appropriate threshold should be selected according to the purpose for developing the species distribution model. Our primary purpose in this research was to find the most suitable habitats for M. o. subsp. officinalis and M. o. subsp. biloba to aid in the conservation and plantation of both subspecies. Therefore, we used a larger threshold (0.7) to transform the continuous probability data into binary in order to ensure that only the most suitable habitats were selected, based on the assumption that areas exceeding the larger thresholds would be better able to support the conservation and plantation.

Nature reserves provide effective refugia for wild plants. In this study, only a small number of habitats, suitable for either M. officinalis subsp. officinalis or M. officinalis subsp. biloba, were located in national nature reserves (Fig. 2). However, due to a lack of data, provincial, municipal and county-level nature reserves, forest parks, scenic spots and geological parks were unable to be included in our conservation area analysis. Therefore, the percentage of suitable habitats which are protected may in fact be higher than 8.4% (M. o. subsp. biloba) and 3.4% (M. o. subsp. officinalis).

Overlay analysis of land cover data and the locations of suitable habitats found that the majority of land in suitable habitats was woodland and forest (Figs 3, 4). However, plantations, including mulberry gardens, orchards and tea gardens were categorised within the “woodland and forest” category in the used datasets. It is thus difficult to know the exact proportion of plantation forests in the woodland and forest category. A considerable amount of land in suitable habitats was also used for agricultural land, residential land and mining and industry land (Figs 3, 4). Similarly, a large amount of land in suitable habitats was grassland (Figs 3, 4). In southern China, most grasslands are the result of secondary vegetation growth following deforestation. Therefore, we can conclude that human activities likely cause significant disturbances in habitats suitable for M. o. subsp. officinalis and M. o. subsp. biloba.

The geographical distributions of plants are heavily influenced by numerous environmental factors including climate, hydrology, soil, human activity and other factors (Gaston 2003). In our study, only climatic and terrain variables were used to assess habitat suitability. This limitation may have introduced some uncertainties into our assessments of habitat suitability. In addition, the latitudes and longitudes of some specimen collection sites were not directly recorded by GPS, being instead inferred from the names of the villages and towns near the collection sites which were recorded as background information for each specimen. Towns and villages are usually distributed in flat areas with low altitudes, which may have resulted in a tendency to assign lower altitudes than the actual altitudes of the sites from which the specimens were collected.

The habitat suitability assessment, conducted in this study, provides a scientific basis for the selection of priority protected areas in the conservation of M. o. subsp. biloba and M. o. subsp. officinalis. In addition to habitat suitability, the genetic diversity of endangered species needs to be considered when selecting priority areas for conservation (Cires et al. 2013). Genetic diversity plays a pivotal role in the maintenance of species populations (Beardmore 1983, Antonovics 1984, Yu et al. 2011). As we develop strategies for the conservation of threatened species, it is necessary to gain further insight into the genetic diversity in individual populations and that existing between populations (Hogbin and Peakall 1999; Yu et al. 2011). Previous studies have found that genetic diversity and population size are not always significantly correlated and the genetic diversity of populations of different sizes should be fully evaluated (Zhang et al. 2010). In conservation, proper attention should be paid to populations occupying smaller areas which have high genetic diversity, as they could be used to effectively reduce the loss of genetic resources (Zhang et al. 2010; Yu et al. 2011).

Magnolia officinalis subsp. officinalis and M. officinalis subsp. biloba are closely related subspecies. However, there are marked differences in the geographical distributions of these two subspecies (Liu 1996). Few studies have investigated the characteristics and distributions of habitats suitable for these two subspecies. We mapped suitable habitats for both subspecies in China in this research (Fig. 2). Based on this information, we conclude that these identified suitable habitats should be selected as conservation priority areas in China to conserve Magnolia species or more natural reserves should be established to conserve these suitable habitats. This is very important for the in situ conservation of these subspecies. Moreover, the findings of our study can assist in the selection of suitable areas in which to cultivate these subspecies.

China is the ancestral home of the Magnoliaceae, with more than 40% of Magnoliaceae species having originated in southwest China (Wang and Jiang 2001; Cicuzza et al. 2007). In the past few decades, the survival of numerous Magnoliaceae species has been severely reduced by the activity of humans and many are greatly threatened by over-utilisation (for timber harvesting, bark, flower bud collection etc.) and habitat fragmentation (Wang and Jiang 2001; Cicuzza et al. 2007; Cires et al. 2013). According to the “China Biodiversity Red List,” 71.7% of Magnoliaceae species are defined as “Threatened” (Ministry of Ecology and Environment of China and Chinese Academy of Sciences 2013). Therefore, it is necessary to properly evaluate habitat suitability and genetic diversity for species in this family and to provide scientific support for their conservation.

In this study, we found that the environmental variables, which influence species distributions, are different for each subspecies. The distribution of M. officinalis subsp. officinalis was primarily determined by variation in minimum temperatures, while the distribution of M. officinalis subsp. biloba was primarily determined by variation in precipitation.

We identified the habitats suitable for both subspecies and found that the two subspecies have distinct habitat preferences. Compared to M. o. subsp. biloba, M. o. subsp. officinalis is found in more northerly areas, grows at higher altitudes and is able to survive in areas experiencing greater fluctuations in ambient temperature, lower extreme temperatures, less precipitation and greater fluctuations in precipitation.

The results of this analysis could provide useful information to support the in situ conservation of both M. o. subsp. officinalis and M. o. subsp. biloba and could aid in the selection of cultivation sites.

In the future, field survey data on the distribution of Magnolia species should be included in the assessment of habitat suitability to offset the deficiencies with regard to the specimen data. Genetic diversity assessment should be performed, together with habitat suitability assessment to provide stronger scientific support for the conservation of Magnolia species.