Predicting Impacts of Climate Change on Suitable Distribution of Critically Endangered Tree Species Yulania zenii in China

1. Introduction

The spatial distribution of plants is affected by climate, and accordingly their geographical range will alter with climate change over time (Dyderski and Pawlik, 2020). The historical fluctuations in global climate have significantly influenced the current distribution patterns of most species (Veloz et al, 2012). As early as the last interglacial period, approximately 12,000-14,000 years ago, when global temperature was around 2°C higher than pre-Industrial Revolution levels, numerous plant species survived climate change by migrating to previously unoccupied habitats (Turney and Jones, 2010). The Middle Holocene is the latest typical great warm period, in which its average annual temperature exceeds that of the present day. During the late Holocene Megathermal, most plant populations migrated from low to high latitude (Feng et al., 2022). In the last few decades global warming is becoming more pronounced, thus exacerbating such impact on plants’ distribution. According to the recently released Sixth Assessment Report of the United Nations Intergovernmental Panel on Climate Change (IPCC), the global temperature is expected to rise more than 1.5°C by 2030 (Kanitkar et al., 2023). Therefore, the global warming may alter the environmental conditions and population sizes for plants, thereby triggering a shift in distribution range (Ren et al., 2020).

Climate change plays a pivotal role in influencing plants’ distribution pattern on a regional scale (Qin et al., 2017). Generally, endangered tree species are vulnerable to climate change compared with those widely distributed. This is primarily attributed to the restricted habitats and small population size of such threatened trees, making them poor in adaptation in the face of climate change (Vincent et al., 2020). Furthermore, these species usually have discontinuous or highly fragmented habitats, resulting in their difficulties in migration, and even rendering them at risk of extinction in the context of climate change. For instance, Zhao et al. (2020) employed the MaxEnt model to predict the suitable distribution of the endangered Carpinus tientaiensis, and pointed out that most of its core distribution will drastically shrink in the next several decades (Zhao et al., 2020). Yan and Zhang (2022) stated that, for the endemic and endangered Parrotia subaequalis, its suitable habitats in China may become more fragmented under future climate conditions (Yan and Zhang, 2022). As a result, endemic tree species may be more vulnerable to climate change.

Yulania zenii (W. C. Cheng) D. L. Fu, endemic to China, is a deciduous tree in the Magnoliaceae family. This tree is hysteranthous, and usually blooms from March to April every year. It has nine nearly spoon-shaped pink or white tepals for each flower, with aromatic odor in spring, and it has cylindric aggregate follicles purplish red turning brown in autumn (Figure 1). Its twig is green when young, and becomes purple when old. Its adult tree has a tall and straight trunk, up to nearly 30 m in height (Wu et al., 2008). Therefore, this Chinese magnolia has high value in ornamental horticulture. Meanwhile, it can serve as a valuable material for elucidating molecular phylogeny and floral evolution in Magnoliaceae (Chu, 2021). Over the past several decades its wild populations have been increasingly threatened by reproductive barrier and anthropogenic interference, especially by tourism development in subtropical eastern China (Zhang, 2007). Hence, it has been listed as the second-grade species in the List of National Key Protected Wild Plants since 2021 (https://www.forestry.gov.cn/). In 2023 it was listed as one of the most critically endangered China's 100 species, i.e. Plant Species with Extremely Small Populations in China (Zhang and Yi, 2023). Furthermore, it has been classified as a "critically endangered" (CR) species in the IUCN Red List (https://www.iucnredlist.org/).

Generally, Yulania zenii is considered to be very limited in geographical distribution. According to Flora of China (Vol.7), it only occurs in Baohua Mountain of Jurong City, Jiangsu Province, eastern China (Wu et al., 2008). Hereafter, most researchers hold the same opinion over its distribution, such as Yin, 2013, Liu, 2015, Zhang et al., 2022. Even in 2023, Jin et al. stated that this species was just distributed in Mt. Baohua in their monograph "National Key Protected Wild Plants of China" (Jin et al., 2023). However, recent studies suggest that the wild populations of Yulania zenii are also present in other provinces. For example, this species was found in North Luoxiao National Forest Park in 2022 (Peng, 2022) and Shending Mountain Provincial Forest Park in 2023 (https://www.forestry.gov.cn/). Both locations are situated in Yueyang City, Hunan Province. The other example is from mountainous area of northwestern Ruichang City, Jiangxi Province (https://www.jiujiang.gov.cn/). There were 11 individuals of Yulania zenii reported therein. It seems that this species currently appears in at least three Chinese provinces including Jiangsu, Jiangxi in eastern China and Hunan in central China. Therefore, we think that its actual distribution range should be much wider than its recorded distribution in China.

Species distribution models (SDMs) serve as a crucial instrument in studying the potential impact of climate change on species distributions. For a certain species, SDMs can predict its suitable area across space and time by integrating its occurrence records with corresponding environmental information (Yackulic et al., 2013). SDMs can be divided into two categories: mechanistic and correlative models. Mechanistic models typically rely on data pertaining to the predicted species' life history, functional traits, or physiological responses to environmental fluctuations. However, it is challenging to obtain such data, and furthermore they may not always be representative on most occasions. In contrast, correlative models primarily utilize the known presence (or non-presence) point data of species. Nowadays, correlative models are widely applied due to the accessibility of these data of species distribution and their environment in relative to mechanistic models (Shabani et al., 2016). Currently, correlative models which are widely used include the Maximum entropy model (MaxEnt), Generalized Additive model (GAM), Generalized linear model (GLM), Random Forest model (RF) and so on (Sillero N et al., 2021). Among them, MaxEnt stands out as one of the most widely utilized SDMs because of accurate prediction, simple operation, small sample sizes, and effective noise reduction (Kong et al., 2019; Rathore and Sharma, 2023). Particularly, this model has now been extensively used for predicting the potential distribution of rare and endangered plants, their key environmental factors, as well as suitable habitat selection for introduction and cultivation under future climate scenarios (Hills et al., 2019; Lu et al., 2022). Even so, each model has its own advantages and disadvantages based on different principles and algorithms. More recently, ensemble models have been developed to predict the potential distribution of endangered species (Kiser et al., 2022).

In this study, we first collected data on the distribution points of Yulania zenii and related environmental factors. Then, we used Biomod2 for pre-modeling. Subsequently, we employed the selected model to predict its potential suitable area in China. The objective of this study is to determine its current potential distribution, its response to various climate scenarios in the past and future. More specifically, this study is aimed to: (1) To identify the key environmental factors affecting the spatial distribution of Yulania zenii; (2) To forecast its potential ranges under past, current and future climate scenarios, and determine its centroid migration routes; (3) To provide corresponding countermeasures for the conservation and management of Yulania zenii.

2. Materials and Methods

2.1. Species Occurrence Data

The data regarding the wild distribution of Yulania zenii were primarily obtained through the following approaches: (1) Investigating in field. In the past three years, we carried out comprehensive surveys for Yulania zenii wild populations in Anhui, Jiangsu, Jiangxi, Zhejiang, and other provinces of eastern China to acquire their spatial localities; (2) Searching through resource sharing platforms. These include the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/), the National Specimen Information Infrastructure (NSII, http://nsii.org.cn/), the Chinese Virtual Herbarium (CVH, https://www.cvh.ac.cn/), and the Plant Photo Bank of China (PPBC, http://ppbc.iplant.cn/); (3) Consulting published literature and relevant reports (Peng, 2022). We used the key words of Yulania zenii’s specific name, Latin name and its synonym (i.e. Magnolia zenii) in Flora of China, provincial floras, and related checklists to gather its distribution points. In this way, we initially collected 32 natural distribution records of Yulania zenii.

After eliminating erroneous and duplicate distribution records, we employed the tool in the SDMs toolbox (version 2.6), which is called Spatially Rarefy Occurrence Data for SDMs, to ensure that each 1 km × 1 km grid contained only one distribution point (Radosavljevic and Anderson, 2014). This approach to rarefying distribution point is supposed to align with the resolution of the environmental data, aiming to minimize the number of distribution points existing spatial autocorrelation and prevent overfitting of the model (Brown et al., 2017). In the end, we acquired the latitude and longitude data of 12 distribution points of Yulania zenii (Figure 2; Table S1). To facilitate data preparation for modeling, the distribution data file was converted into the CSV format.

2.2. Environmental Variables

38 selected environmental factors were classified into three categories: climate, terrain, and soil. For climate data, we opted for bioclimatic factors that had great biological significance (Poirazidiset al., 2019). The past period (the Last Interglacial period, approximately 12,000-14,000 years ago, and the Middle Holocene, around 6000 years ago), the current and future periods (specifically, the 2050s and 2070s) were all considered in this study. 19 bioclimatic factors of these periods were downloaded from the Worldclim (https://www.worldclim.org/) (version 1.4). Then, we unified the resolution to 30s (1 km × 1 km), primarily to guarantee the precision during modeling. Given that a single climate model is unrepresentative in predicting future climate scenarios, we chose to use a combined model that integrates multiple climate models (Chen et al., 2020). Accordingly, the bioclimatic data for the future periods were obtained by calculating the equally-weighted average values of three global climate models: the Beijing Climate Center Climate System Model version 1.1 (BCC-CSM1-1), the Community Climate System Model version 4 (CCSM4), and an Earth system model based on the Model for Interdisciplinary Research on Climate (MIROC-ESM). In the context of future climate, there are four typical Representative Concentration Pathways (RCPs) for greenhouse gases: RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5. These pathways represent different scenarios of climate change, ranging from the lowest to the highest emission scenario. Furthermore, they have been widely employed in examining the species' response to climate change (Zhang et al., 2021). Both RCP4.5 and RCP6.0 represent intermediate and stable climate change scenarios, with RCP4.5 being considered more significant than RCP6.0 (Moss et al., 2010). Therefore, we reduce our selection to three commonly utilized RCPs: RCP 2.6 (representing a moderate emission scenario), RCP 4.5 (signifying a medium and stable emission scenario), and RCP 8.5 (indicating a high emission scenario). Ultimately, climate data were collected for six emission scenarios encompassing two future periods: 2050s RCP 2.6, 2050s RCP 4.5, 2050s RCP 8.5, 2070s RCP 2.6, 2070s RCP 4.5, and 2070s RCP 8.5.

Topographic data comprised elevation and slope. In general, topographic factors change so slowly over time that they can be negligible (Stanton et al., 2012). Therefore, they were incorporated into the model as static variables to enhance accuracy of projected results. The elevation data were downloaded from the WorldClim (https://www.worldclim.org/), with spatial resolution of 30s. The slope data were extracted from DEM (Digital Elevation Model) Data downloaded from the National Earth System Science Data Center (http://www.geodata.cn). We also downloaded the Chinese soil dataset (version 1.2) from the National Qinghai-Tibet Plateau Scientific Data Center (http://www.tpdc.ac.cn/zh-hans/). Finally, 17 types of topsoil data (0-30 cm) were selected from the dataset as the soil data used in this study.

Considering that the Last Interglacial Period and the Middle Holocene represent two significant paleoclimatic epochs, during which the Earth's environment underwent profound transformations, only 19 bioclimatic factors were chosen in these two epochs for subsequent modeling (Yang et al., 2022). As for the current and future periods, we still selected three distinct types of environmental data: climate, topography, and soil.

The environmental data from the aforementioned three types were uniformly standardized under the WGS1984 coordinate system. The function of Extract by mask and Clip in the ArcGIS 10.8 software was applied to make sure that the data were only confined to the territory of China. Subsequently, the resampling tool was employed to harmonize the resolution of all data to a level of 30s. Concurrently, Pearson correlation analysis was used to mitigate collinear interference among correlated environmental factors, ensuring that redundant information did not contaminate the model prediction process (Sillero and Barbosa, 2021). This approach significantly enhanced the precision of prediction outcomes. The specific operation steps were outlined as follows: Initially, the distribution data of Yulania zenii and environmental data were input into the model for a preliminary simulation, which allowed us to obtain the original contribution rate of each environmental factor. Afterwards we utilized the Spatial Analyst tool in ArcGIS 10.8 to extract the value of each environmental factor precisely at all distribution points. Lastly, we conducted a test in R 4.3.1 to determine the Pearson correlation coefficient (r) between these environmental factors. Environmental factors with low contribution rate among these factors with correlation coefficient | r | ≥0.8 were subsequently eliminated from the analysis (Kiser et al., 2022). Table 1 presents the environmental factors for subsequent modeling in various periods and their corresponding contribution rates.

2.3. Modeling Process

Biomod2 is a multi-model ensemble platform that relies on ten widely used modeling techniques (Zhao et al., 2021). We employed Biomod2 to model the potential distribution of Yulania zenii in the current period. By utilizing 75% of the distribution points as the training set, we conducted 10 repeated operations to obtain the AUC and TSS values for each model. By doing so, this practice can enable us to accurately assess its performance. We then noticed that MaxEnt performed much better than each of the others (Table 2).

Prior to formal modeling, it is imperative to select suitable values for the Regularization Multiplier (RM) and Feature Class (FC) within the MaxEnt model (Lu et al., 2022). The regularization multiplier ranged from 0.5 to 4.0, with increments of 0.5, resulting in eight regularization multipliers. The MaxEnt model encompassed five distinct feature types: Linear (L), Quadratic (Q), Hinge (H), Product (P), and Threshold (T) (Zhou et al., 2023). After optimizing the MaxEnt model, the final parameter combination was determined to be RM=1 and FC=LQ. Subsequently, the modeling operation was executed in MaxEnt 3.4.4. 25% of the distribution points were randomly chosen as the test set while the remaining 75% served as the training set. To guarantee the precision of prediction results, we employed the Bootstrap method to replicate the calculation ten times. Then we selected "Cloglog" as the output mode and "ASC" as the preserved format.

To evaluate model performance, we utilized the area under the receiver operating characteristic curve (AUC) and the true skill statistic (TSS) (Allouche et al., 2006). The AUC value typically ranged from 0 to 1, and a value closer to 1 indicated higher accuracy. It can be categorized as follows: failing (0.5-0.6), poor (0.6-0.7), fair (0.7-0.8), good (0.8-0.9), and excellent (0.9-1.0) (Singh et al., 2021). TSS values varied between -1 and +1, with values closer to 1 indicating superior performance, while values closer to or below 0 suggested inferior performance. Models can also be divided into the five groups in terms of TSS: excellent (TSS > 0.8), good (0.6-0.8), fair (0.4-0.6), poor (0.2-0.4), and failing (TSS < 0.2) (Wang et al., 2023). The overall performance of the model was assessed by calculating the average AUC and TSS values obtained from 10 replicates.

2.4. Geospatial Data Analysis

The results of 10 average operations generated by MaxEnt were imported into ArcGIS 10.8 for visualization. Given the endangered status of Yulania zenii, we employed the approach of maximizing the sum of sensitivity and specificity (max SSS) to set the threshold for suitable area (Liu et al., 2013; Xu et al., 2021). This threshold selection was considered to be highly effective when using model with presence-only data to divide suitable area into different levels (Liu et al., 2013). In light of the threshold of the max SSS (0.2639), we categorized the potential distribution of Yulania zenii into four levels: unsuitable (0.0-0.27), low suitable (0.27-0.51), moderately suitable (0.51-0.76), and highly suitable (0.76-1.00) (Lu et al., 2022). Next, we calculated the suitable area of each type.

Centroid migration can characterize distribution changes of species under different climatic scenarios. We employed the SDMtoolbox in ArcGIS 10.8 to simulate the changing situation of this species' centroid migration under nine climate scenarios (Brown et al., 2017). Moreover, we also determined their direction and distances in different periods.

3. Results

3.1. Model Performance

We calculated the AUC and TSS values of each model except GAM, which could not run due to its few occurrence records (Luo et al., 2017). Among the nine models (Table 2), MaxEnt had the minimum standard deviation, exhibiting its excellent stability. More importantly, only MaxEnt simultaneously satisfied the criteria of AUC > 0.9 and TSS > 0.8, thus indicating its superior performance compared to the other models within Biomod2. Consequently, MaxEnt was ultimately chosen for subsequent modeling.

Using the MaxEnt model, we calculated its AUC and TSS values following 10 simulations under various climate scenarios (Table 3). On the whole, the average AUC value was 0.9820, and the average TSS value was 0.9085, suggesting that the MaxEnt performed well.

3.2. Main Environmental Factors

We determined the percentage contribution of selected environmental factors using MaxEnt (Table 1). Among the key environmental factors affecting the distribution of Yulania zenii at present, the top three were Bio2 (Mean diurnal range), Bio15 (Precipitation seasonality of variation coefficient), and Elevation. Their contribution rates were respectively 32.9%, 21.1%, and 14.8%, with the cumulative proportion 68.8%. In the Last Interglacial Period, the key environmental factors were Bio15 (43.1%), Bio3 (Isothermality, 27.1%), and Bio2 (25.2%), amounting to 95.4%. During the Middle Holocene, the key environmental factors were Bio15 (39.9%), Bio2 (32.1%), and Bio3 (19.4%), respectively; their total contribution rate was up to 91.4%. Under different climate scenarios, Bio2 and Bio15 were recognized as common leading factors.

When the existence probability was greater than 0.51, it corresponded to the moderately and highly suitable area, which was conducive to the growth of Yulania zenii. We opted to utilize the response curve generated by MaxEnt to reveal the relationships between environmental factors and existence probability in the current period (Figure 3). The existence probability of this species showed a trend of first increasing and then decreasing with the change of mean diurnal range (Bio2), and reached a maximum at 6.88 °C. As shown in Figure 3a, it was optimal for the survival of Yulania zenii when the range was 4.87-9.06 °C. Its existence probability gradually decreased with increasing of the precipitation seasonality of variation coefficient (Bio15), and its adaptation range was 20-64 (Figure 3b). The suitable condition for elevation factor ranged from 0 to 500 m (Figure 3c). And the existence probability dropped sharply as the elevation rose.

3.3. Current Potential Suitable Distribution

The current suitable area for Yulania zenii was primarily concentrated in southern Jiangsu, southern Anhui, northern Zhejiang, eastern Hunan, central Hubei, and northern Jiangxi (Figure 4); additionally, few suitable area was forecasted in Sichuan and Chongqing in southwestern China. Overall, the suitable distribution was predominantly concentrated in eastern and central China. Moderately suitable area was relatively continuous while highly suitable area was noticeably fragmented (Figure 4). For Yulania zenii, its total suitable area was 14.68 × 10⁴ km², taking only 1.53% of China's total territory. Furthermore, its highly suitable area only accounted for 0.55% (Table 4).

3.4. Potential Suitable Distribution in the Past

During the Last Interglacial, the suitable habitats for Yulania zenii were primarily concentrated in eastern Hunan, southern Anhui, southern Jiangsu, northern Zhejiang, and the junction area of Jiangxi, Anhui, and Hubei in China (Figure 5). On the whole, its suitable distribution exhibited a continuous pattern. However, the highly suitable area displayed fragmented distribution in certain parts of Zhejiang, Anhui, and Hunan provinces. The total suitable area amounted to 23.58 × 10⁴ km², which increased by 60.63% compared with the current. The highly suitable area occupied 11.42 × 10⁴ km², with a sharp increase of 117.52% in relative to the current (Table 4).

In the Middle Holocene, the suitable habitats were mainly distributed in southern Anhui, southern Jiangsu, northern Jiangxi, eastern Hubei, central and eastern Hunan, the coastal areas of Zhejiang, the junction of Zhejiang and Anhui provinces, the junction of Hunan and Hubei provinces, and the junction of three provinces including Jiangxi, Hubei, and Anhui (Figure 5). Compared with the Last Interglacial, the suitable distribution of the Middle Holocene exhibited a more continuous pattern, with less fragmentation. The total suitable area amounted to 35.68 × 10⁴ km², showing an increase of 143.05% compared to present. Moreover, the highly suitable area was 14.73 × 10⁴ km², with a dramatic increase of 180.57% compared to the current (Table 4).

In brief, from the Last Interglacial to the Middle Holocene and subsequently to the current, the suitable area for Yulania zenii initially increased and then decreased, and the same went for its highly suitable area. The suitable area decreased significantly in the current compared to the two past periods. Furthermore, the current habitat became increasingly fragmented due to climate change.

3.5. Potential Suitable Distribution in the Future

The suitable area were primarily concentrated in southern Jiangsu, southern Anhui, northern Zhejiang, eastern Hunan, central Hubei, and northern Jiangxi in the future, which largely aligned with the current suitable distribution. However, suitable area reduced in these provinces with varying degrees. In addition, we noticed that the highly suitable area expanded in the coastal area of Zhejiang Province under three scenarios, namely 2050s RCP4.5, 2050s RCP8.5, and 2070s RCP8.5 (Figure 6).

Under various emission scenarios, the average suitable area amounted to 11.71 × 10⁴ km², indicating a decrease of 20.26% compared with the current. Among six future scenarios, the suitable area reduced least in the 2050s RCP4.5, with a decline of only 2.04% relative to the current. Conversely, the suitable area underwent the largest loss in the 2070s RCP4.5, decreasing by nearly 40% by comparison with the current. Furthermore, the highly suitable area reduced by varying degrees. The mean highly suitable area under six future scenarios amounted to 3.88 × 10⁴ km², which decreased by 26.07% compared with the current condition. The reduction was the least (10.67%) in the 2050s RCP4.5, while it was the highest (47.81%) in the 2070s RCP4.5 (Table 4).

In the future periods, there was less reduction of the suitable area in the 2050s (16.83%) than in the 2070s (23.68%) (Table 4). Additionally, the rate of decline in the suitable area differed under various emission scenarios in the same period. In the 2050s, the suitable area exhibited the smallest decrease under RCP4.5, whereas the largest decrease was observed under RCP2.6. However, the suitable area decreased the least under RCP8.5 and the most under RCP4.5 in the 2070s (Table 4).

Overall, the suitable habitats for Yulania zenii exhibited a decreasing trend in the future (2050s and 2070s), with severe habitat fragmentation.

3.6. Centroid Migration under Different Scenarios

The centroid coordinate of Yulania zenii in the current (31.177°N, 111.656°E) was situated in Yuan'an County, Yichang City, Hubei Province. From the Last Interglacial (30.041°N, 115.156°E) to the Middle Holocene (31.567°N, 111.271°E) and then to the current, the centroid migrated northwestward by 407.69 km at first and subsequently southeastward by 56.64 km. Under RCP2.6, the centroid migrated southeastward with a distance of 85.62 km in the 2050s (30.845°N, 112.467°E), while it shifted 142.53 km to the southeast in the 2070s (30.592°N, 112.986°E). Under RCP4.5, it migrated 51.97 km southeastward in the 2050s (30.900°N, 112.096°E), and 132.31 km to the southeast in the 2070s (31.001°N, 113.031°E). Under RCP8.5, the centroid moved northeastward with a distance of 29.65 km in the 2050s (31.417°N, 111.792°E) and 57.51 km in the 2070s (31.303°N, 112.243°E).

Collectively, the centroid exhibited a sinuous migratory pattern, initially shifting northwestward from the Last Interglacial to the Middle Holocene, and subsequently southeastward towards the current. From the current to future, the migration direction did not show the consistent trend. Specifically, the centroid shifted southeastward under RCP2.6 and RCP4.5, whereas it migrated northeastward under RCP8.5 (Figure 7).

4. Discussion

4.1. Model Evaluation

Ten models from Biomod2 package were employed to preliminarily simulate the potential geographical distribution of endangered Yulania zenii in China. Only MaxEnt model had the AUC > 0.9 and simultaneously TSS >0.8 (Table 2), exhibiting superior performance. Therefore, MaxEnt was selected for final modeling.

Then, we removed collinearity of the 38 environmental factors and optimized the model parameters (i.e. FC and RM) to ensure that the MaxEnt model had high accuracy and reliability (Table 1). The results showed that the mean values of both AUC and TSS were all greater than 0.9 under the nine climate scenarios (Table 3), indicating that the MaxEnt performed well. In addition, each of them had a very low standard deviation.

Besides, the forecasted current distribution is in line with the known occurrence records of Yulania zenii. Therefore, we used the optimized MaxEnt to predict the potential areas under past, current and future climate scenarios, respectively.

4.2. Key Environmental Factors

Bio2, Bio3 and Bio15 were the main environmental factors in two paleoclimatic periods. During the Last Interglacial period, the contribution rates of the three factors were 25.2%, 27.1% and 43.1%, respectively. In the Middle Holocene, their contribution rates were 32.1%, 19.4% and 39.9%, respectively. Under the current climate conditions, the main environmental factors were Bio2, Bio15 and Elevation, with the contribution rates of 32.9%, 21.1% and 14.8%, respectively (Table 1). This indicated that among the three types of environmental factors (climate, terrain and soil), climate may play a more significant role in limiting Yulania zenii’s distribution. Furthermore, Bio2 and Bio15 were identified as the key factors during the past and current periods. Shi et al. (2021) simulated the future suitable distribution of Magnolia wufengensis using three types of environmental factors and concluded that climate was the main influencing factor (Shi et al., 2021), which is consistent with our results.

In the current, the key environmental factors restricting the distribution of Yulania zenii included Bio2, Bio15 and Elevation. When the mean diurnal range (Bio2) is 4.87-9.06 °C, the precipitation seasonality of variation coefficient (Bio15) varies between 20 and 64, and the altitude ranges from 0 to 500 m, it is beneficial to the survival of this species (Figure 3).

In general, the diurnal temperature range is the largest in the low latitudes (average 12 °C), followed by the middle latitudes (7-9 °C), and the smallest in the high latitudes (3-4 °C) (Kong, 2020). At present, the distribution area of Yulania zenii is concentrated the three provinces of Jiangsu, Jiangxi, and Hunan (Figure 2). Their geographical range is between 24°29 ' N and 35°08' N, which is within the middle latitude region. Therefore, the mean diurnal range herein corresponded to the suitable range for this species. As shown in Figure 3b, with the increase of the variation coefficient of precipitation seasonality (Bio15), the existence probability of Yulania zenii showed an obvious decline. Yulania zenii is vulnerable to Bio15, and therefore it grows well in such area with little variation coefficient of precipitation seasonality. As far as altitude is concerned, Yulania zenii is largely distributed at low altitudes. An example is its wild population from Mt. Baohua, Jurong City, Jiangsu Province, in which it mainly occurs in hilly areas with an altitude of 220 m (Wu et al., 2008). Another example is its wild population from Ruichang City, Jiangxi Province, with the elevation of about 400 m. The third example is from Mt. Shending in Miluo City, Hunan Province, where its wild individuals have been recorded in recent years. And this tree only reaches up to 464.3 m of elevation in distribution. Therefore, the altitude of Yulania zenii populations from known sites is consistent with our predicted elevation range (0-500 m) (Figure 3c).

Therefore, we have identified the key environmental factors affecting the distribution of Yulania zenii for the first time, and further determined their corresponding optimal range. Namely, this tree species is suitable to grow in the low altitude areas with small variation in the mean diurnal range and precipitation seasonality.

4.3. Current Suitable Area of Yulania zenii

The MaxEnt model predicted, for the first time, the current suitable area of Yulania zenii was 14.68 × 10⁴ km² (Table 4), only accounting for 1.53% of China's total territory. These suitable areas were mainly distributed in southern Jiangsu, southern Anhui, northern Zhejiang, eastern Hunan, central Hubei and northern Jiangxi (Figure 4). At present, it is generally recognized that this species only occurs in Mt. Baohua in Jurong City, Jiangsu Province (Jin et al., 2023; Yu et al., 2023). This indicates that the actual distribution of Yulania zenii is larger than the known in China.

For one thing, there are many similarities among closely related species from Yulania in ecological and morphological characteristics (Wang, 2003). This makes it difficult to distinguish Yulania zenii from other similar species in field because this tree usually has short florescence in subtropical forests with high canopy density (Chu, 2021). Furthermore, such a species often has small populations which are restricted in segregated habitats within mountainous regions due to human interference. For the other thing, this species has bright red seeds with arils at maturity, which can attract frugivorous birds to disperse its seeds (Zhang, 2007). Accordingly, this may efficiently expand its range.

4.4. Suitable Area Change in the Past and Future

Yulania zenii had 23.58 × 10⁴ km² of suitable area in the Last Interglacial, and it extended to 35.68 × 10⁴ km² in the Middle Holocene. Compared to the current period, its suitable area in these two historical periods increased by 60.63% and 143.05%, respectively (Table 4). Hence, its suitable area of the past is significantly larger than that of the current. During the Last Interglacial period, the climate was warm and arid with low precipitation (Yan et al., 2022). This may constrain the growth of this species. The Middle Holocene is the last great warm period, and its climate is warm and humid, with high precipitation (He et al., 2022). During this period it may be more conducive for Yulania zenii to grow and reproduce, thus its suitable area increased significantly.

On the contrary, Yulania zenii had 11.71 × 10⁴ km² of the mean suitable area in the six future climate scenarios (2050s and 2070s), and it decreased by 20.26% on average compared to the current distribution (Table 4). The average suitable area was 12.21×10⁴ km² in the 2050s while it was 11.20 × 10⁴ km² in the 2070s. Therefore, Yulania zenii will significantly decline in suitable area in the future scenarios. In the coming future, global temperature is constantly rising, coupled with frequent extreme weather events and intensive human activities (Wang et al., 2023), thereby the suitable distribution of Yulania zenii is expected to continue shrinking.

By and large, the centroid of Yulania zenii will migrate to the southeast under the future climate scenarios (Figure 7), which is consistent with the whole migrating direction of Allium mongolicum in China in the future (Lang et al., 2023). This may be attributed to such a tree trait that Yulania zenii prefers to grow in a warm and humid habitat (Wang, 2020).

Overall, Yulania zenii has been shrinking in suitable area from the past to the current till the future. Our result is in line with those of other endangered tree species such as Pseudotaxus chienii (Zhang et al., 2020), Semiliquidambar cathayensis (Ye et al., 2020), and Lonicera oblata (Wu et al., 2021). Based on our MaxEnt modeling, the whole population of Yulania zenii was relatively widely distributed in the past (Figure 5), then it dramatically contracted currently (Figure 4), and further became more fragmented in the future (Figure 6). Therefore, we believe that climate change have taken a negative effect on Yulania zenii in distribution, especially in terms of suitable area and habitat integrity.

4.5. Conservation Implications for Yulania zenii

Currently, most studies have long held that Yulania zenii only occurs in Mt. Baohua, Jiangsu Province, eastern China (Yin, 2013). However, our findings demonstrate that its current potential suitable area is much larger than the known range. More recently, this species is found to present in the northern part of Mt. Luoxiao, which is located in Yueyang, Hunan Province, central China (Peng, 2022). Therefore, we propose to carry out a comprehensive survey on its wild populations in central and eastern China, particularly in southern Anhui, central Hubei, eastern Hunan, southern Jiangsu, northern Jiangxi, and northern Zhejiang. Furthermore, it is reported that Yulania zenii has higher genetic diversity at species level, which is detected by inter-simple sequence repeat (ISSR) markers, than other endangered tree species with limited geographical distribution (Chen and Nan, 2016). However, the ISSR analysis sampled just one population from Mt. Baohua, and accordingly such a practice of sampling seems unlikely to reflect the actual situation of this species. Hence, given its wide distribution range across various provinces in China, it is necessary to conduct extensive sampling to reveal its genetic diversity and structure in the coming future.

Secondly, we find that climate change has a negative impact on the distribution of Yulania zenii. As early as 2012, it was regarded as one of the most-at-risk 120 species in China which were described as Plant Species with Extremely Small Population (PSESP) (Yang et al., 2020). In 2023, it is listed on the updated national PSESP checklist with 100 species (Zhang and Yi, 2023). Although its potential distribution is larger than known, our modeling results indicate that its main distribution areas are separated from each other in different provinces of China (Figure 4). To make things worse, its suitable habitat may become more fragmented under future climate scenarios. Therefore, climate change should be taken into account to expand its population size when making the conservation plan of Yulania zenii in the future.

We also find that the key environmental factors affecting its distribution are Bio2, Bio15 and Elevation. According to the modeling results, for Yulania zenii its suitable area may be the low altitude regions with small mean diurnal range and little variation of seasonal precipitation. For Yulania zenii, its adult trees usually blossom in early spring, and have short florescence, lasting at most 15 d. Soon afterwards some carpels may be abortive during embryonic development, resulting in distorted aggregate follicles (Figure 1c). This is probably because its fruit setting rate is susceptible to precipitation. The fruit setting rate is low while the climate is dry; in contrast, it is high while the precipitation is abundant (Wang, 2020). Therefore, the leading climate factors predicted by MaxEnt (e.g. Bio2, Bio15) are in agreement with the fact that this species is mainly distributed in the warm and humid subtropical mountains in China. Accordingly, this provides valuable information on reasonable management and planting practices of Yulania zenii, particularly in climate.

There are 25 species of Yulania in the Magnoliaceae family worldwide, which occur in temperate and subtropical Asia as well as North America (Chen et al., 2020). There are 18 species of Yulania in China, which exist in the Yangtze River basin and its southern regions (Chen et al., 2020). Among them, nine species are currently under threat in light of IUCN Red List Category and Criteria. Specifically, Yulania zenii and Yulania sinostellata are classified as critically endangered (CR). Yulania dawsoniana, Yulania pilocarpa and Yulania viridula are as endangered (EN). Yulania amoena, Yulania liliiflora, Yulania sargentiana and Yulania campbellii are as vulnerable (VU) (Qin, 2020). Except Yulania campbellii, the other eight Yulanias are all endemic to China. In this study, we select Yulania zenii as the representative one. Our findings indicate that climate change may take an adverse effect on its suitability of habitats from past to current till future. Our MaxEnt modeling, for the first time, reveals the response of Yulania zenii to climate change during different periods. Therefore, this study can provide reference for other endangered Yulania species in China under the condition of climate change.

5. Conclusions

In this study, we first applied Biomod2 for preliminary modeling, and chose MaxEnt for formal modeling. We used its distribution data and three types of environmental variables (climate, terrain and soil) to project the potential distribution of the endangered Yulania zenii in China. Our results showed that the main environmental factors affecting its distribution were Bio2 (Mean diurnal range), Bio15 (Precipitation seasonality of variation coefficient) and Elevation, and that climatic factors may play a vital role in limiting its geographical distribution. For the first time we determined its current suitable area of 14.68 × 10⁴ km², accounting for only 1.53% of China's total territory. This distribution area is larger than known, which is mainly distributed in eastern and central China. During two historical periods (the Last Interglacial and Middle Holocene), Yulania zenii had much larger area than at present. Under the future scenarios, its suitable area will averagely decrease by 20.26%. Our findings confirm that from the past to current till future climate change has a negative influence on this species, especially in terms of suitable area and habitat integrity. Our study is conducive to the conservation and management of Yulania zenii, and can also provide valuable information for other endangered Yulania species in China.