Climate change affects the suitability of Chinese cherry (Prunus pseudocerasus Lindl.) in China

Abstract

The rapid development of Prunus pseudocerasus related industry has increasingly contributed to rural vitalization in China. This study employed a biomod2 ensemble model, utilizing environmental and species occurrence data from 151 P. pseudocerasus germplasm wild/local samples, to predict potential geographical distribution, suitability changes, climate dependence, and ecological niche dynamics. The optimized maximum entropy (MaxEnt) model yielded the most accurate predictions. The climate variables with the greatest impact on suitability were precipitation of warmest quarter and mean diurnal temperature range. The total potential suitable area for P. pseudocerasus was approximately 2.78 × 10⁶ km², increasing with CO₂ concentration. The highly suitable area was primarily concentrated in basin terrains, plateaus, and plains of Sichuan Province. The current centroid in Lichuan exhibited gradual latitudinal and longitudinal movement. The predicted (2090s) ecological niche trends of P. pseudocerasus varied under different pathways and periods, with higher CO₂ concentration associated with lower niche overlap. The CO₂ emission concentration in the SSP246 scenario emerged as the most suitable climate model. Climate change is driving both the expansion of geographical distribution and the contraction of overlapping geographical distribution areas of P. pseudocerasus. These findings provide a theoretical basis for wild resource conservation, site selection for production, and introduction of allopatry for P. pseudocerasus.

Core

The climate variables with the most significant impact on suitability were precipitation of warmest quarter and mean diurnal temperature range. The total potential suitable area for P. pseudocerasus was approximately 2.78 × 10⁶ km², exhibiting an increase in correlation with CO₂ concentration. The highly suitable areas were predominantly concentrated in the basin terrains, plateaus, and plains of Sichuan Province. Climate change is driving both the expansion of geographical distribution and the contraction of overlapping geographical distribution areas of P. pseudocerasus.

Introduction

China has abundant germplasm resources of wild Cerasus species worldwide (Wu et al. 2003; Jiu et al. 2023, 2024a; Wang et al. 2023). Chinese cherry (Prunus pseudocerasus) originated in the middle and lower reaches of the Yangtze River in China and is distributed throughout the Yangtze River Basin. It is extensively cultivated in Sichuan, Guizhou, Chongqing, Henan, Shandong, Anhui, and Gansu provinces, particularly on sun-exposed mountain slopes and ravine sides (Chen et al. 2013). According to the Flora of China, P. pseudocerasus has been cultivated in China for over two millennia, making it one of the most significant and ancient cultivated fruit trees (Wu et al. 2003). P. pseudocerasus is highly regarded by consumers for its aesthetically pleasing tree shape, pleasant floral fragrance, vibrant fruit color, and organic nature. Consequently, it plays a crucial role in China's rural tourism industry (Chen et al. 2015, 2016).

Limited research has been conducted on the geographical distribution and ecological characteristics of these plants to enhance their production (Zhu et al. 2017, 2019). As a deciduous fruit tree, temperature plays a critical role in determining the dormancy transition of P. pseudocerasus. Consequently, climate change significantly impacts its distribution, and the expansion of suitable planting areas will substantially influence producers' income and industry development (Olsen, 2010). In recent years, with the growth of agricultural tourism and agriculture, P. pseudocerasus has garnered increasing attention due to its high ornamental value, early phenological period (flowering from February to March, and maturation from April to May), strong adaptability, and disease resistance. Consequently, the area under cherry cultivation has progressively expanded (Zhang et al. 2018; Jiu et al. 2024b).

The relationships between plants' geographical distribution patterns, climate change, and predictions of distribution changes are crucial for understanding basic ecology and biogeography, and have become a prominent research focus (Elith et al. 2009). Globally, plants’ responses to hydrothermal conditions can effectively predict their geographical distribution limits. Historical climatic conditions have shaped the current distribution pattern of P. pseudocerasus, which continues to be influenced by ongoing climate change. Species distributions at a regional scale are primarily affected by climate, with a combination of hydrothermal variables reflecting fundamental climatic characteristics (Woodward et al. 1987; Thuiller et al. 2011; Nadine et al. 2012). Various climatic changes exert different promoting or inhibitory effects on plants, resulting in corresponding variations in suitable areas (Camallec et al. 2024). The global climate is undergoing significant changes, substantially impacting plant geographical distribution (Kristine et al. 2013). By the end of the twenty-first century, global average surface temperature is predicted to increase by 0.3 to 4.8 ℃, accompanied by significant changes in precipitation patterns. Such dramatic climate changes may alter future geographical distribution patterns of species, exacerbate biodiversity reduction and germplasm resource loss, and potentially accelerate species extinction rates (Céline et al. 2012). Climate change significantly affects the areas, growth, development, and yield of suitable crops, as well as the prevalence of diseases and pests, and overall plant systems (Wu et al. 2023).

Mounting evidence suggests that climate warming has significantly altered agricultural climate resources in China, leading to shifts in the geographical distribution and spatial patterns of numerous crops. In Japan, climate change has modified the phenology of flowering cherry populations. Research has explored the impact of climate change on the initial four phenological stages of sweet and sour cherry trees grafted on seedlings in a temperate-continental climate. These stages include bud swelling, budburst, beginning of flowering, and end of flowering. The study also analyzed the implications for farmers. The researchers hypothesized that the onset of these phenological stages occurs earlier due to climate change (Yang et al. 2020; Primack et al. 2009; Paltineanu et al. 2020).

Climate change has significantly impacted species distribution and abundance in recent years (Thuiller et al. 2005; Alexander et al. 2015; Rumpf et al. 2019). To address research objectives, numerous scholars have introduced species distribution modeling methods (Phillips et al. 2006; Chen et al. 2011). These digital models employ algorithms to quantitatively analyze potential distribution areas or specific requirements of target species, utilizing known species distribution information and environmental variables (Phillips et al. 2006). Species distribution models (SDMs) are constructed based on specific algorithms to predict possible species distribution areas or ecological conditions of particular requirements. These models have been widely applied in various fields, including predictions of suitable area patterns and pests, as well as biodiversity monitoring (Kumar et al. 2014; Lake et al. 2020; Fang et al. 2021). SDMs have been extensively utilized in studies of species distribution, ecology, biodiversity, climate change, conservation biology, and horticultural crops (Wan et al. 2014, 2016; Duan et al. 2016; Searcy et al. 2016). For instance, Rivera-Parra et al. (2020) simulated suitable distribution areas for high-quality tea trees in Ecuador, identifying appropriate regions for Ceylon and Nilgiris. Li et al. (2020) employed the maximum entropy (MaxEnt) model and high-performance liquid chromatography to comprehensively analyze the suitability of three Coptis chinensis species in China. Their research indicates that the annual precipitation range and isothermality are the most significant environmental variables affecting the suitable distribution of Chinese goldthread and tea.

The ecological niche is likely a primary determinant of changes in habitat selection and spatial distribution patterns of species. Niche width indicates the range of ecological adaptation under varying carbon dioxide (CO₂) concentration pathways and periods for P. pseudocerasus. Niche overlap can elucidate the similarity in ecological requirements under different CO₂ concentrations and periods for P. pseudocerasus. Consequently, ecological niche analysis can predict the impact of climate change on community succession in P. pseudocerasus (Lasky et al. 2015). While the low information density of finite data points impedes the acquisition of high-precision spatial distribution information on temperature and precipitation, appropriate interpolation methods can effectively evaluate discrete site data to obtain precise spatial distribution information of regional precipitation (Zhang et al. 2010).

Employing the P. pseudocerasus fruit tree as the subject of investigation, this study utilized the MaxEnt model and a geographic information system (GIS) to simulate the suitable distribution of the tree across China. The research forecasted changes in the suitable area, area of occupancy, ecological niche, and distribution range of P. pseudocerasus based on climate data from various temporal periods. These findings offer a theoretical foundation for identifying suitable cherry cultivation areas and provide data-driven support for wild resource conservation, cultivation practices, agricultural industry distribution, refined planting management, and strategies to address future climate change impacts.

Results

Analysis of current geographical distribution and evaluation of model accuracy

The geographical distribution (151 sample points) of P. pseudocerasus is depicted in Fig. 1a and Table S1. Based on specimen and literature records, P. pseudocerasus primarily occurs in the Chinese provinces of Sichuan, Shandong, Yunnan, Zhejiang, and Guizhou. Its horizontal geographical distribution spans from 23.448°–37.356°N and 101.540°–121.828°E. The species' range extends from Yantai City, Shandong Province, in the east and north, to Liangshan Yi Autonomous Prefecture, Sichuan Province, in the west, and Mengzi City, Yunnan Province, in the south. P. pseudocerasus exhibits a substantial vertical distribution range (3–3488 m), with 58 collection sites recorded at altitudes exceeding 1000 m (Fig. S1a). The slope range varies from 0.007°–6.855°, and the aspect ranges from 2.251–357.022 (Fig. S1b, c).

Fig. 1

Geographical distribution evaluation indices of individual predictive models and accuracy. a Geographical distribution and local P. pseudocerasus resources. b AUC model accuracy test of eight climate models. c TSS model accuracy test of eight climate models. d MaxEnt model accuracy tests by fractional value. e MaxEnt model accuracy tests by sensitivity. f–i Linear regression analysis of Shandong temperature (f), Sichuan temperature (g), Shandong precipitation (h), and Sichuan precipitation (i). j Correlation analysis of model climate data and observed climate data. k Classification labels of the confusion matrix. l–o Confusion matrices of Shandong temperature (l), Shandong precipitation (m), Sichuan temperature (n), and precipitation (o). In b–c, the grey dots represent the data for each distribution. In j, ** indicates significance of the correlation at the 0.01 level (2-tailed). In f–o, temperature and precipitation data are calculated as monthly average values, with red and blue lines representing the fitted curves of observed and predicted data, respectively

The model accuracy test results, based on the area under the receiver operating characteristic (ROC) curve (AUC) and total sum of squares (TSS), are presented in Fig. 1b, c and Table S2. MaxEnt exhibited the highest values for both metrics (0.890 and 0.915), leading to its selection for modeling. Following the modeling process, the mean AUC value of the test data was 0.939 (Fig. 1d, e). This value exceeded the random test AUC (0.500), indicating that the ROC curve demonstrated good predictive performance for cherry distribution. To validate the accuracy of the WorldClim data and MaxEnt model, we compared real observed data with model-predicted data for Sichuan Province and Shandong Province, which had the highest number of records (35 and 30, respectively). Observed data from 1979–2015 included monthly average temperature and precipitation. We employed three analysis methods (linear fitting, correlation analysis, and confusion matrix) to verify the accuracy of these two climate indicators between observed and predicted data. Linear regression analysis revealed similar temperature data for Shandong and Sichuan provinces (Fig. 1f, g), while precipitation data showed low similarity between the two provinces (Fig. 1h, i). Correlation analysis for Shandong and Sichuan provinces indicated that Pearson's, Spearman's, and Kendall's tau-b correlations were all significant at the 0.01 level (Fig. 1j). The hierarchical labels of the confusion matrix are shown in Fig. 1k. The confusion matrix demonstrated high accuracy (95.71%) for Shandong Province's temperature data (Fig. 1l and Table S3), whereas precipitation data accuracy was only 45.84% (Fig. 1m). Both temperature and precipitation data for Sichuan Province exhibited low accuracy (Fig. 1n, o).

Dominant environmental variables limiting P. pseudocerasus distribution

We screened the climate data from the distribution records and combined the principal component analysis (PCA) results to verify the data's conformity to a normal distribution. Components with eigenvalues > 1 were extracted. As shown in Table S4, the first six principal components summarized most of the original data. According to the eigenvalue and contribution rate of each eigenvector, the first principal component had the highest contribution rate, determined by Bio12, Bio1, and Bio6; the main variables were associated with hydrothermal interactions. The second principal component's contribution rate was 26.49%, determined by Bio4, Bio7, and Bio10, with the main variables associated with temperature. The third principal component's contribution rate was 10.34%, determined by Bio8, Bio1, and Bio19, also primarily associated with temperature. Three-dimensional loading plots (Fig. 2a) indicated that the environmental variables in different periods were almost on the same confidence ellipse. To avoid model overfitting, using the contribution rate of each variable obtained by PCA as one basis for screening environmental variables, we performed Spearman rank correlation testing for all data (Fig. 2b). After excluding highly correlated variables, eight variables were selected for modeling (Bio2, Bio8, Bio14, Bio15, Bio18, elevation, slope, and aspect). Furthermore, we employed Jackknife's test to assess the relative importance of the selected eight variables (Fig. 2c). The environmental variables with the highest contribution rates were Bio18, Bio2, and Bio14. These three environmental variables likely represent the dominant factors affecting the suitable habitat distribution of P. pseudocerasus.

Fig. 2

Climatic variable selection and trend analysis. a Principal component three-dimensional load diagram. b Correlation matrix heat map of 22 environment variables. c Jackknife test of the relative importance of environmental variables for P. pseudocerasus in China. d Response curves of Bio2, Bio8, Bio14, Bio15, and Bio18. e Changes in important climate variables for SSP126-2090s and SSP585-2090s. In a, each ball in the three-dimensional coordinate system represents an independent sample, with different colors denoting different environment variables. In i-n, the dashed line indicates outliers with a factor of 1, the box range represents standard deviation, the solid line in the middle of the box denotes the median, the white squares indicate averages, and the grey dots represent the data for each distribution

The response curves depicted in Fig. 2d-h illustrate the relationship between predicted suitability and selected variables, as well as the interdependencies induced by correlations among variables. The optimal climatic conditions for P. pseudocerasus were identified as follows: 5–7℃ (Bio2), 20–25℃ (Bio8), 10–13 mm (Bio14), 60–80 mm (Bio15), and 480–520 mm (Bio18). Notably, all climate variables except Bio2 exhibited a similar pattern characterized by an initial increase followed by a decrease. Bio2 displayed a more complex trend. The suitability peaked when Bio2 reached 6℃, decreased between 6–8℃, slightly increased at 10.3℃, and then rapidly declined, approaching zero. Additionally, we examined the five climate variables used in the modeling under two extreme scenarios for the 2090s. Given the significant role of Bio1 in biological characteristics, we also analyzed its changes (Fig. 2i). Most variables, except Bio14, showed an upward trend, albeit with varying patterns (Fig. 2j-n). Further analysis of climate variable changes at each data point revealed some variations (Table S5). Generally, the following trends were observed: Bio8 increased with CO₂ concentration, except in Gejiawu and Fengtong villages in Zhejiang Province. Bio14 decreased with CO₂ concentration, but increased in Datun village (Sichuan Province), Lanzi village (Shandong Province), and Yujiahe village (Jiangxi Province). Bio18 increased with CO₂ concentration, except in Honghuagang district and Yongyi village (Guizhou Province), and Lujiao town (Chongqing City).

Predicted expansion and northeast shift of suitable areas of P. pseudocerasus

The current potential distribution area of P. pseudocerasus encompasses Liaoning, Shandong, Tianjin, Henan, Anhui, Sichuan, Guizhou, Yunnan, Shaanxi, Hunan, Hubei, Gansu, Guangxi, and Tibet. This range extends beyond that described in the Flora of China (Fig. 3a). Figure 3b illustrates the changes in suitable areas across different regions during the 2050s, 2070s, and 2090s. The suitable areas consistently expand over time (Table S6). To provide a clearer representation of future changes in P. pseudocerasus suitability, we analyzed changes in suitable projections for provinces, cities, and autonomous regions (Fig. 3c and Table S7). Currently, the total suitable projected area is 2.78 × 10⁶ km². Yunnan, Sichuan, and Hunan provinces possess the largest total suitable projected areas, while Inner Mongolia has the smallest, at zero. Across different pathways, Yunnan Province exhibited the largest suitable projected area in the 2050s, 2070s, and 2090s. However, Sichuan Province consistently demonstrated the largest high suitable projected area across all stages and pathways. Moreover, nearly all regions displayed three levels of suitability change in each period. Hunan, Qinghai, Inner Mongolia, and Hubei provinces showed insufficient areas of increased suitability, while areas of decreased and unchanged suitability were substantial. Hong Kong, Macao, and Taiwan provinces exhibited no suitability changes, although Taiwan Province experienced slight suitability decreases and remained unchanged near the north–south section of the Taiwan Mountains (Fig. 3b). Across all periods, decreased and unchanged suitabilities primarily concentrated in two regions: Shandong, northern Jiangsu, eastern Henan, and the Sichuan Basin, Yunnan Guizhou Plateau. Calculation of the sum of the three types of suitable areas (high, medium, and low) indicated that the SSP246 CO₂ emission concentration model was the most suitable climate model for P. pseudocerasus.

Fig. 3

Potential distribution of the P. pseudocerasus cultivars. a Current potential distribution of P. pseudocerasus. b Potential distribution of P. pseudocerasus across various periods and SSPs pathways. c Heat map illustrating projected suitable areas for each province

The study compared potential distribution patterns and analyzed their spatial transformation under various future scenarios (Fig. 4a). The largest suitable gain areas were observed in the shared socioeconomic pathway (SSP585-2090s) and SSP245-2070s (Fig. 4b and Table S8). These gain rates were primarily concentrated in eastern Heilongjiang Province and southwestern Gansu Province, respectively. Conversely, SSP370-2090s and SSP585-2090s predicted the most significant loss of suitable areas, mainly in eastern Yunnan, central Guizhou, and southeastern Sichuan. The most extensive stable areas were found in SSP370-2050s and SSP126-2050s (83.49% and 82.20%, respectively), with relatively wide and fixed distributions (Fig. 4a, c). Notably, the change in suitable areas correlated with CO₂ concentration, following a similar trend across all periods. However, the proportion of suitability change areas in each period gradually increased, suggesting that P. pseudocerasus may have some capacity to adapt to global warming. These suitability changes resulted in significant geographical variations. The changed areas were dispersed throughout the region, while stable areas remained relatively fixed, located primarily east of the Taihang Mountains, around the Qinling Mountains, and between the Hengduan, Daba, and Wuling Mountains (Fig. 4a).

Fig. 4

Suitable distribution and gravity center movement trajectory of the highly suitable zone under different SSPs pathways and periods. a Spatial distributions of suitability changes. b Heat map of projected area changes. c Proportions of gain, stability, and loss. d Gravity center movement trajectory of the highly suitable zone. In d, the arrow color indicates the magnitude of the distance (map data from http://bzdt.ch.mnr.gov.cn/)

Utilizing ArcGIS spatial analysis capabilities and MaxEnt simulation results, this study examined the migration patterns of the centroid (geometric center) of suitable cherry areas in China under climate change scenarios. The centroid of the highly suitable area was identified in Lichuan City (Enshi Tu and Miao autonomous prefecture, Hubei Province), near the Hubei-Chongqing border (Fig. 4d). The center of gravity trajectory varied across different CO₂ emissions pathways and time periods. The most significant displacement occurred from the current position to SSP585-2050s (215.14 km), with the center of gravity shifting northeast to Zhuxi County (Shiyan City, Hubei Province). Despite variations in the moving trajectory under different pathways and periods, the predominant direction of movement was consistently northeast. Notably, the trajectories for SSP245-2050s–SSP245-2070s and SSP245-2070s–SSP245-2090s exhibited unique patterns, even moving contrary to other trajectories, possibly due to irregular climate changes during these intervals. The longitude and latitude ranges were 30.319°N–33.537°N and 108.515°E–112.046°E, respectively, with the center of gravity traversing Hubei, Chongqing, and Henan provinces. In conclusion, the analysis indicates a general trend of the center of gravity moving towards higher latitudes and longitudes (Table S9).

Climate dependence of P. pseudocerasus suitability

Climate change impacts the suitable areas for species distribution, with its influence being multifaceted, encompassing both climatic and species-specific factors. Given the comprehensive nature of climate change, to better elucidate the specific impact of one or more variables, we conducted correlation analyses between the projected areas of different suitability levels (obtained through prediction) and the climate variable data for each sample point. The results are presented in Fig. 5. The correlations varied among unsuitable, low-, medium-, high-, and total suitable areas with different climatic variables. Unsuitable areas demonstrated sensitivity to temperature variation, consistently showing negative correlations. Bio1, Bio6, Bio7, Bio8, Bio9, and Bio11 exhibited significant correlations with unsuitable areas, as did Bio3, Bio5, and Bio10. Among precipitation variables, only Bio18 showed a significant negative correlation with unsuitable areas. Low-suitable areas positively correlated with most variables. Temperature variables Bio1, Bio2, Bio5, Bio6, Bio8, Bio9, Bio10, and Bio11 showed highly significant correlations with low-suitable areas, while Bio4 demonstrated a significant correlation. Precipitation variable Bio18 was significantly correlated, and Bio12, Bio13, Bio16, and Bio18 showed significant correlations with the low-suitable area. The total suitable area also exhibited sensitivity to temperature variations, with Bio1, Bio6, Bio8, and Bio9 showing highly significant correlations, and Bio3, Bio5, and Bio10 being significantly correlated. Notably, medium- and high-suitable areas did not demonstrate significant correlations with any single climatic variable (Fig. 5a). Grey relational analysis (GRA) was employed to analyze the correlation between climatic variables and suitable areas across all periods to determine which climate variable exerts the greatest impact on the suitable area. The influences of 18 climate variables on the un-, low-, medium-, high- and total suitable areas were similar (0.851–0.984). Only the weight of Bio6 was markedly lower than those of other climate variables, with values of 0.677 for unsuitable, 0.680 for low-suitable, 0.675 for medium-suitable, 0.687 for high-suitable, and 0.679 for total suitable areas. Bio6 represents the minimum temperature of coldest month; this value suggests that the effect of low temperature on the adaptability of P. pseudocerasus was not as pronounced compared to that of other climate variables (Fig. 5b). The mean results of the weight for the un-, low-, medium-, high- and total suitable areas also indicated that Bio9 (0.886) and Bio11 (0.884) were slightly lower than the other variables. The remaining variables were similar (all > 0.940), except for Bio6 (0.680) (Fig. 5c). To investigate whether climate variables in different periods and sample points affect the results of correlation analysis, we analyzed the coefficient of variation (CV) of all climate variables (Fig. 5d). The CV of Bio6 was the largest among all climatic variables (1.667), it indicates a large north–south span, and also proved that low temperature was not the critical variable in adaptability of P. pseudocerasus.

Fig. 5

Climate dependence of P. pseudocerasus suitability. a Correlation analysis between 19 climatic variables and the five suitability area categories. b GRA of 19 climatic variables and the five suitability area categories. c Mean value of the weight for the five suitability area categories. d CV of all climatic variables. In a, ** denotes significance of the correlation at the 0.01 level (2-tailed). In b, US: Unsuitable area, LS: Low-suitable area, MS: Medium-suitable area, HS: High-suitable area, TS: Total suitable area. The green circles indicate precipitation variables, while green squares represent temperature variables

Ecological niche comparisons

The ecological niche represents a species' capacity for resource utilization and adaptability, while ecological overlap indicates the state of resource competition between species. The geographical distribution of P. pseudocerasus is partially determined by its fundamental niche, which represents the environmental conditions necessary for its survival. Figure 6 illustrates the ecological niche dynamics of P. pseudocerasus in China, comparing climate niche spaces between native and invasive areas. Under each pathway, SSP585-2090s exhibited a low degree of ecological niche overlap (Fig. 6a). The ecological niches of SSP370s in invasive areas expanded compared to those in native areas. Figure 6b presents an alternative verification method for ecological niches, utilizing five climate variables for PCA to demonstrate more intuitively the ecological characteristics, revealing that different climate variables influence the niche across various pathways and periods. Table 1 presents the niche width results, indicating an increase in niche width across all pathways and periods (inverse concentration and uncertainty). These findings suggest a gradual increase in ecological adaptation. Figures S2 and S3 depict the ecological niche for the 2050s and 2070s periods. Schoener's D value greater than 0.5 indicated higher ecological niche overlap between periods under different CO₂ emission pathways. Thus, from the perspective of ecological niche overlap, no significant changes were observed between periods. However, the degree of niche overlap decreased with increasing CO₂ concentration. At the highest CO₂ concentration, the niche overlap index was only 0.696 (Table S10). The null hypothesis of ecological niche similarity based on bioclimatic variables of native and invasive areas was not rejected. Figure 6c demonstrates that the null hypothesis of ecological niche similarity based on different CO₂ emission pathways was not rejected (P = 0.04762, P < 0.05). As a species native to China, the evolution of its ecological niche aligned with expectations. Both D (D > 0.7) and I (I > 0.9) indices suggested a high degree of ecological niche overlap between different pathways and years. The ecological niches of different pathways and years may be similar but not identical.

Fig. 6

Ecological niche comparisons of P. pseudocerasus under different SSPs pathways in 2090s. a Ecological niches of P. pseudocerasus under SSP126–SSP585. b PCA load diagram for principal climate variables of ecological niches. c Ecological niche similarity across different pathways in 2090s. In a, Red arrows denote Schoener's D, blue indicates ecological niche overlap, green represents unfilling, and red signifies expansion. In c, the red arrow indicates the centroid of the realized ecological niche for each species

Table 1 Niche width between different pathways and periods of P. pseudocerasus

Degree of climate anomaly and spatial interpolation inversion

The mean multivariate similarity across the four future climate change scenarios ranged from 6.57 to 17.24. The SSP126-2050s climate scenario exhibited the highest multiple similarity and the lowest degree of climate anomaly, while the SSP585-2090s scenario demonstrated the lowest multiple similarity value and the highest degree of climate change (Figs. S4 and S6). As CO₂ emission concentrations increased for the same time period, the mean multivariate similarity value decreased (Figs. S4-6 and Table S11). At the locality level, the range spanned from −0.25 to 54.23. The lowest value was observed in Taihe County, Fuyang City, Anhui Province (SSP585-2090s), while the highest was recorded in Qinzhou District, Tianshui City, Gansu Province (SSP126-2090s). These differences primarily stemmed from variations in climate, altitude, geographic location, and terrain. In the modern suitable habitats for P. pseudocerasus, the most dissimilar variables were the mean temperature of the wettest quarter, mean diurnal range, precipitation in the warmest quarter, precipitation seasonality, and precipitation in the driest month. Precipitation variables dominated the vast majority of the adaptive area, followed by temperature variables. Terrain variables accounted for a negligible portion of the adaptation area.

The spatial interpolation results for P. pseudocerasus in Shandong and Guizhou are presented in Fig. S7. The timescale data were transformed into a spatial scale to illustrate the spatial distribution of various climate variables. All climate variables demonstrated similar current and future distribution patterns, with a notable spatial migration strongly correlated with climate change (Tables S3 and S5). We computed the error value and mean square error value for the current period and the 2090s (Fig. S7, Tables S12, 13). The mean square error range for most variables fell between −10 and 15 (Fig. S8). Notably, under the SSP585 pathway, Bio18 exhibited the largest error value among the five climate variables in Shandong Province for the 2090s. This discrepancy may be attributed to the uneven distribution of the distribution records used in the interpolation. In summary, Guizhou Province displayed smaller interpolation errors, while Shandong Province demonstrated greater irregularity, potentially due to its more complex geographical environment.

Discussion

Potential distribution area of P. pseudocerasus and reliability of MaxEnt prediction

This study predicted potentially suitable areas for P. pseudocerasus and analyzed the effects of 22 environmental variables on their growth and distribution. Yunnan, Sichuan, and Hunan provinces emerged as the most suitable distribution areas in China (Fig. 1a), aligning with records in Flora of China and other sources. Taiwan Province exhibited low- and medium-suitability areas among the current potentially suitable regions after simulation. Although reliable geographical information on P. pseudocerasus in Taiwan was not found in relevant literature, this finding corroborates the actual situation (Wu et al. 2003). The model accuracy test method employed in this study was the AUC method, with values ranging from 0.5 to 1. Values approaching 1 indicate higher accuracy. The AUC value of the test and training data was 0.927 and 0.928, respectively, suggesting high reliability of the selected geographic information (Phillips et al. 2006). As P. pseudocerasus is a large woody plant, it possesses a wider niche and higher degree of conservation compared to herbaceous and dwarf shrubs. Consequently, this study utilized the MaxEnt model to simulate the potential distribution of P. pseudocerasus. However, MaxEnt has limitations in considering biological interactions, tending towards a conservative niche. Overly complex models may present challenges in ecological interpretation. To mitigate this error, sample deviation processing and variable correlation analysis were employed to enhance the model's prediction accuracy. Nevertheless, the default parameters in MaxEnt were not adjusted according to species characteristics, potentially introducing errors and uncertainties that warrant further investigation to refine the model.

Model verification is a crucial final step in the modeling process, as accuracy directly determines data reliability (Deng et al. 2022). In this study, model accuracy was assessed using the ROC response curve. Actual climate data from the two provinces with the most extensive species distribution were utilized to compare similarity and accuracy with the model data. The precipitation data accuracy for Shandong and Sichuan provinces was low (< 50%) based on the confusion matrix (45.8% and 38.8%, respectively). These results may be attributed to the confusion matrix being a fuzzy-index evaluation method. Additionally, the classification standard of precipitation levels by the China Meteorological Administration, used as the label in this study, may contribute to this outcome (Deng et al. 2016). Precipitation is influenced by complex variables, and models capable of accurately predicting future precipitation have not yet been developed, despite precipitation being a critical environmental factor affecting species distribution and migration. Nevertheless, based on high-precision temperature data, linear regression analysis, and correlation analysis, we conclude that the model simulation results are more accurate.

Environmental variables influence the geographical distribution of P. pseudocerasus

Environmental conditions significantly influence plant growth and development, serving as primary determinants of species distribution (Thuiller et al. 2007). Following PCA and correlation analysis of 22 environmental variables, eight variables were identified as having the most substantial impact on the geographical distribution of P. pseudocerasus. These variables encompassed temperature, precipitation, and geographical factors. Previous MaxEnt-based studies on species distribution in China have demonstrated that these eight climate variables also exerted significant influence on the geographical distribution of other plant species, including Ulmus lamellosa, Scutellaria baicalensis, and Dalbergia cultrata Graham ex Benth (Yan et al. 2017; Liu et al. 2019; Xu et al. 2020 ; Zhao et al. 2018). Isothermality exhibited a positive correlation with mean diurnal temperature range (Bio2) and a negative correlation with the annual temperature range (Bio7). Optimal suitability was observed when the mean temperature of the wettest quarter (Bio8) ranged from 20–25℃, indicating that extreme temperatures, whether high or low, can diminish adaptability (Zhou et al. 2023). Excessively high temperatures adversely affect pistil and stamen development in P. pseudocerasus, potentially leading to abortion (Li et al. 2010a, b), while excessively low temperatures reduce cell viability, delay development, and may cause frostbite (Wang et al. 2019; Guy et al. 2008). Interestingly, the CV of minimum temperature of coldest month (Bio6) was the largest among all climate variables, it indicates that P. pseudocerasus was less sensitive to low temperature and more cold-tolerant. The phenological period of P. pseudocerasus precedes that of sweet cherry, with flowering typically occurring from February to April, new shoot growth peaking from February to May, and fruit ripening from April to May. However, due to significant geographical variations, the phenological period can vary considerably (Li et al. 2010a, b; Wang et al. 2022). Most plants undergo a warm-cold-warm stage during their growth and development periods. Consequently, temperature variations (mean diurnal range, daily temperature difference, and annual temperature difference) likely represent one of the most critical variables affecting plant growth and development (Fadón et al. 2019).

The contribution rates of precipitation in the wettest month (Bio18) and driest month (Bio14) were the most significant among all climate variables. These variables influence the geographical distributions of Alternanthera philoxeroides, Chromolaena odorata, and Lantana camara (Yan et al. 2020; Saranya et al. 2021). Adequate water facilitates plant growth and development during peak periods. However, P. pseudocerasus exhibits poor waterlogging tolerance. Excessive soil water content can compromise tree vigor and potentially lead to mortality. Furthermore, excessive precipitation during maturation directly results in fruit cracking (Cui et al. 2009). Fruit cracking in P. pseudocerasus has not received substantial attention, although it suggests that excessive rainfall reduces environmental adaptability (Suran et al. 2019). Analysis of precipitation in the potential growth areas revealed that despite differences in topography and climate types in the Sichuan Basin, Yunnan-Guizhou Plateau, North China Plain, northern Jiangsu Plain, and Taiwan Mountains, common characteristics were observed. Notably, more than 50% of the annual rainfall is concentrated in the summer, coinciding with fruit maturity (Chen et al. 2015, 2016) and generally aligning with the biological habits of P. pseudocerasus. Climate variability is evident, with precipitation seasonality (Bio15) having the third-highest water variable contribution rate, correlating strongly with the rainfall characteristics of the current potential distribution areas. The variance positively correlated with seasonal variations in precipitation fluctuations. Within a certain range, greater variance indicates better adaptability.

Climate variables significantly influence geographical distribution patterns by regulating the growth and development of species. Hydrothermal variables, in particular, play a crucial role in shaping ecological adaptation. These findings align with studies conducted on other species (Lu et al. 2018; Zhao et al. 2018). Climate change primarily disrupts the fit between species' adaptive traits and their environment. However, plants have evolved mechanisms to cope with climate-related disturbances (Aitken et al. 2008). In the case of P. pseudocerasus, further research is necessary to investigate the application of climate scenario simulations using various climate system models and to conduct comparative analyses of habitat differences in Chinese cherry populations.

Changes in distribution patterns and protection and management of P. pseudocerasus

The analysis of potential distribution areas in the current and future periods revealed that the present highly suitable areas for P. pseudocerasus primarily encompass the eastern Himalayas (southeast of the Tibet Autonomous Region), Hengduan Mountains (at the junction of eastern Sichuan, Yunnan, and Guizhou), Shandong, northern Jiangsu, and eastern Henan. The simulations for suitable areas in the 2050s, 2070s, and 2090s demonstrated minimal changes compared to the geographical distribution and current potential suitable areas. These areas remained predominantly concentrated around the Qinling Mountains-Huaihe River, with highly suitable regions primarily at the eastern and western extremities. The suitable regions encompass diverse terrains, including the Sichuan Basin, Yunnan-Guizhou Plateau, North China Plain, and Jianghuai Plain. These areas provide adequate atmospheric and thermal conditions, and their numerous mountains offer a relatively stable environment for plants during this period of climate fluctuation (Jordi et al. 2011).

Geographical and climatic variables, including landform, altitude, precipitation, and temperature, revealed the high niche diversity of P. pseudocerasus in China and its capacity to adapt to diverse environments with significant variations in climate, terrain, and altitude. The findings indicate the robust climate adaptability of P. pseudocerasus. Analysis of potential distribution areas across three different periods demonstrated that the total suitable area did not change substantially, possibly due to the absence of significant climate change. Most species exhibit niche conservatism and thus tend to inhabit areas with the most favorable climatic conditions. Consequently, a lesser degree of climate change would result in less pronounced alterations in species distribution (Parmesan et al. 1996). The potentially suitable habitat for this species has experienced considerable reduction in areas with minor topographic and geomorphic changes.

In recent years, global warming, driven by increased greenhouse gas emissions and a marked rise in extreme weather events, has significantly impacted the habitats of numerous organisms (Barnaud et al. 2010; Hirabayashi et al. 2022; Wu et al. 2022). Globally, species migration generally exhibits a trend towards higher altitudes and latitudes. However, the SDM results indicate that the potential future distribution area of P. pseudocerasus does not demonstrate a clear pattern of migration to higher altitudes. This observation may be attributed to either minimal climate change in its habitat or the species' long evolutionary history and successful adaptation to its current ecological niches (Shi et al. 2014). Consequently, future research should incorporate additional variables, such as soil characteristics, growth patterns, and anthropogenic influences, to enhance model accuracy and provide a more comprehensive projection of the future distribution of P. pseudocerasus.

To enhance the classification of climate variable changes across different pathways and periods, we conducted a cluster analysis of 151 distribution records based on five selected climate variables. The results are presented in Fig. S9. The Ward method was employed to categorize the distribution records into five distinct groups. The clustering outcomes for the current and 2090s periods demonstrated remarkable consistency in terms of group numbers, distances, and members. Notably, geographical proximity did not necessarily correlate with similar clustering results for all sampling points. Some geographically distant sampling points exhibited close clustering distances. While sampling points within the same province were partially clustered, others were dispersed, particularly in Sichuan and Shandong. This pattern suggests two key insights: firstly, it indicates high terrain and climate diversity in these regions; secondly, it points to significant climate change occurring in these areas.

P. pseudocerasus has an extensive cultivation history in China, yet only a limited number of wild resources have been commercially cultivated for an extended period. The slow promotion rate in China is attributed to the challenges associated with long-distance transportation after the color-changing stage. To develop commercial P. pseudocerasus cultivars, increased protection should be afforded to wild resources, and efforts should be made to artificially domesticate high-quality wild accessions. Elucidating the ecological characteristics of P. pseudocerasus will facilitate its introduction into suitable areas for adaptive cultivation, thereby promoting and preserving P. pseudocerasus resources in China.

Conclusion

The MaxEnt model provided geographical distribution information of P. pseudocerasus and corresponding environmental variable data. This study predicted that the primary potential distribution areas for cherry cultivation in China are Yunnan, Sichuan, and Hunan provinces. These regions are considered suitable for commercial cultivation. Additionally, most areas of Shandong Province and Dalian City were identified as highly suitable for P. pseudocerasus growth. Among the climate models examined, the SSP246 CO₂ emission concentration scenario yielded the largest total suitable area for P. pseudocerasus. While the 22 environmental variables selected in this study offer valuable insights, they may not comprehensively represent all climatic factors influencing the geographical distribution of P. pseudocerasus. Hydrothermal conditions, abiotic factors (such as light, air, and soil), species interactions, and migration ability significantly impact species distribution. Future research should incorporate a broader range of biological and abiotic factors during model construction to enhance prediction reliability. This study elucidates the response of P. pseudocerasus geographical distribution patterns to climatic variables, providing a scientific foundation for the conservation, development, and utilization of wild germplasm resources.

Materials and methods

Species geographical distribution data

The geographical distribution records of P. pseudocerasus were compiled and analyzed based on the checklist of SDMs. Initially, the geographical distribution data of P. pseudocerasus were sourced from multiple databases: the Chinese Virtual Herbarium (http://www.cvh.org.cn), National Specimen Information Infrastructure (NSCII, http://www.nsii.org.cn), Teaching Specimen Resource Sharing Platform (http://mnh.scu.edu.cn), Global Biodiversity Information Facility (https://www.gbif.org/), and Flora of China (http://www.iplant.cn/foc). Subsequently, duplicate distribution records and those lacking detailed geographical information were eliminated. The longitude and latitude of distribution records at different locations were then converted into actual projected linear distances. Lastly, distribution records with intervals < 5 km were removed. The resulting 151 distribution records of P. pseudocerasus are presented in Table S1.

Environmental variables

Raster files for 22 bioclimatic variables (19 climate and three digital elevations) with a resolution of 5 arc minutes were obtained from WorldClim (v2.1, https://worldclim.org/). Future climate and digital elevation data were acquired from the WorldClim database (http://www.worldclim.org//; accessed 2023). The Climate Model Intercomparison Project 6 (CMIP6), which offers significant advantages in climate simulations, was selected. Shared Socioeconomic Pathways (SSPs) were chosen for climate data during the 2050s, 2070s, and 2090s. Atmospheric CO₂ content is the primary factor influencing global temperature. CMIP6 was selected due to its high accuracy in predicting changes in CO₂ concentration. It is a universal climate system model applicable to most terrestrial plants. Human emissions of CO₂ and other greenhouse gases are the main drivers of climate change (Mitchell et al. 2000). Using the provincial administrative division vector map as the base map, data for the 22 variables for the three periods were cropped in ArcGIS to obtain the climate layer for the corresponding period. To determine the geographical distribution of suitable habitats for a species, it is crucial to identify the environmental variables affecting species growth. The environmental variables and their definitions are listed in Table 2.

Table 2 The 22 environmental variables of this study

Analysis of dominant environmental variables

The ArcGIS software's Extract Multi-values to Points tool was employed to extract environmental variable information from the 151 distribution records. Due to the correlation among environmental variables, it was necessary to select and remove certain variables to prevent overfitting of the species distribution model. PCA and Spearman rank correlation tests were conducted using SPSS v.24.0 software on the environmental variables. Variables with high variance contribution rates were identified, and correlation coefficients were compared to eliminate variables with high correlations (|R|> 0.7). In instances where two variables exhibited a high correlation, the variable with the higher contribution rate was retained.

Suitable habitat distribution model

Biomod2 (BIOdiversity MODelling2), utilizing distribution point data and environmental variables, predicted suitable habitat distribution. After evaluating all models' accuracy through AUC and TSS values, the optimal model was selected for further analyses. The contribution of each environmental variable was determined using PCA (Akomolafe et al. 2020). Response curves were established to obtain logistic relationships between distribution probability and environmental variables. The ROC curve was employed to assess the model's prediction accuracy, with the random test percentage set at 25, model replicates at 10, and remaining parameters at default values (Segurado et al. 2004). ArcGIS software was utilized for post-processing simulation results, including hierarchical displays, area statistics, and centroid calculations (Yan et al. 2021). Following predictions and analysis using the MaxEnt model, the potential distribution areas of P. pseudocerasus were categorized into four grades using the natural segmentation method. The suitable habitat class of P. pseudocerasus was determined based on the P value, where P < 0.2 indicated unsuitable areas, 0.2 ≤ P < 0.5 denoted low-suitable areas, 0.5 ≤ P < 0.7 represented medium-suitable areas, and P ≥ 0.7 signified high-suitable areas (Xu et al. 2019; Sharifian et al. 2021).

Model accuracy test

The model evaluation was conducted by inputting the acquisition point distribution data and environmental variable layers into MaxEnt. The accuracy of the suitable area distribution predicted by MaxEnt was assessed using the ROC curve, from which the AUC was derived. The AUC values range from 0.5 to 1, with values closer to 1 indicating higher accuracy (Swets et al. 1988). The interpretative standards for AUC values were as follows: 0.50–0.60 (fail), 0.61–0.70 (poor), 0.71–0.80 (fair), 0.81–0.90 (good), and 0.91–1.00 (excellent). The confusion matrix analysis was performed using MATLAB R2022b, employing the following calculation formulae:

$$\left\{\begin{array}{l}\mathrm{YValidatio}\;\mathrm n\;=\;\mathrm{data}.\;\mathrm{VarName}2;\\\mathrm{YValPred}=\mathrm{data}.\;\mathrm{VarName}2;\\\mathrm{Plabels}=\mathrm{categorical}(\mathrm{YValPred});\\\mathrm{Alabels}=\mathrm{categorical}(\mathrm{YValidat}\;\mathrm{ion});\\\mathrm{plotconfus}\;\mathrm{ion}(\mathrm{alabel}\;\mathrm s,\;\mathrm{plabels});\end{array}\right.$$

The extreme values derived from temperature and precipitation data were categorized into five classes for subsequent analysis. Historical climate data were obtained from various weather stations, while model climate data were acquired from the National Tibetan Plateau / Third Pole Environment Data Center (https://data.tpdc.ac.cn/home).

Measurement of ecological niche and spatial interpolation

To enhance understanding of the influence of CO₂ concentration and temporal factors on environmental variables, we integrated the MaxEnt modeling results to obtain spatial interpolation data. Thin plate smoothing spline was employed to invert spatiotemporal changes in environmental variables. The ecospat package in R Studio 4.1.3, utilizing PCA-env and centroid shift, overlap, unfilling, and expansion schemes, was implemented to analyze ecological niches (Broennimann et al. 2012; Guisan et al. 2014; Di et al. 2017). Relevant bioclimatic variables were used to perform PCA-env. Subsequently, the ecospat package was utilized to conduct climate niche similarity tests, performed bidirectionally with 1000 repetitions. Enmtools1.3 was employed to determine niche differences between various CO₂ concentration pathways and periods, including niche overlap and width. In the ecological niche overlap test, Schoener's D and Hellinger's distance (I) were utilized to quantify similarities between the ecological niches of the three periods and to compare the degree of ecological niche overlap in their native areas (Xian et al. 2023). The values of D and I ranged from 0 (no overlap) to 1 (> 0.6 indicating significant overlap).

Analyses of multivariate environmental similarity surface and most dissimilar variable

Utilizing the environmental variables in the current potential distribution area (101.5395°–121.8275°E, 23.4476°–37.3556°N) of P. pseudocerasus as the reference layer, a multivariate environmental similarity surface was employed to assess the degree of climatic anomaly in the current suitable habitats under past and future climate change scenarios. The most dissimilar variable was analyzed to identify the key factors driving changes in the potential geographical distribution. The multivariate environmental similarity surface calculated the similarity (S) between environmental variables under future climate conditions and contemporary environmental variable point sets to determine the extent of environmental change in the distribution area. The environmental variable with the lowest S-value (highest degree of abnormality) at a point was considered the least similar variable at that location. Environmental variables exhibiting the highest levels of anomalies were likely to be the primary factors causing shifts in the distribution of suitable habitats. This analysis was conducted using the ‘density.tools.Novel' command in the ‘maxent.jar' file.

Migration of the center point of suitable area

The SDMTool package in R software was employed to analyze the trends in suitable habitat regions and the geometric center positions of these areas for the current period, as well as the 2050s, 2070s, and 2090s. The suitable habitat was considered as a single entity and reduced to a vector particle, utilizing the centroid position changes to reflect the size and direction of the species' suitable habitat. Subsequently, the centroid was tracked using various SDMs to examine its position across different time periods and under varying climatic conditions. This approach facilitated the evaluation of the suitable zone's migration distance in terms of latitude and longitude coordinates.

Statistical analyses

Statistical analysis was conducted using Microsoft Excel 2020. SPSS v.24.0 was employed for the PCA and best-fit analysis, while Pearson's method was utilized for linear and partial correlation analyses. PCA was performed to explore correlations among variables and 151 distribution records, determine maximum or minimum correlation directions, and achieve data compression to enhance result accuracy (Abdi et al. 2010). Data fitting was conducted using a linear regression equation based on the discrete form of the data (Zhao et al. 2022). Additionally, GRA was employed to analyze 19 climate variables to determine which had a greater impact on suitable areas. To more intuitively understand the relationship between correlated climate variables and suitability, we analyzed the variability of climate variables across different periods to demonstrate the impact of fluctuations on suitability. Microsoft Excel 2010 was used for CV analysis. Ward's method, Euclidean distance, and sum of distances cluster analyses were performed using Origin Pro 2022 software (Xie et al. 2020). Analysis of variance and Duncan's multiple range test were conducted to determine significant (P = 0.05) differences between groups. The ggcor package in R software was employed to analyze the correlations of variables and projected areas. Trends and figures within the data were graphically depicted using OriginPro 2022 and RStudio. All data used for correlation analysis and PCA underwent data standardization.