Distributional Response of the Rare and Critically Endangered Ilex nanchuanensis to Climate Change in East Asia

1. Introduction

Global climate change alters the temperature and precipitation patterns on the Earth’s surface [1,2,3,4,5,6], which in turn affects plant growth and development patterns, morphological structures, and geographic distribution patterns[7,8,9,10,11,12] , ultimately endangering the diversity and stability of global ecosystems[13,14,15] . Global climate change, dominated by climate warming[16,17] , is one of the greatest threats to, and challenges for, ecosystem balance[18,19,20,21]. The Sixth Assessment Report on Climate Change published by the Intergovernmental Panel on Climate Change (IPCC) in 2021 states that the global surface temperature increased by 0.99°C from 2001 to 2020 compared to the previous century [22,23], and global warming is expected to continue to 1.5°C by 2030—2052[24]. Thomas et al. (2004) found that approximately 37% of species globally will be at risk of extinction and habitat loss by 2050 as a result of climate warming[25], while Urban (2015) found that the risk of species extinction accelerates with increasing temperatures[26], especially for some rare and endangered species with small population sizes, high habitat specificity, and severe habitat fragmentation[27,28]. Daskalova et al. (2020) found that rare and endangered species with small population sizes and limited geographic distributions will be the first to go extinct[29]. In addition, mountain ranges are crucial for global sustainable development because of their specific geographic features and diverse climates[30,31]. Global climate change has caused serious harm to species living in mountainous regions[32,33,34,35,36], and such species are more sensitive to climate change and at greater risk of species extinction[37,38,39,40]. This is especially true for the native plants of low-altitude mountain areas, which will be limited in terms of their vertical dispersal, leading to a greater risk of extinction[41]. Therefore, analyzing the response of rare and critically endangered species in low-altitude mountain areas to climate change can not only clarify the dominant climatic factors affecting the geographical distributions of species, but also be crucial for vegetation restoration and ecological conservation efforts.

Species distribution models (SDMs) are widely used in the fields of endangered species conservation and species invasion[42,43,44,45], and can efficiently and accurately analyze the effects of global changes on species distributions and diversity patterns[46]. SDMs are effective tools for predicting the effects of climate change on the distributions of endangered species in the context of ecological conservation and management[47]. Several widely used SDMs emerged in the late 20th century, among which the MaxEnt (Maximum Entropy) model not only has stable predictive performance, but also has higher accuracy[48,49]. By studying the distribution data of 226 species, Elith and Graham et al. (2006) concluded that the MaxEnt model simulated the best results[50]. Meanwhile, Prevéy, J. S. et al. (2020) used the MaxEnt model to predict the potential fitness zone and phenology of Vaccinium membranaceum (American lingonberry) and found that the simulated results were in good agreement with observations[51]. The MaxEnt model is an ecological niche model based on the theory of maximum entropy[52], and is the preferred model for predicting the potential distributions of critically endangered species because of its simplicity, high degree of accuracy, low sample size requirement, and high stability compared to other modeling methods [53,54,55].

Ilex nanchuanensis is a species of holly tree endemic to China and is classified as critically endangered in the IUCN’s Red List of Threatened Species , standard conservation level. It is mainly found on Jinfo Mountain (the highest peak of the Dalou Mountains in Nanchuan District, Chongqing, China) at altitudes of 600–800 m. It has high economic and ecological value and plays an important role in the ecological balance and water conservation of the region; plus, it is an important indicator species for the state of the ecological environment[56,57]. To date, little research has been reported on I. nanchuanensis, and in particular, its response to climate change remains unclear. Therefore, using the MaxEnt model to simulate the response of the geographic distribution of I. nanchuanensis to climate change will not only help to clarify the causes of habitat loss and potential species extinction, but also the effects of strengthened climate change on its potential geographic distribution, thus providing a theoretical and practical basis for conservation management of this species as well as other rare and endangered species.

In this study, based on current distribution data for I. nanchuanensis and climatic factors in different periods, the MaxEnt model, combined with ArcGIS software, was used to simulate the potential habitat of the I. nanchuanensis during the Last Interglacial (LIG), the Last Glacial Maximum (LGM), the current period, and two future periods (2050s and 2070s), to clarify the response of this cherished and critically endangered species to climate change, as well as provide a theoretical and practical basis for conservation management of not only this species but also other rare and endangered species throughout the world living in fragile habitats.

2. Materials and methods

2.1. Species distribution data

Specimen record points were obtained by searching the literature, resource survey reports, and herbaria (e.g., the China Virtual Herbarium, Global Biodiversity Information Facility database, China Biological Species List 2022, Flora of China, Plant Science Data Center)[45,58], removing duplicate records and sample record points lacking geographic information, and using Google Maps software to 5 km. The spatially overlapping data points within the range were removed and the coordinate system was converted to latitude and longitude information, and 34 samples of I. nanchuanensis were obtained with clear geographic information centered on Nanchuan District, Chongqing, China [46,59].

2.2. Bioclimatic data acquisition and screening

The bioclimatic data were obtained from the WorldClim database (http://www.worldclim.com (accessed on Oct. 25, 2022)) at a resolution of 2.5 arc-min[60,61]. The data contain four sections: the LIG period (120,000–14,000 years ago), the LGM period (22,000 years ago), the current period (1970–2000), and two future periods (the 2050s and 2070s). They include 19 climatic factors, such as monthly precipitation, monthly mean temperature, monthly maximum temperature, and monthly minimum temperature (Figure 1). To avoid the influence of correlation between climatic factors on the simulation results, according to the climatic factor screening method of Yang and Kushwaha et al. (2013), multicollinearity analysis was used in SPSS19.0 to detect the correlation between variables[62]. Environmental variables with correlation values <0.8 were retained, and their contribution rates were analyzed. Six climatic factors were finally screened to construct the model (Table 1): seasonal temperature variation (Bio4); hottest monthly minimum temperature (Bio6); annual mean temperature range (Bio7); wettest monthly precipitation (Bio13); seasonal precipitation (Bio15); and driest quarterly precipitation (Bio17).

With reference to the IPCC’s assessment reports, four greenhouse gas emission scenarios (Representative Concentration Pathways, RCPs) were employed: the low emission (RCP2.6) scenario; two moderate emission scenarios (RCP4.5 and RCP6), and the high emission (RCP8.5) scenario [63,64]. Among the various state-of-the-art climate models, it was found that BCC-CSM1.1 (Beijing Climate Center, Climate System Model, version 1.1) was more consistent with the characteristics of climate change in China, while MIROC-ESM (Model for Interdisciplinary Research on Climate, Earth System Model) and CCSM4 (Community Climate System Model, version 4) were used for the LGM period [42].

2.3. MaxEnt model operation and evaluation

Geographical distribution data (in CSV format) and bioclimatic data (in ASC format) were imported into the MaxEnt model for calculation, while the bioclimatic data needed to be modeled separately with the geographical distribution data of the species according to the different periods. Following the method of Moreno et al. (2011), the bioclimatic factors in the model were set as continuous variables, and then bioclimatic data for the LGM, LIG, and future periods (2050s and 2070s) were input in turn to analyze the geographic distribution of I. nanchuanensis for each of these periods[65]. Response curves were constructed by selecting the ‘Create response curves’ option. To ensure the accuracy and objectivity of the data, the MaxEnt model should be constructed by randomly selecting 75% of the species’ geographic distribution as training values, while the remaining 25% of the data should be tested by 100 replications to test the model and create ROC (receiver operating characteristic) curves. The remaining model parameters were operated according to the default values of the software[66].

The results of any model simulation need to be evaluated comprehensively. After the MaxEnt model was established, the area under the ROC curve was selected for evaluation in order to verify the accuracy of the evaluation results, and the AUC (area under curve) value was used as an important indicator to measure the predictive ability of the evaluation model. Note that, because the AUC value is the MaxEnt model’s own procedure, the simulation process does not require additional thresholds, is involved in the assessment of simulation predictions of the potential geographic distribution of numerous species, and has a high degree of accuracy[67,68]. The AUC takes values in the range [0, 1], and values closer to 1 prove that the model’s results are more accurate [69]. In general, when the AUC value is >0.7, it proves that the model’s simulation results are good, as well as indicating potentially suitable geographic distribution areas for the species[70].

2.4. Classification of habitat suitability

The results of the MaxEnt model simulations were collated and analyzed based on ArcGIS software and converted to Raster format, and the spatial analysis tool was used to classify and visualize the species’ suitability zones and to obtain a map of the potentially suitable geographic distribution zones of I. nanchuanensis in different historical periods. The results of the MaxEnt model were interpreted in the form of the fitness probability P of the species occurring in the area to be tested, with P taking values in the range [0, 1]. In addition, based on the actual distribution of I. nanchuanensis and the methodology reported by the IPCC, the suitable areas for the species were classified into four types: non-suitable (P < 0.049), low suitability (0.049 ≤ P < 0.196), medium suitability (0.196 ≤ P < 0.441), and high suitability (0.441 ≤ P < 0.833). Finally, the areas of different fitness classes of I. nanchuanensis in the different historical periods were calculated using the area statistics function in the spatial analysis tool[71].

2.5. Analysis of centroid migration in suitable distribution areas

The center of mass is the equilibrium point of moment for species in a region during a historical period, and can effectively represent the spatial aggregation and migration of species over a relative period of time. The Toolbox software in the SDM can capture the location of the species’ center of mass and analyze the changes in the migration of the species’ center of mass in different historical periods by combining geographical information data such as latitude and longitude[72,73]. The methods were as follows: (1) The results of the MaxEnt model simulating the potential geographic distribution of I. nanchuanensis were transformed into a vector binary, and the suitability of the suitability zone was defined according to the suitability probability P, with P ≥ 0.5 being the total suitability zone and P < 0.5 being the non-suitable zone. (2) Using the Zonal Geometry statistical tool in the Spatial Analysis Tool, the plasmodesmata were selected in the Geometry Type option to capture the specific coordinates of the placenames in the species-suitable zone that were captured. (3) The placenames of I. nanchuanensis were analyzed at different times and under different climatic conditions, and the migration direction and distance of I. nanchuanensis were finally assessed[46].

3. Results

3.1. Accuracy of model analysis

In order to assess the accuracy of the MaxEnt model, the area under the ROC curve, i.e., the AUC value, was selected as an indicator to measure and assess the predictive ability of the model, wherein the closer the AUC value is to 1, the better the simulation results are[53]. In this study, among the simulation results for the LIG period, the LGM period, and the current period, the AUC value was 0.995 (Figure 2), indicating that the accuracy of the MaxEnt model simulation in this study was extremely high and the confidence level was excellent. It also shows that both CCSM4 and MIROC-ESM are optimal in the LGM period.

3.2. Dominant environmental factors

Environmental variables influence the potential distribution of the species at different times. As can be seen from Table 2, there are six environmental variables that are most important in influencing the potential distribution of I. nanchuanensis in different historical periods: Bio4 (standard deviation of seasonal variation in temperature), Bio6 (minimum temperature in the hottest month), Bio7 (annual mean temperature range), Bio13 (precipitation in the wettest month), Bio15 (seasonal coefficient of variation in precipitation), and Bio17 (precipitation in the driest season).

Using Jackknife to test the percentage of regularized training gain of environmental variables on the potential distribution of I. nanchuanensis, we found (Figure 3) that the environmental factors in the LIG period, LGM period (CCSM4 and MIROC-ESM), and current period had the same weighting, from highest to lowest: Bio6, Bio15, Bio17, Bio4, Bio7, Bio13.

Analysis of the probability response curves of environmental variables and species distributions revealed that Bio6, Bio15, Bio17, and Bio7, all four environmental factors, had single-peaked, one-way response curves with probabilities of the presence of major environmental factors >0.5, which remained largely consistent with the range of environmental factor parameters for highly suitable habitats. The analysis showed (Figure 4) that the optimum range for Bio6 was −30°C to 25°C, for Bio15 it was 50–65 mm, for Bio17 it was 70–130 mm, and for Bio7 it was 230°C–310°C. According to the Jackknife test environmental variable method of simulation detection, the environmental variable that achieves the highest formalized training gain when only one variable is present during the simulation is Bio6. Therefore, Bio6 has the most effective simulation information.

3.3. Suitable species distribution areas

The MaxEnt model was used to simulate the distribution of I. nanchuanensis during the LIG period, LGM period (CCSM4 and MIROC-ESM), and current period, and it was found that its distribution range varied significantly during different periods (Figure 5). The potential distribution of I. nanchuanensis was classified into four classes using the natural break method (grey for unsuitable, green for marginally suitable, yellow for moderately suitable, and red for highly suitable distribution areas). This species of holly is endemic to China, but the large map presented an incomplete map of China, so a small map was used to refine the U-shaped line. As shown in Figure 5D, I. nanchuanensis is currently distributed in China mainly in Nanchuan District, Wulong District, Chongqing Municipality, Wuping Town, Fengdu County, Chongqing Municipality, Nanchuan City, Yichang City, and Wufeng City, Hubei Province, mainly between 29.00728°N and 30.17125°N, with the highly suitable area being mainly in Nanchuan District, Chongqing Municipality.

The MaxEnt model simulation revealed that the distribution range of I. nanchuanensis can be better captured at present. The total suitable distribution area was 37.722223×10⁴ km² ; the highly suitable distribution area was 5.708333×10⁴ km²; the moderately suitable distribution area was 11.98785×10⁴ km²; and the marginally suitable distribution area was 20.02604×104 km². As shown in Table 3, comparison with the current distribution area of the species reveals that the highly suitable area decreased by 14.317707×10⁴ km² and 6.279517×10⁴ km² compared to the marginally suitable area and moderately suitable area, respectively, while the total suitable area of the species increased compared to the LIG period, CCSM4 LGM period, and MIROC-ESM LGM period, with an increase of 2.348959×10⁴ km². The total area suitable for this species has increased by 2.348959×10⁴ km², 5.373267×10⁴ km² and 3.309025×10⁴ km², respectively, and it can be seen from Figure 5 that the area suitable for each distribution class of I. nanchuanensis has shown an increasing trend compared to the LIG, the CCSM4 LGM, and MIROC-ESM LGM periods.

The geographical distribution of I. nanchuanensis in two future period (the 2050s and 2070s) was also studied using the MaxEnt model, combined with four climate scenarios (RCP2.6, RCP4.5, RCP6, and RCP8.5) under global climate change to simulate the potential geographical distribution range of I. nanchuanensis and its suitable distribution area for different classes (Table 3, Figure 6). It was found that the potential distribution of I. nanchuanensis will continue to shrink in the future, especially in the 2070s, when the total suitable area, highly suitable area, and moderately suitable area will be reduced to varying degrees compared to the current one, especially under RCP6, with a maximum reduction of 9.04688×10⁴ km², and under RCP2.6, with a minimum reduction of 4.354159×10⁴ km². The maximum reduction in distribution area under RCP6 was 17.779519×10⁴ km², and the minimum reduction under RCP4.5 was 3.937494×10⁴ km². In the highly suitable distribution area, a maximum reduction of 8.055556×10⁴ km² was found under RCP6, a minimum reduction of 1.286459×10⁴ km² under RCP2.6, and a small increase under RCP4.5. In the moderately suitable distribution area, a maximum reduction of 7.390623×10⁴ km² was found under RCP6, and a minimum reduction of 1.286459×10⁴ km² under RCP8.5. In the marginally suitable distribution area, a maximum reduction of 7.390623×10⁴ km² was found under RCP6, and a minimum reduction of 0.9375×104 km² under RCP8.5.

3.4. Suitable distribution under future climate scenarios

The MaxEnt model was used to predict the distribution of I. nanchuanensis under different climate scenarios for the future periods, and the area and distribution of I. nanchuanensis in unsuitable, unchanged, expanded, and contracted areas can be seen in Table 4 and Figure 6. The area in unchanged areas will not change much compared to the current period, while the area in expanded areas will increase. Compared to the potential distribution area of I. nanchuanensis in the current period, its distribution range in Sichuan, Chongqing, Hubei, Guizhou, southern Shaanxi, Jiangsu, and Fujian in China will increase in the 2050s (Figure 6A–D), while it will decrease significantly in the 2070s (Figure 6E, F).

3.5. Migration trends of the geometric center of suitable habitat

The SDM-Toolbox in ArcGIS 10.6 software was used to analyze the distribution area of I. nanchuanensis in the current period and under different climate scenarios in the future periods to obtain the mass migration trajectory and direction (Figure 7). The distribution area of I. nanchuanensis in the current period is located in Nanchuan District, Chongqing, China (29.12382309°N, 107.2505951°E), while in the LIG and LGM periods it is still in Nanchuan District but at (29.02701706°N, 107.146697°E) and (29.02930604°N, 107.164464°E), respectively (Figure 6a). Although the distribution of I. nanchuanensis in the above three periods was in Nanchuan District, Chongqing, in the current period it showed a trend of migration to the northeast compared with that in the LIG and LGM periods (59229 m and 54900 m, respectively) (Figure 6b). It was found that the distribution of I. nanchuanensis in China in the future (2050s and 2070s) will migrate to the northwest compared to the current period. By simulating the migration trajectory of I. nanchuanensis placenames under different climate scenarios in the future periods, it was found that it migrated 86072 m, 122,762 m, 125,080 m and 166,544 m under RCP 2.6, RCP 4.5, RCP 6 and RCP8.5, respectively, in the 2050s compared to the current period, and in the 2070s by 94,894 m, 101,163 m, 153,836 m and 150,276 m (Figure 6B).

4. Discussion

4.1. Response of species distribution to climate change

Climate change is a major driver of species distributions[74,75], and species distributions can also reflect climate change characteristics[76,77,78]. Species change seasonal activity patterns, abundance, and migration in response to challenges such as extreme events, heat stress, and water scarcity brought about by climate change[79]. Studies have found that global climate change is already causing species to migrate to higher altitudes and latitudes[80,81]. Ilex nanchuanensis is a critically endangered species endemic to China, with only one distribution area of <100 km², and has experienced habitat degradation or even loss due to climate warming[16,17]. In the present study, I. nanchuanensis was found to be migrating northwards in both the current and future periods, demonstrating the stress of climate warming imposed on this species, with the main climatic factors affecting its geographical distribution found to be temperature and precipitation. Similarly, Ning et al. (2021) found that temperature and precipitation were the most important environmental factors affecting the geographical distribution of Pinus armandii[75], while You et al. (2018)found that temperature-related variables affected the geographical distribution of Rhodiola more profoundly than precipitation[82].

The use of the MaxEnt model to simulate the distribution ranges of species in different historical periods based on current species distribution data and climatic factors can help maintain ecosystem balance and protect rare and critically endangered species[47,83]. It was found that current global temperatures have increased by 5°C–12°C compared to the LGM period and that the current glacial area in China has decreased by 8.4 times compared to the LGM period, and these changes have affected the geographical distribution patterns of species[84]. In the present study, we found that the distribution range of I. nanchuanensis in the current period was larger than that of both the LGM period and LIG period, which was related to the species’ biological characteristics such as its preference for warmth and tolerance of shade. Willis et al. (2004) found that mid-latitude and high-latitude regions were affected by glacial cover during the Quaternary period, while China was located at a relatively low latitude and the complex topography and mountains were able to effectively block the glacial cover[85,86], confirming the findings of this study that I. nanchuanensis migrated to the northeast during the current period[87]. Cornejo-Romero et al. (2017) found that the range of drought-tolerant plants was less during the current period than during the LIG, which is contrary to the results of this study owing to the fact that the biology of drought-tolerant plants is very different from that of subtropical plants[88].

4.2. Potential distribution area changes and conservation of Ilex nanchuanensis under future climate scenarios

Climate change is not only forcing species to migrate but also accelerating the risk of extinction[89]. In this study, we found that the distribution range of I. nanchuanensis in the future periods (2050s and 2070s) is likely to be reduced compared to the current period, and it keeps moving towards the northwest owing to extreme phenomena such as local increases in temperature and water shortages[80,81]. In addition, the distribution range of I. nanchuanensis is projected to shrink in the future periods, especially the 2070s, with the total suitable area being reduced the most under RCP6 and the least under RCP4.5, which is similar to the findings of other studies on rare and endangered species in China, such as ginkgo (Ginkgo biloba) and Taiwan fir [90,91]. This is closely related to the phenomenon of climate warming in the future, with studies finding a global temperature increase of 1.3°C by 2050 and a minimum global temperature increase of 2.2°C by 2070[5,92]. This study found that the optimum range for the hottest monthly minimum temperature (Bio6), the most important climatic factor affecting I. nanchuanensis, is −30°C to 25°C, and that extreme weather and global warming have made it difficult for the species to survive better locally, which is one of the main reasons for the migration of I. nanchuanensis to the northwest in the coming period. Also, I. nanchuanensis migrated the furthest under RCP6 in the 2070s, confirming the threat of rising temperatures and increased greenhouse gas emissions to the distribution status of the species, as well as the global ecosystem[93]. Therefore, there is a need to reduce greenhouse gas emissions to protect global species diversity and ecosystem balance. This study also found that the best ranges of climatic factors affecting the geographical distribution of I. nanchuanensis were 70–130 mm for the driest seasonal precipitation (Bio17) and 230°C–310°C for the mean annual temperature range (Bio7). Meanwhile, I. nanchuanensis not only prefers a warm climate but also has a certain degree of cold tolerance, which is another reason for its migration to the northwest in the coming period. To better protect the species, research on drought stress and salt stress should be a focus. Therefore, we recommend the construction of a core reserve in an area centered on Nanchuan District of Chongqing, China, to improve habitat conditions and the conservation of germplasm resources using techniques such as hybrid seeding and genetics to increase the genetic diversity of wild populations.

5. Conclusions

This study simulated the geographical distribution of I. nanchuanensis using the MaxEnt model, and the following conclusions were obtained: (1) The climatic factors affecting the geographical distribution of I. nanchuanensis are Bio4, Bio6, Bio7, Bio13, Bio15, and Bio17 (see Table 1 for definitions). (2) Ilex nanchuanensis migrated to the northeast in historical periods but is projected to migrate northwest in the future. (3) The total suitable area for I. nanchuanensis was projected to reduce the most under RCP6 in the 2070s (17.779519×10⁴ km²), and the least (3.937494×10⁴ km²) under RCP4.5. To address the threat of global climate change to I. nanchuanensis, a protected area centered on Nanchuan District in Chongqing, China, should be established, and germplasm resources protected and cultivated, with core, buffer, and experimental areas scientifically and rationally delineated. In addition, greenhouse gas emissions need to be reduced, and future greenhouse gas emissions need to be controlled within the scenario of RCP4.5.