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.

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₃ 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.

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.

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.

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.

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 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.

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).

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.

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.

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.

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.

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.

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.

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.