Global potential distribution of Paeonia lactiflora and its climate-driven shifts: insights from an enhanced MaxEnt model integrating soil and solar radiation variables

TL;DR

This study predicts the future global habitat of Paeonia lactiflora using a detailed climate model, showing expansion and shifts due to climate change.

Contribution

The study introduces an enhanced MaxEnt model integrating soil and solar radiation variables to predict Paeonia lactiflora's distribution.

Findings

The model shows excellent performance with AUC = 0.945 and TSS = 0.762.

Future climate scenarios project a 25–45% expansion of suitable habitat with a northeastward shift.

High-emission scenarios cause contraction in southern core habitats and require conservation planning.

Abstract

Paeonia lactiflora Pall. is a globally important medicinal perennial whose habitat suitability remains poorly known beyond China. Using an enhanced MaxEnt model integrating 45 climatic, soil, and solar radiation variables, we predicted its current and future global distribution based on 833 spatially thinned occurrence records and 12 low-collinearity predictors. The model performed excellently (test AUC = 0.945 ± 0.001; TSS = 0.762 ± 0.018). Precipitation of the warmest quarter (bio18), mean temperature of the coldest quarter (bio11), temperature seasonality (bio4), and November solar radiation (srad11) were the dominant drivers. Currently, total suitable habitat is centered in East Asia, central Europe, and northeastern/midwestern USA. All future scenarios (SSP2-4.5 and SSP5-8.5, 2041–2060 and 2061–2080) project about 25–45% expansion of total suitable area, accompanied by a consistent northeastward centroid shift of highly suitable habitat (up to ~1,234 km under SSP5-8.5 2061–2080). Late-century high-emission conditions cause localized contraction of core habitat in southern margins. P. lactiflora is likely to benefit from moderate warming, but high-emission pathways will drive major reorganization and degradation of traditional production areas, necessitating strengthened conservation in current strongholds and proactive planning in emerging northern regions.

Figures

Distribution of P. lactiflora, with each blue dot indicating a recorded occurrence point.

Heatmap of Pearson correlation coefficients among the initial 45 environmental variables.

Optimized methodological workflow of the MaxEnt species distribution modeling process for P. lactiflora.

The ROC curves of the MaxEnt model for P. lactiflora.

Jackknife test of environmental variables for P. lactiflora.

The permutation importance and percent contribution.

Response curves for bio18, bio11, srad11 and bio4.

Current potential distribution of P. lactiflora.

Future climate–driven potential distribution of P. lactiflora.

Changes in the area of different suitability classes for P. lactiflora under current and future climate scenarios.

Spatial patterns of suitable habitat gain and loss for P. lactiflora under four future climate scenarios.

Shift map of the center of gravity in the suitable area of P. lactiflora under future climate scenario.

Tables

Table 1. ENMeval results comparing default and selected parameter settings.

Type	FC	RM	AUCtrain	AUCdiff.avg	OR10.avg	ORMTP.avg	ΔAICc
Selected	LQH	1.5	0.952	0.047	0.253	0.001	0.000
Default	LQHP	1	0.949	0.068	0.412	0.004	46.37

Table 2. Environmental factors and their contribution rates after screening.

Variable	Description	Unit
bio_4	Temperature Seasonality (Standard Deviation * 100)	°C
bio_11	Mean Temperature of Coldest Quarter	°C
bio_15	Precipitation Seasonality (Coefficient of Variation)	mm
bio_18	Precipitation of Warmest Quarter	\
bio_19	Precipitation of Coldest Quarter	\
Elevation	altitude	m
Slope gradient	Gradient of the terrain	°
drainage	Drainage class	\
t_silt	Topsoil Silt Fraction	% wt.
t_clay	Topsoil Clay Fraction	% wt.
srad07	Solar radiation in July	KJ/m2/day
srad11	Solar radiation in November	KJ/m2/day

Table 3. AUC value and TSS value of each component module in the study.

Type	current	SSP245-2050	SSP245-2070	SSP585-2050	SSP585-2070
AUC	0.946	0.945	0.946	0.945	0.944
TSS	0.782	0.735	0.769	0.764	0.762

Table 4. Area changes of suitable habitats under four future scenarios relative to current conditions.

Status	SSP2-4.5-2050	SSP2-4.5-2070	SSP5-8.5-2050	SSP5-8.5-2070
Stable	1104.11	1113.61	1118.58	1096.60
Gain	234.79	422.02	399.68	666.83
Loss	163.02	153.40	148.49	170.20

Full text

Introduction

1

Paeonia lactiflora Pall., a perennial herbaceous species in the family Paeoniaceae, has been used medicinally for over two millennia (Shagjjava et al., 2016). Its earliest recorded medicinal use appears in the “Prescriptions for Fifty-two Diseases,” a medical text unearthed from the Mawangdui Han Tomb in Changsha, China. For more than 2,000 years, P. lactiflora and its processed root products—Paeoniae Radix Alba and Paeoniae Radix Rubra—have played a central role in traditional Chinese medicine. According to the Chinese Pharmacopoeia, the primary authentic production regions (dao-di areas) are located in northeastern China, northern China, Shaanxi Province, and northern Gansu Province, although the species is now widely cultivated throughout the country.

P. lactiflora exhibits relatively strong environmental adaptability and typically occurs on grassy slopes, forest margins, and in thickets within temperate regions. Its growth and development are strongly influenced by environmental factors, particularly light availability, soil moisture, and temperature (Zhang et al., 2019). Suitable ecological conditions are essential for the accumulation of its bioactive compounds. Chemical studies have shown that P. lactiflora is rich in monoterpene glycosides (primarily paeoniflorin), tannins, flavonoids, polysaccharides, and other active constituents, which underpin its significant medicinal and economic value (He and Dai, 2011).

Despite its long history of clinical use and well-documented efficacy, the environmental factors driving the formation of its authentic production regions (dao-di areas) and the extent of genetic diversity in germplasm resources remain poorly understood (Huang, 2011). Moreover, wild populations of P. lactiflora are increasingly threatened by habitat loss and overharvesting, posing serious challenges to germplasm conservation and sustainable utilization (Yadav et al., 2024).

Species distribution models (SDMs) have become indispensable tools for assessing species’ geographic patterns and their responses to environmental change, with wide applications in conservation biology, medicinal plant resource management, and biogeography. Among the various SDM approaches, MaxEnt (maximum entropy) stands out for its robust performance when using presence-only data and relatively small sample sizes, making it particularly well-suited for studying threatened or poorly documented medicinal plants (Phillips et al., 2004). The algorithm is grounded in the principle of maximum entropy, which provides a statistically rigorous framework for generating predictive distributions from incomplete information.

Proposed by Jaynes (1957), the principle of maximum entropy states that, among all probability distributions consistent with the available constraints, the one with the highest entropy is the least biased and most appropriate (Shore and Johnson, 1980). In the context of species distribution modelling, entropy reflects predictive uncertainty: higher entropy corresponds to a more uniform (less environmentally constrained) predicted distribution, whereas lower entropy indicates stronger environmental filtering (Banavar et al., 2010).

Unlike many traditional modelling approaches that rely on restrictive assumptions about species–environment relationships (e.g., linearity or predefined response curves), MaxEnt imposes minimal a priori constraints. It uses only empirically derived relationships between known occurrence records and environmental predictors, avoiding the extrapolation of untested ecological assumptions. This assumption-lean, data-driven framework makes MaxEnt especially valuable for modelling poorly surveyed or narrowly distributed species, including many medicinal plants with limited occurrence data (Lissovsky and Dudov, 2021).

In practice, MaxEnt constructs a probability distribution over the study area that maximizes entropy (i.e., is as uniform as possible) while remaining consistent with the empirical constraints derived from known occurrences. For a given plant species, these constraints typically include the observed ranges and variances of environmental variables (e.g., temperature and precipitation) at presence localities (Elith et al., 2011). The resulting distribution is the least committed to unobservable patterns, effectively balancing ecological signal against sampling noise and bias.

A key strength of MaxEnt is its ability to capture nonlinear responses and interactions among predictors without assuming specific functional forms (e.g., quadratic or threshold responses). Because it relies solely on empirically supported constraints, the model naturally accommodates complex, non-additive relationships that are common in real species–environment associations (Halvorsen et al., 2016).

ArcGIS provides robust tools for spatial data processing, standardization, and visualization, making it widely used in species distribution modelling (Wong and Lee, 2005; Scott and Janikas, 2009). In this study, it served three primary functions: (1) pre-processing and formatting of occurrence records and environmental layers to meet MaxEnt input requirements; (2) derivation of additional topographic and bioclimatic variables; and (3) post-modelling visualization, including production of habitat suitability maps, calculation of suitable area, and centroid analysis of highly suitable regions (Yan et al., 2020; Zeng et al., 2021; Xu et al., 2024). This integrated GIS–MaxEnt workflow substantially streamlines data preparation, model implementation, and result interpretation.

Climate change is profoundly affecting the growth, phenology, and geographic distribution of plant species worldwide (Li et al., 2022). Although previous studies have used MaxEnt and Biomod2 to model the potential distribution of P. lactiflora within China (Bi et al., 2022; Wang et al., 2024), the species also maintains naturalized or cultivated populations in parts of East Asia, Europe, and North America (Shagjjava et al., 2016). To date, no study has assessed its global habitat suitability or projected climate-driven range shifts at a worldwide scale.

This study addresses these gaps by: (1) compiling a comprehensive global occurrence dataset; (2) incorporating a broader suite of environmental predictors, including soil properties and solar radiation in addition to climatic variables; and (3) evaluating both current suitability and future range dynamics under multiple shared socioeconomic pathways (SSPs) and time periods. These improvements are expected to yield more robust and biologically realistic predictions than regional models that rely solely on climatic data.

Materials and methods

2

Acquisition of global P. lactiflora geographic distribution data

2.1

Global occurrence records of P. lactiflora were downloaded from the Global Biodiversity Information Facility (GBIF.org (17 July 2025) GBIF Occurrence Download https://doi.org/10.15468/dl.vpzb7g). We retained only records with valid geographic coordinates collected between 1970 and 2025, yielding an initial set of 1,434 presence points. Records lacking coordinates, duplicate entries, or obvious georeferencing errors (e.g., located in the ocean or at coordinate origin 0,0) were removed during initial cleaning (Radosavljevic and Anderson, 2014).

To reduce spatial autocorrelation and sampling bias caused by clustered collections, the dataset was spatially thinned using the “Spatially Rarefy Occurrence Data” tool in SDMtoolbox 2.0 for ArcGIS (Isaac et al., 2020). A thinning distance of 20 km was applied, resulting in a final set of 833 spatially independent occurrence records (Figure 1). This thinning distance was chosen to approximate the spatial resolution of the environmental layers used (ca. 5 arc-min, ≈10 km at the equator) while retaining sufficient records for robust modelling.

Distribution of P. lactiflora, with each blue dot indicating a recorded occurrence point.

Processing of environmental variable data

2.2

The distribution of herbaceous perennials in the Paeoniaceae is strongly governed by moisture availability, temperature extremes, solar radiation, and soil physical and chemical properties (Shagjjava et al., 2016). To capture these drivers at a global scale, we compiled an initial set of 45 environmental layers representing climatic, edaphic, and radiation conditions.

Current (1970–2000) bioclimatic variables (19 variables), monthly precipitation, solar radiation (srad), and shortwave radiation downward flux (swd) were obtained from WorldClim version 2.1 at 5 arc-minute resolution (~10 km at the equator) (Fick and Hijmans, 2017). Elevation was derived from the GMTED2010 dataset and resampled to the same grid. Soil data were extracted from the Harmonized World Soil Database v1.2 (HWSD). Soil layers were held constant across future periods, a standard assumption in global-scale SDM studies (Slessarev et al., 2019). This practice is justified because soil properties (e.g., texture, organic carbon, pH) form and change over centennial to millennial timescales through processes such as weathering, organic matter accumulation, and erosion, which are far slower than the decadal-to-centennial pace of projected climate change under CMIP6 scenarios (Meyer et al., 2025). High-resolution global soil data for future conditions are currently unavailable, and dynamic soil models are computationally intensive and not yet sufficiently validated for broad-scale projections (Ni and Vellend, 2024).

Nevertheless, this assumption introduces limitations. Soil properties may shift indirectly under prolonged climate change through altered vegetation dynamics, erosion rates, or microbial activity, potentially modifying habitat suitability in ways not captured here. In particular, in arid or high-latitude regions where soil formation is slow but climate-driven degradation (e.g., desertification or permafrost thaw) is rapid, future suitability could be overestimated or underestimated (Huang et al., 2017). Future refinements could incorporate dynamic soil–climate feedbacks or scenario-based soil projections when such data become available. To maintain consistency, future soil layers were assumed to remain unchanged, a standard practice in global-scale SDM projections given the slow rate of soil property change relative to climate. From the HWSD, we retained ten topsoil (0–30 cm) attributes known to influence root development and nutrient uptake in medicinal Peony species: drainage class, USDA texture class, sand fraction, clay fraction, pH (water), organic carbon, CaCO_3_ content, electrical conductivity, bulk density, and cation exchange capacity of the clay fraction.

Future climate layers (2041–2060 and 2061–2080) were sourced from the BCC-CSM2-MR model under two Shared Socioeconomic Pathways: SSP2-4.5 (intermediate emissions) and SSP5-8.5 (fossil-fueled development, high emissions) within CMIP6. These periods are hereafter referred to as the 2050s and 2070s, respectively. Future climate layers (2041–2060 and 2061–2080) were sourced from the BCC-CSM2-MR model under SSP2-4.5 and SSP5-8.5 within CMIP6 (Wu et al., 2019). BCC-CSM2-MR was selected for its strong performance in simulating East Asian monsoon precipitation and temperature seasonality—key drivers of P. lactiflora distribution—particularly in the native range and analogous temperate zones (Xin et al., 2020). While multi-model ensembles are increasingly recommended to reduce GCM-specific uncertainty, the use of a single GCM is common in regional-to-global medicinal plant SDMs when computational constraints or regional performance validation is prioritized. Therefore, prioritizing this regionally validated, high-performing single GCM provides more reliable projections at the regional scale under computational constraints (Cianfrani et al., 2010; Araújo et al., 2019).

Nevertheless, this choice introduces potential model-specific bias, particularly in projections of precipitation extremes or temperature variability. Future refinements could incorporate a multi-GCM ensemble (e.g., 3–5 models with high skill in temperate Asia and Europe) to quantify uncertainty and increase robustness of range-shift predictions.

All layers were clipped to a global extent (−180° to 180° longitude, −60° to 90° latitude), resampled to a common 5 arc-minute grid using bilinear interpolation, and converted to ASCII format using SDMtoolbox 2.0 for ArcGIS. Masking was applied to exclude permanent ice and water bodies.

Establishment of the MaxEnt model

2.3

To identify the key environmental factors shaping the global distribution of P. lactiflora, the MaxEnt algorithm (version 3.4.4) was employed (Phillips et al., 2025). To avoid overfitting and ensure ecologically realistic and statistically robust predictions, we first optimized the regularization multiplier (RM) and feature classes (FC) using the ENMeval R package (v2.0.0) (Muscarella et al., 2014; Bao et al., 2022). The modelling protocol was designed to maximize robustness and minimize overfitting. All 833 spatially thinned occurrence records and the environmental layers were imported into MaxEnt. ENMeval evaluated 48 parameter combinations (RM from 0.5 to 4.0 in 0.5 increments; FC combinations: linear, quadratic, hinge, product, threshold, and all subsets) using 10-fold cross-validation and the Akaike Information Criterion corrected for small sample size (AICc). The combination with the lowest AICc (RM = 1.5, FC = linear + quadratic + hinge) was selected as optimal. The optimal settings (FC = LQH, RM = 1.5) outperformed the default MaxEnt settings (FC = LQHP, RM = 1) in validation AUC, overfitting metrics, and AICc (see Table 1 for detailed comparison). The model was configured with 10-fold cross-validation, whereby the data were randomly partitioned into a 75% training subset and a 25% test subset in each fold (Wei-Yao et al., 2019). This approach provides an objective assessment of predictive performance on independent data and reduces the risk of inflated evaluation metrics.

To further enhance reliability, ten replicate runs were performed with different random seeds. A preliminary model using the full set of 45 environmental variables was first executed, and jackknife analysis was applied to estimate the individual contribution of each predictor (Wolter, 2007). Because strong multicollinearity was detected among the initial 45 variables, collinearity was subsequently diagnosed using ENMTools (Warren et al., 2021). Pairs of variables with a Pearson correlation coefficient |r| ≥ 0.80 were identified, and in each pair, the variable with the lower permutation importance in the preliminary model was removed. The complete Pearson correlation matrix heatmap for the initial 45 environmental predictor variables is shown in Figure 2. The color scale ranges from blue (−1, strong negative correlation) to red (+1, strong positive correlation), with white indicating no correlation. Absolute values ≥ 0.80 are considered highly collinear (threshold for removal) and are visually prominent in red/blue. Variables are labeled along the axes. This matrix was calculated across the entire global study area using ENMTools. This iterative screening process ultimately retained 12 non-redundant, biologically meaningful predictors.

Heatmap of Pearson correlation coefficients among the initial 45 environmental variables.

The final MaxEnt model was then rerun using only these 12 selected variables with the optimized settings (RM = 1.5, FC = linear + quadratic + hinge) while maintaining the same 10-fold cross-validation settings, auto-feature selection, regularization multiplier (β=1.5), and 10,000 background points. The definitive percent contribution and permutation importance of each retained predictor are reported in Table 2.

Data processing of prediction results and visual analysis

2.4

The MaxEnt logistic output (suitability values ranging from 0 to 1) was reclassified using a consistent threshold for binary and quantitative analyses. For presence/absence mapping, range-change detection, area statistics of suitable habitat, and centroid shift calculations, we applied the threshold that maximized training sensitivity plus specificity (MaxSSS; mean = 0.312 ± 0.017 across 10-fold cross-validation). This threshold is widely recommended for presence-only data because it balances omission and commission errors without requiring true absence information (Liu et al., 2005; Jiménez-Valverde and Lobo, 2007). Pixels exceeding 0.312 were classified as suitable.

To visualize relative suitability gradients and facilitate ecological interpretation, we additionally applied the Jenks natural breaks algorithm in ArcGIS to divide the continuous suitability values into four classes: unsuitable, marginally suitable, moderately suitable, and highly suitable (Chen et al., 2013). Current and future suitability layers were overlaid using the raster overlay tool to quantify spatial patterns of habitat expansion, contraction, and stability under each emission scenario–time-period combination. Shifts in the centroid of suitable habitat (>0.312) between present and future projections were calculated with SDMtoolbox 2.0.

To enhance clarity and reproducibility, the overall methodological workflow is presented in Figure 3. This diagram details the sequential process from GBIF occurrence data collection and spatial thinning, environmental variable preparation and collinearity filtering, parameter tuning using ENMeval, final MaxEnt modeling and evaluation, suitability reclassification (using MaxSSS for binary analysis and Jenks for gradient visualization), to post-processing analyses including area quantification, centroid shift calculation, and projections under SSP scenarios.

Optimized methodological workflow of the MaxEnt species distribution modeling process for P. lactiflora.

Results

3

Accuracy of MaxEnt model analysis

3.1

Model performance was evaluated using two widely accepted metrics: area under the receiver operating characteristic curve (AUC) and the True Skill Statistic (TSS) (Allouche et al., 2006; Chen et al., 2013; Carrington et al., 2023; Yoon and Lee, 2023).

The AUC measures the model’s ability to discriminate presence records from background points, independent of any specific threshold. Values >0.8 indicate good performance, whereas values >0.9 are considered excellent. In the present study, the mean test AUC across the 10 cross-validation folds was 0.945 ± 0.001 (SD), demonstrating outstanding discriminatory power (Figure 4; Table 3).

The ROC curves of the MaxEnt model for P. lactiflora.

The TSS complements AUC by incorporating a specific threshold, thereby directly assessing classification accuracy while remaining independent of prevalence. It ranges from −1 to +1, with values >0.7 generally regarded as evidence of high predictive reliability suitable for conservation and management applications (Bradie and Leung, 2017). The mean TSS obtained here was 0.762 ± 0.018, calculated in R v4.3.2 using the package ‘dismo’.

Collectively, the high AUC and TSS values confirm that the final MaxEnt model is both highly accurate and robust, making it well-suited for mapping the current and future global distribution of suitable habitat for P. lactiflora.

Assessment of primary environmental determinants governing the geographic distribution of P. lactiflora

3.2

The relative influence of each environmental predictor was assessed using MaxEnt’s built-in jackknife procedure on regularized training gain and test AUC (McIntosh, 2016). The four dominant variables collectively accounted for the majority of model explanatory power (Figure 5). Precipitation of the warmest quarter (bio18), mean temperature of the coldest quarter (bio11), solar radiation in November (srad11), and temperature seasonality (bio4) emerged as the most important drivers of global habitat suitability for P. lactiflora (Figure 6).

Jackknife test of environmental variables for P. lactiflora.

The permutation importance and percent contribution.

Percent contribution reflects the increase in regularized gain attributed to each variable during model fitting, whereas permutation importance measures the drop in test AUC when the values of that variable are randomly shuffled on withheld data; the latter is generally regarded as a more reliable indicator of true predictive relevance. In the final model, bio18 exhibited the highest permutation importance, followed by bio11, bio4, and srad11 (Figure 6, Table 2).

Marginal response curves for these four key predictors are shown in Figure 7. These curves illustrate how predicted suitability changes as each environmental variable is varied while all others are held at their average sample value, thereby revealing the realized niche boundaries of the species along each gradient.

Response curves for bio18, bio11, srad11 and bio4.

The response curves reflect the physiological tolerance and ecological niche preferences of a species. According to Shelford’s Law of Tolerance, each species possesses an ecological amplitude—that is, a tolerance range for every environmental factor, including a minimum threshold, an optimal range, and a maximum limit (Lynch and Gabriel, 1987). The response curve serves as a visual representation of this principle within the framework of mathematical modeling.

Precipitation of the warmest quarter (bio18), which represents water availability during the main growing season, was the most influential predictor. Suitability increased sharply when bio18 exceeded approximately 200 mm and reached its highest values above 280–350 mm, reflecting the species’ strong association with summer-monsoon climates and its sensitivity to growing-season drought.

Mean temperature of the coldest quarter (bio11) exerted a clear constraining effect on the northern and altitudinal range limits. Predicted suitability was highest when bio11 fell between −15 °C and 7 °C, with a pronounced optimum around −10 °C to 2 °C. Values below −15 °C or consistently above 7 °C resulted in a rapid decline in suitability, confirming the critical role of moderate winter cold for dormancy and spring regrowth.

Temperature seasonality (bio4) further restricted the species to regions with moderate annual thermal variation. Suitability declined markedly when bio4 exceeded approximately 850 (standard deviation × 100), effectively excluding continental interiors characterized by extreme temperature fluctuations.

November solar radiation (srad11) also emerged as an important driver. Optimal conditions occurred between roughly 4,000 and 12,500 kJ m⁻² day⁻¹, with suitability decreasing at both lower values (typical of higher latitudes in late autumn) and higher values (low-latitude regions).

Secondary contributions were made by elevation, soil drainage class, and precipitation of the coldest quarter, which primarily refined local-scale habitat quality within the broader climatic envelope defined by the four dominant variables.

Potential geographical distribution and habitat evaluation

3.3

Spatial pattern of the current potential distribution of P. lactiflora

3.3.1

Under present-day climate conditions (1970–2000), the MaxEnt model predicts the potential suitable distribution areas of peony globally, as shown in Figure 8. This extensive suitable range spans three continents and reflects both the species’ native East Asian distribution and the climatic analogues that support its successful cultivation and naturalization elsewhere.

Current potential distribution of P. lactiflora.

For descriptive purposes, suitability classes were derived using the Jenks natural breaks classification. Highly suitable habitat (logistic value >0.44) covers 5.53 million km² and forms three distinct core areas: the traditional Chinese cultivation heartland and native range in eastern China (especially Shandong, Henan, Jiangsu, Anhui, and surrounding provinces), together with the Korean Peninsula and central Honshu, Japan; a broad belt across central Europe, centered on Germany, Poland, the Czech Republic, Austria, and adjacent lowlands; and a smaller but well-defined region in the northeastern and midwestern United States, where the species has long been cultivated and occasionally naturalized.

Moderately suitable habitat (0.22-0.44), totaling 7.06 million km², surrounds these high-suitability cores and expands continuously into neighboring temperate zones. In China, it extends northward into Northeast China, westward into parts of North and Southwest China, and southward along lower-elevation corridors. In Europe, it includes most of France, the United Kingdom, Belgium, the Netherlands, and southern Scandinavia. In North America, the band covers much of the Upper Midwest and parts of the Mid-Atlantic states.

Marginally suitable habitat (0.07-0.22) accounts for the remaining 13.72 million km² and occupies transitional climatic zones at the outer edges of the species’ realized niche. These areas include the semi-arid northwest of China and the fringes of the Qinghai–Tibet Plateau, the Mediterranean littoral (northern Italy, coastal Balkans, Greece, and western Turkey), as well as the northern Pacific coast of the United States and interior valleys of the western Cordillera. Although these regions currently support only limited cultivation or occasional escape, they share climatic features (particularly adequate warm-season precipitation and moderate winter cold) that place them within the broader ecological amplitude of the species.

Visualization analysis of potential habitat distribution prediction for P. lactiflora under future climate conditions

3.3.2

Under the two Shared Socioeconomic Pathways examined (SSP2-4.5 and SSP5-8.5) and the two future time periods (2041–2060 and 2061–2080), the MaxEnt projections consistently indicate a substantial northward and poleward expansion of suitable habitat for P. lactiflora (Figure 9).

Future climate–driven potential distribution of P. lactiflora.

For visualization of future habitat change trends, the Jenks natural breaks classification was applied. Compared with the current total suitable area of 26.31 million km², all future scenarios exhibit marked increases (Figure 10). In the 2041–2060 period, the total suitable area rises to 32.97 million km² under SSP2-4.5 and 33.74 million km² under SSP5-8.5. By 2061–2080, the extent continues to grow, reaching 33.74 million km² under SSP2-4.5 and 38.27 million km² under SSP5-8.5 — representing a maximum increase of nearly 45% relative to present-day conditions.

Changes in the area of different suitability classes for P. lactiflora under current and future climate scenarios.

Highly suitable habitat shows more nuanced dynamics. It expands modestly in the mid-century (2041–2060) under both scenarios, but under the high-emission SSP5-8.5 pathway during 2061–2080 it contracts slightly by approximately 0.28 million km² compared with the 2041–2060 period, even as total suitable area continues to grow. This suggests that extremely high temperatures and altered seasonality under SSP5-8.5 begin to erode core habitat quality in some traditionally favorable regions by the late 21st century.

The newly suitable areas are predominantly located at higher latitudes: western Russia, southern Scandinavia, southeastern Canada, and interior portions of the northeastern United States become incorporated into the moderate- and high-suitability zones. Meanwhile, the original core regions in eastern China, the Korean Peninsula, Japan, and central Europe largely retain their high suitability under SSP2-4.5, whereas under late-century SSP5-8.5 some localized degradation is projected in the southernmost portions of these historic ranges.

Changes in suitable habitat areas under different climate conditions

3.3.3

Overlay analysis in ArcGIS was used to compare current suitable habitat with projections under the four future climate scenario–period combinations. To ensure accurate quantification of habitat gain and loss, suitability was binarized using the MaxSSS threshold in the analyses presented in this section. The resulting patterns of habitat gain, loss, and stability are displayed in Figure 11, while the corresponding area changes are quantified in Table 4.

Spatial patterns of suitable habitat gain and loss for P. lactiflora under four future climate scenarios.

All four scenarios project a net expansion of suitable habitat, but the magnitude and spatial expression differ markedly (Figure 9, Table 4). The greatest overall increase occurs under SSP5-8.5 for 2061–2080, where newly suitable area reaches its maximum and habitat loss, although present, remains considerably smaller than gains. This results in the largest net gain among all scenarios.

In contrast, the two SSP2-4.5 scenarios (both 2041–2060 and 2061–2080) and the mid-century SSP5-8.5 scenario exhibit more moderate and balanced dynamics, with smaller areas of both gain and loss. Habitat turnover is therefore lowest under SSP2-4.5, indicating relatively stable range configuration despite overall northward expansion.

The pronounced habitat loss under late-century SSP5-8.5 primarily affects lower-latitude portions of the current range (southern China, parts of the Mediterranean fringe, and the southern United States), driven by excessive warming and altered precipitation regimes. Conversely, the most extensive gains are concentrated at the northern margins, particularly in western Russia, southern Scandinavia, and southeastern Canada, where previously unsuitable cooler regions become climatically favorable.

These results highlight that, while P. lactiflora is likely to benefit from net range expansion in a warming world, the highest-emission pathway in the late 21st century introduces the greatest spatial reorganization and localized degradation of historically optimal habitat.

Spatial shifts in the centroid of suitable habitat zones

3.3.4

Centroid analysis provides a concise summary of the overall direction and magnitude of range displacement under climate change (Fordham et al., 2012). In this study, the current centroid of suitable habitat (>0.312) for P. lactiflora is located in western Romania (approximately 45.3°N, 22.0°E). Under all future scenarios, the centroid consistently shifts northeastward, crossing Moldova and reaching eastern Ukraine (Figure 12).

Shift map of the center of gravity in the suitable area of P. lactiflora under future climate scenario.

Despite the shared northeastward trajectory, the distance travelled varies markedly between emission pathways. Under the SSP5-8.5 scenarios, the centroid migrates considerably farther than under SSP2-4.5, with the longest displacement occurring in SSP5-8.5 2061–2080 (approximately 1,234.41 km from the present-day position). In contrast, centroid movement under both SSP2-4.5 time slices remains more modest (approximately 1011 km).

This pronounced difference indicates that higher-emission pathways impose stronger environmental pressure, forcing the core of the most favorable habitat to shift farther poleward and eastward in search of cooler and wetter conditions that match the species’ current climatic optimum. The shorter migration distances projected under SSP2-4.5 suggest that moderate mitigation could substantially reduce the geographic reorganization required for the persistence of high-quality habitat.

Discussion

4

Reliability and limitations of MaxEnt model prediction

4.1

The MaxEnt model developed here exhibited excellent performance (mean test AUC = 0.945 ± 0.001; mean TSS = 0.762 ± 0.018), confirming its strong discriminatory power and robustness for predicting the global distribution of P. lactiflora. These values compare favorably with those reported in other high-quality SDM studies of medicinal plants and temperate perennials, providing confidence that the identified environmental drivers and projected range shifts are biologically credible.

Nevertheless, several limitations inherent to the MaxEnt framework and our implementation should be acknowledged (Lissovsky et al., 2021). First, predictive accuracy remains contingent on the quality, resolution, and temporal relevance of the environmental layers. Although we used widely accepted, high-quality datasets (WorldClim 2.1, HWSD, etc.), their ~10 km resolution inevitably smooths local topographic and microclimatic variation that can be critical for a montane–lowland species such as P. lactiflora. Second, MaxEnt assumes that the realized niche is primarily shaped by the supplied abiotic predictors and that species–environment relationships remain relatively stable over time (Phillips and Dudík, 2008; Bradie and Leung, 2017). In reality, biotic interactions (e.g., competition with invasive species, herbivory, or mutualisms) and anthropogenic factors (e.g., land-use change, urbanization, and historical cultivation) also exert substantial influence on local presence and abundance. These processes were not explicitly incorporated and may lead to overestimation of suitability in heavily modified landscapes (particularly in Europe and eastern North America) or underestimation in regions where human management has historically facilitated persistence outside the strict climatic niche (Wisz et al., 2013; Gwitira et al., 2014).

Despite these constraints, the consistency between our modelled core areas and the known centers of traditional cultivation and naturalization lends strong support to the overall reliability of the projections. Future refinements could usefully integrate dynamic land-use layers, dispersal constraints, or ensemble approaches across multiple algorithms to further reduce uncertainty.

Ecological impacts of key environmental variables on P. lactiflora distribution

4.2

Jackknife analysis and permutation importance identified four variables that collectively dominate the realized niche of P. lactiflora: precipitation of the warmest quarter (bio18), mean temperature of the coldest quarter (bio11), temperature seasonality (bio4), and solar radiation in November (srad11). The overriding influence of these predictors highlights the species’ specialization for a temperate monsoon climate characterized by abundant summer rainfall, moderate winter cold, limited annual thermal amplitude, and adequate late-autumn irradiance.

The paramount importance of bio18 underscores the critical dependence of P. lactiflora on high moisture availability during the main growing and flowering season. This requirement aligns closely with the summer-monsoon precipitation regime of its native East Asian range and explains both its absence from Mediterranean winter-rainfall regions and its successful cultivation in similarly summer-wet areas of central Europe and eastern North America (Chetvertak et al., 2025).

The strong constraining effect of bio11 reflects the necessity of a pronounced but not extreme cold period for proper dormancy induction and subsequent spring regrowth, consistent with the species’ well-documented chilling requirement and frost tolerance (Pescador et al., 2018; Marković et al., 2022). Similarly, the negative response to high temperature seasonality (bio4) indicates low tolerance for continental climates with large annual thermal fluctuations, restricting the species to oceanic or sub-oceanic temperate zones (Menzel and Sparks, 2006).

November solar radiation (srad11) emerged as an unexpectedly influential photoperiod- and energy-related cue. High suitability within a relatively narrow late-autumn irradiance window (≈4,000–12,500 kJ m⁻² day⁻¹) suggests that sufficient light penetration before leaf senescence is essential for carbohydrate translocation to roots, thereby supporting overwinter survival and the following year’s flowering (Zhao et al., 2015).

These findings substantially deepen our understanding of the ecological amplitude of P. lactiflora and provide physiologically interpretable thresholds that can guide both cultivation practices and conservation planning under climate change. Nonetheless, the model necessarily simplifies reality; factors not explicitly included—such as soil microbial communities, atmospheric pollutants, or fine-scale edaphic heterogeneity—may further modulate local realized niches and warrant investigation in future studies.

This study demonstrates that integrating soil properties and solar radiation variables into MaxEnt significantly enhances explanatory power for temperate herbaceous perennials, particularly by capturing fine-scale edaphic constraints and photoperiodic cues often overlooked in purely climatic SDMs. The emergence of November solar radiation as a dominant driver highlights the potential of radiation data to refine niche predictions in mid-latitude species, where autumn light availability critically influences resource storage and overwintering success. These findings suggest that future SDM applications for medicinal plants and other temperate herbs should routinely incorporate solar radiation and soil layers to improve mechanistic understanding and predictive accuracy under climate change scenarios.

Global-scale ecological mechanisms underlying distribution patterns

4.3

The distribution of P. lactiflora at a global scale reflects a suite of interacting ecological mechanisms that operate across continental climate gradients, biogeographic barriers, and species-specific traits. As a temperate herbaceous perennial, P. lactiflora exemplifies how global climate patterns shape plant niches through interactions between abiotic drivers and biotic processes.

Precipitation during the warmest quarter (bio18) acts as the primary limiter, reflecting the species’ dependence on summer monsoon regimes typical of East Asia. This mechanism aligns with broader global patterns where herbaceous perennials in humid temperate zones rely on seasonal moisture pulses for rapid growth and reproduction, distinguishing them from woody species that tolerate greater drought through deeper roots (Wang et al., 2020; Xing et al., 2024). At continental scales, this creates disjunct distributions, with climatic analogues in central Europe and North America enabling naturalization, but biogeographic barriers (e.g., oceans, deserts) prevent natural dispersal from the native range, underscoring human-mediated introduction as a key mechanism for global expansion.

Winter temperature (bio11) and seasonality (bio4) further constrain poleward limits, enforcing a chilling requirement for dormancy break—a common mechanism in temperate perennials to synchronize growth with seasonal cycles and avoid frost damage (Lubbe et al., 2021). Globally, this interacts with latitudinal gradients, where increasing temperature variability in continental interiors excludes the species, favoring oceanic climates with buffered thermal regimes (Gillespie and Volaire, 2017). November solar radiation (srad11) introduces a photoperiodic dimension, likely regulating autumn resource allocation to roots, enhancing overwintering success in mid-latitude zones but limiting equatorward expansion where daylength cues are less pronounced (Petterle et al., 2013).

These abiotic mechanisms are modulated by biotic interactions at global scales, such as competition with invasive species in newly suitable areas or facilitation by soil microbes in authentic production regions—factors our model indirectly captures through edaphic variables but may underestimate in projections (Wisz et al., 2013). Climate change amplifies these dynamics, as seen in projected northeastward shifts driven by poleward migration of isothermal lines—a mechanism observed in many herbaceous taxa (Chen et al., 2011). However, dispersal limitations and habitat fragmentation could hinder realization of these projections, potentially leading to extinction debt in trailing edges (Tilman et al., 1994). Overall, P. lactiflora’s global patterns underscore how climate envelopes interact with evolutionary history and anthropogenic factors, offering a model for predicting medicinal plant vulnerabilities in a changing world.

Trends in the distribution of P. lactiflora under future climate scenarios

4.4

Our projections reveal a consistent expansion of suitable habitat for P. lactiflora under all examined scenarios, accompanied by a clear northeastward shift of the centroid of highly suitable areas. These twin patterns—range enlargement coupled with poleward and eastward displacement—are typical of many temperate East Asian perennials responding to warming and changing precipitation regimes.

The magnitude of both expansion and centroid shift is markedly greater under SSP5-8.5 than SSP2-4.5, particularly by 2061–2080. The longer migration distances projected under the high-emission pathway reflect stronger climatic pressure, forcing the core of optimal habitat to track cooler and wetter conditions farther into higher-latitude regions of eastern Europe and western Russia.

This climate-driven reorganization has dual implications for conservation and sustainable utilization. On one hand, the substantial northward extension creates new opportunities for cultivation in areas that are currently too cool, potentially broadening the resource base for Paeoniae Radix. On the other hand, localized loss of highly suitable habitat at the southern and western margins of the current range—especially under late-century SSP5-8.5—risks reducing population stability and genetic diversity in traditional production regions, while increasing competitive interactions with resident species in newly colonized northern zones.

These findings underscore the urgent need for proactive, forward-looking conservation strategies. Priority actions should include: (i) establishing monitored ex-situ collections from across the current native and cultivated range, (ii) identifying and protecting future climatic refugia identified in this study, and (iii) developing assisted migration or new cultivation guidelines for emerging high-suitability areas. Timely implementation of such measures will be essential to safeguard both wild germplasm and the long-term security of this economically and culturally important medicinal resource in a warming world.

Strategic suggestions for the sustainable use of P. lactiflora germplasm resources

4.5

The findings of this study provide clear, actionable guidance for the long-term conservation and sustainable utilization of P. lactiflora germplasm in a changing climate.

Firstly, in regions currently identified as highly suitable—particularly eastern China, the Korean Peninsula, central Europe, and the northeastern United States—conservation measures should be significantly strengthened. Establishing or expanding nature reserves, implementing stricter land-use regulations, and reducing habitat fragmentation will help protect genetically diverse wild and cultivated populations that have co-evolved with local conditions over centuries.

Secondly, for areas projected to become highly suitable in the future (especially western Russia, southeastern Canada, and southern Scandinavia under SSP5-8.5), proactive planning and ecological restoration are recommended. Creating connectivity corridors, improving soil conditions where needed, and conducting pilot assisted-migration programs using provenances from trailing-edge populations can facilitate successful range expansion and prevent future genetic bottlenecks.

Lastly, given the consistent northeastward shift of the suitability centroid, long-term dynamic monitoring networks should be established across the entire current and future range. At the same time, international cooperation among East Asia, Europe, and North America must be enhanced through shared databases, joint field surveys, and coordinated research initiatives to track real-time responses and refine adaptive management strategies effectively.

Implementing these integrated, forward-looking measures will be essential to ensure both the ecological persistence and continued medicinal availability of this globally significant species throughout the 21st century.

Conclusions

5

This study employed an optimized MaxEnt model coupled with ArcGIS spatial analysis to assess the global habitat suitability of P. lactiflora under present-day and future climate conditions. The excellent model performance (mean AUC = 0.945, TSS = 0.762) confirms its reliability for large-scale ecological forecasting.

Precipitation of the warmest quarter, mean temperature of the coldest quarter, temperature seasonality, and November solar radiation emerged as the dominant constraints on the species’ realized niche, collectively explaining its restriction to temperate summer-wet climates with moderate winters and stable annual thermal regimes.

Currently, highly suitable habitat is concentrated in East Asia, central Europe, and the northeastern and midwestern United States. All future scenarios (SSP2-4.5 and SSP5-8.5, 2041–2060 and 2061–2080) project a substantial expansion of suitable area, indicating that P. lactiflora is likely to benefit from moderate warming in the coming decades. However, under the high-emission SSP5-8.5 pathway by late century, core habitat quality declines in parts of the current range, accompanied by pronounced northeastward displacement of the suitability centroid.

These results highlight both opportunities and risks: while new cultivation regions will open at higher latitudes, traditional production areas face increasing climatic stress. Effective long-term conservation and sustainable utilization of P. lactiflora germplasm will therefore require (i) reinforced protection of existing high-suitability zones, (ii) proactive preparation of emerging suitable regions, and (iii) enhanced international monitoring and collaboration. Future refinements should prioritize incorporation of dynamic land-use layers, dispersal kernels, and biotic interactions (particularly competition and herbivory) to better distinguish realizable from fundamental niche space.