We collected data on P. heterophylla distribution through the end of 2023 from the National Specimen Information Infrastructure (www.nsii.org.cn/2017/home.phpNSII), Chinese Virtual Herbarium (www.cvh.ac.cn), Plant Photo Bank of China (ppbc.iplant.cn), and literature indexed in the Chinese National Knowledge Infrastructure database. The “spThin” package in R (R Core Team, 2021) was used to delete data that were redundant or from locations within 10 km of a major sampling site in order to minimize variations due to sampling location (Aiello Lammens and Akçakaya, 2017). A total of 330 distribution records were retained in the final analysis ( Supplementary Table S2 ; Figure 1 ).

These data from previous sampling were supplemented with triplicate sampling of P. heterophylla from 44 plots in Anhui, Shandong, Fujian, Jiangsu and Guizhou during the harvesting season between July and August ( Supplementary Table S3 ; Figure 1 ).

The environmental data used in this study came from the Worldclim version 2.1 (http://www.worldclim.org/) for 19 environmental variables, which were selected based on their relevance to plant survival and growth (Fick and Hijmans, 2017; Viljanen et al., 2024). Data were downloaded for the following periods: current period (1970–2000), future 2040s (averages for 2021-2040), 2060s (averages for 2041-2060), 2080s (averages for 2061-2080), and 2100s (averages for 2081-2100) (Debella-Gilo and Etzelmüller, 2009). In addition, The data of 55 ecological factors was acquired from the following sources: the soil data were obtained from the World Soil Database (http://www.crensed.ac.cn/), which was established by FAO and IIASA, terrain factor data were collected from the international scientific data service platform (http://datamirror.csdb.cn/), climatic data is sourced from the World Climate Database (http://www.worldclim.org/), which interpolates recorded data from various meteorological stations worldwide between 1950 and 2000. It comprises 43 climate type data (including 12-month rainfall and average gas temperature, as well as 19 comprehensive climate factors), 8 soil type data (such as soil pH, cation exchange capacity, sand content, clay content, soil subcategory, effective water content grade, soil texture classification, and organic carbon content), three topographic data (elevation, slope degree and slope direction), and vegetation types ( Supplementary Table S4 ).

Predicted future climate data for China were downloaded from the Beijing Climate Center Climate System Model (BCC-CSM2-MR) from the sixth phase of the Coupled Model Intercomparison Project (Wu T. et al., 2018; Sang et al., 2020; Jin et al., 2021; Kumar and Sarthi, 2021; Wang et al., 2021). Data covering the same time periods were downloaded for the shared socioeconomic pathway (SSP) 1-2.6, corresponding to low compulsion; SSP2-4.5, corresponding to medium radiative forcing; and SSP5-8.5, corresponding to high compulsion. These SSPs were chosen because of their relatively high accuracy and resolution and because they take into account local development (He et al., 2017; Liu et al., 2021). To eliminate collinearity, we performed Pearson’s correlation analysis across environmental variables and Exclude high-correlation environmental variables (|r| ≥ 0.8), for factors with a correlation coefficient exceeding 0.8, only those with greater ecological significance were retained ( Supplementary Figures S1 , S2 ).

To ensure consistency with the data format requirements of the MaxEnt 3.4.1 model for environmental variables, resolution of all environmental variables was resampled to 2.5′(5 km×5 km) using ArcGIS 10.4 software (Esri.(2020).ArcGIS Desktop: Release 10.4. Redlands, CA: Environmental Systems Research Institute.).

Using the Biomod2 package (v4.3.1, https://cran.r-project.org/bin/windows/base/). Nine algorithms were applied in this study: SRE (Surface Range Envelope), RF (Random Forest), FDA (Flexible Discriminant Analysis), CTA (Classification Tree Analysis), GAM (Generalized Additive Model), ANN (Artificial Neural Network), GLM (Generalized Linear Model), MARS (Multivariate Adaptive Regression Splines), MaxEnt (Maximum Entropy) (Zhao Y. et al., 2021). ENMeval optimized feature class setting and β regularization multiplier are used, and the other algorithms use Biomod2 default model setting. The Biomod2 is widely used to predict species distribution and test species distribution model (Huang D. et al., 2023; Huang L. et al., 2023; Suicmez and Avcı, 2023) From the distribution data points in the model, 75% were randomly selected to train the model, and the remaining 25% were used to evaluate model performance in 10 replications (Muscarella et al., 2015). Algorithm performance was assessed in terms of average AUC (the area under the subject operating characteristic curve) and TSS (true skill statistic values) of single and ensemble model (Warren and Seifert, 2011).

The ENMeval package in R software was used to optimize the selective model. The ‘regularization multiplier’ and ‘feature combination’ modules in the ENMeval package within R were used to reduce model complexity and improve model accuracy (Deka and Morshed, 2018). The value of the regularization multiplier ranged from 0.5 to 4, and it was increased in steps of 0.5 for a total of eight frequency doublings. Nine feature combinations were applied: “L”, “LQ”, “H”, “LQH”, “LQHP”, “LQHPT”, “QHP”, “QHPT”, and “HPT”, where L was linear; Q, quadratic; H, hinge; P, product; and T, threshold. Altogether 72 parameter combinations were tested.

The spatial distribution of suitable habitats for P. heterophylla was output as raster data from ArcGIS 10.4.1 (Esri.(2020).ArcGIS Desktop: Release 10.4.Redlands,CA: Environmental Systems Research Institute.) with attributes of the setting layer defined using the “reclassify” tool according to the natural breakpoint method (Guiquan et al., 2023). Habitats were classified as unsuitable if their fitness value was 0-0.1; secondarily suitable, 0.2-0.25; suitable, 0.26-0.5; or optimal, 0.6-1.

In this study, the SDMTool package in R software was used to calculate the change trends in the suitable region and the geometric center positions of the suitable regions in the modern era and the future. We considered the suitable P. heterophylla habitat as a whole and reduced it to a vector particle and used the change of the position of the centroid to reflect the size and direction of the species’ suitable habitat. Finally, The center of mass was tracked using different SDMs to examine the centroid of P. heterophylla during various periods and under diverse climatic conditions, in order to evaluate the migration distance of the suitable P. heterophylla zone in latitude and longitude coordinates (Brown et al., 2017).

A multivariate environmental similarity surface was generated for the current distribution of suitable habitats as well as for future distributions predicted under different climate change scenarios in order to explore the degree of climatic anomaly in current and future distributions. The surface was generated relative to a reference layer based on the environmental variables in the current distribution area of P. heterophylla. The surface indicates the degree of similarity between a set of predictive variables (V1, V2, Vi,…) and a set of reference points. In the reference layer, mini and maxi are the minimum and maximum values of the environmental variable Vi; pi, the value of the environmental variable Vi at a certain point P during the given period; and fi, the percentage of points where Vi< pi (Muscatello et al., 2021). When fi = 0, the surface is 100 (p - mini)/(maxi - mini); when 0< fi ≤ 50, the surface is 2fi; when 50< fi< 100, the surface is 2 (100 - fi); and when fi = 100, the surface is 100 (maxi-100)/(maxi-mini). If the surface is negative, it indicates that the value of at least one variable falls outside the environmental range of the reference point during the given period, which is called the climate anomaly point. When the value of the point is 100, it indicates that the climate environment is completely consistent with the reference layer, and the climate is normal.

Environmental factors that most heavily influenced modeled changes in distribution of suitable habitats were identified based on the “most dissimilar variable” (Zhao G. et al., 2021) defined as the minimum similarity value for the surface across different environmental variables at point P. The most dissimilar variable is also the variable showing the greatest anomaly and therefore a variable likely to strongly influence the geographic migration of species. This variable was determined using the “novel density” (Yao et al., 2016) tool in the maximum entropy model.

A total of 44 samples of P. heterophylla from the sites indicated in ( Figure 1 ) were assayed for heterophyllin B according to the Chinese Pharmacopoeia.

They were assayed for polysaccharide content as follows. An aliquot of Pseudostellariae Radix powder (approx. 0.1 g) was accurately weighed, added to a 100 mL round-bottom flask, dissolved in 70 mL of 80% ethanol, heated in a water bath, refluxed for 30 min, and filtered while still hot. The residue and filter paper were placed in a clean flask, 80 mL of water was added, the suspension was soaked for 1 h in a water bath at 90 °C, it was filtered while it was still hot, the residue was washed three times with hot water (5 mL each time), and the washes were pooled together with the filtrate. The solution was allowed to cool and transferred to 100 mL volumetric flask. An aliquot (2.0 mL) was transferred into a 25 mL test tube and mixed well with 1 mL of 4% phenol solution, immediately after which 5 mL of concentrated sulfuric acid was added dropwise and the solution was shaken well each time. The solution was allowed to stand at room temperature.

To prepare a reference solution for polysaccharides, glucose was dried to constant weight (15.07 mg) in a 25 mL volumetric flask, and water was added to yield a final glucose concentration of 602.8 µg·mL^-1. A series of calibration solutions were prepared by transferring 0.5, 0.8, 1.0, 1.3, 1.5 and 1.8 mL of the reference solution into 25 mL measuring bottles, diluted with water and shaken well.

Absorbance of the reference solutions and Pseudostellariae Radix samples was measured at 489 nm. The calibration curve showed excellent linearity over the glucose concentration range from 3.0140 to 10.8505 µg·mL^-1 [absorbance = (0.0548 glucose concentration) + 0.0042, r = 0.9996]. Polysaccharide content was defined as

where c is the concentration of glucose (µg·mL^-1) among the polysaccharides in the sample, d is the polysaccharide dilution factor, f is the polysaccharide conversion factor (2.38), and w is the sample weight (g).

Associations of ecological factors with content of polysaccharides and heterophyllin B in Pseudostellariae Radix were explored using stepwise regression in SPSS 26.0(IBM, Chicago, IL, USA).

Based on 330 occurrence records and 19 environmental variables, the AUC and TSS values of nine species distribution niche models (SRE, RF, FDA, CTA, GAM, ANN, GLM, MARS, MaxEnt) were compared using the Biomod2 platform. The results are presented in Table 1 . Among the tested models, MaxEnt exhibited higher AUC values (0.969), followed by RF (0.956) and MARS (0.950). Similarly, MaxEnt demonstrated higher TSS values of 0.845 along with ANN (0.808) and RF (0.786). Notably, MaxEnt outperformed all other models with its largest AUC and TSS values indicating superior prediction accuracy. Furthermore, the current distribution predicted by this model showed the most accurate fit to the actual plant distribution ( Figure 2 ). Consequently, we selected the MaxEnt for subsequent analysis. This model was optimized by testing eight values of the regularization multiplier from 0.5 to 4 as well as nine combinations of parameters (L, LQ, H, LQH, LQHP, LQHPT, QHP, QHPT, HPT), giving a total of 72 combinations of multipliers and parameters. Applying a regularization multiplier of 1 and combination of features LQH gave an AICc = 0 as well as an AUC that was 76.38% lower and a 10% training omission rate that was 59.55% lower than with default parameters ( Figure 3 ). This optimized model was used in all subsequent analyses.

We utilized 55 ecological factors to screen for the primary determinants influencing the distribution of P. heterophylla ( Supplementary Table S5 ), and ascertained that all of them were climate-related variables. Consequently, in this investigation, we narrowed down our focus to 19 climatic ecological factors when screening for the principal determinants impacting the distribution of P. heterophylla.

The jackknife test was employed to evaluate the contribution of each ecological factor in the MaxEnt to the predicted distribution of P. heterophylla. Environmental factors with a contribution exceeding 4% were identified as primary factors. Six main factors were determined, collectively accounting for 87% of the distribution ( Table 2 ): Precipitation of driest month accounted for 31%; annual precipitation accounted for 18%; precipitation of driest quarter accounted for 15.8%; mean annual temperature accounted for 11.7%; temperature seasonality (standard deviation *100) accounted for 5.5%; mean temperature of coldest quarter accounted for 5%.

In order to better grow P. heterophylla, further investigations were carried out to obtain the optimum habitat values of the above six key environmental variables (Bio14, Bio12, Bio17, Bio1, Bio4, Bio11). A presence probability greater than 0.5, the same range with optimum habitat, was used to estimate the range of the optimum habitat value for each variable. Maximum entropy modeling for the current period showed increasing P. heterophylla occurrence with higher precipitation of driest month (variable Bio14), with the optimal precipitation ranging from 18.8116 to 41.2911 mm ( Figure 4A ). Probability of the plant’s occurrence exceeded 60% when annual precipitation (Bio12) was between 906 and 1394 mm ( Figure 4B ); precipitation of driest quarter (Bio17), between 64 and 164 mm ( Figure 4C ); and annual mean temperature (Bio1), between 13 and 16°C ( Figure 4D ). The optimal range of temperature seasonality (Bio4) was 651-935 ( Figure 4E ), while that of mean temperature of coldest quarter (Bio11) was 2.04-6.84°C ( Figure 4F ).

The MaxEnt was used to simulate the distribution map of suitable areas for P. heterophylla in the current period ( Figure 5 ). Using 55 environmental variables, our established model accurately predicted the predominant distribution of P. heterophylla optimum habitats in central Guizhou, northeastern Fujian, eastern Anhui, western Jiangsu, and southern Shandong ( Figure 5A ). Additionally, Using 19 climatic factor, we model and predict the optimum habitats of P. heterophylla were mainly distributed in the southwestern region (Guizhou) and eastern regions (Anhui, Zhejiang, Fujian, Jiangsu) of China ( Figure 5B ). The similarity in the prediction of suitable habitat based on 55 ecological factors and 19 climatic factors further supports the notion that climate factors predominantly influence the optimal growth conditions for P. heterophylla.

Maximum entropy modeling predicted that the area of habitat suitable for P. heterophylla would increase in the future, but this increase would be primarily in suitable and secondarily suitable habitats ( Figure 6 ). Under the SSP1-2.6 scenario, the area of habitat showing low suitability for the plant is predicted to increase by 0.033×10⁶-0.103×10⁶ km² during the four periods that we analyzed, while the area of optimal habitat is predicted to decrease by 0.004×10⁶ km² from 2081 to 2100. Under the SSP2-4.5 scenario, in contrast, the area of habitat showing low suitability is predicted to increase by 0.069×10⁶-0.403×10⁶ km² and the area of optimal habitat by 0.078×10⁶-0.135×10⁶ km² during the four periods. Similarly, under the SSP5-8.5 scenario, the area of habitat showing low suitability is predicted to increase by 0.117×10⁶-0.379×10⁶ km² and the area of optimal habitat by 0.059×10⁶-0.691×10⁶ km². In both of these scenarios, the largest increase in habitat of low suitability is predicted to occur in 2061-2080, and the largest increase in optimal habitat is predicted to occur in 2021-2040.

Future average multivariate similarities under different climate change scenarios ranged from 8.40 to 24.98 ( Figures 7A–L ). Scenario SSP1-2.6 showed the highest similarity and weakest climate anomaly, while scenario SSP5-8.5 showed the lowest similarity and strongest anomaly. The variables predicted to be most influential on future plant distribution, regardless of scenario, were annual mean temperature, minimal temperature in the coldest month, annual precipitation, precipitation in the driest quarter, and precipitation in the driest month ( Figures 8A–L ). The temperature-related variables were most influential over a much wider area than precipitation-related ones.

Relationships of the most influential ecological factors with levels of polysaccharides and heterophyllin B in Pseudostellariae Radix from 44 sites across China ( Table 3 ) were explored using step regression. Two regression equations emerged, indicating significant correlations between content of polysaccharides or heterophyllin B with several ecological factor.

(where Y₂ is Heterophyllin B content, X₁ is Bio6, X₂ is Bio2, X₃ is Bio3, X₄ is Bio4, X₅ is Bio12, X₆ is Bio 16. R² = 0.467, P = 0.000< 0.05.)

These equations show that the contents of two active compounds are mainly related to temperature. The ecological factors that impact polysaccharides include Bio 6, Bio1, Bio2, Bio3, Bio4, Bio12, Bio14 and Bio16. Among these factors, the primary influencers of polysaccharides are Min Temperature of Coldest Month (bio6), Annual Mean Temperature (bio1), Mean Diurnal Range (bio2), and Isothermality (Bio2/Bio7) (*100) (bio3). The ecological factors that influence Heterophyllin B also encompass Bio 6、Bio2、Bio3、Bio4、Bio12 and Bio16. Among these factors, the key determinants affecting Heterophyllin B are Min Temperature of Coldest Month (Bio6), Mean Diurnal Range (Bio2), and Isothermality (Bio2/Bio7) (*100) (bio3).

The spatial analysis function of ArcGIS was used to explore the future distribution of polysaccharide content of Pseudostellariae Radix under different climate change scenarios ( Figures 9A–L ). Regardless of the scenario, high polysaccharide content was predicted for P. heterophylla in Hebei, Shaanxi, Shandong, Jiangsu, Anhui, Zhejiang, Jiangxi, Henan, northeast Hubei, northwest Hunan, and Sichuan ( Figures 9A–L ). These areas are predicted to be larger than the current distribution ( Figures 9M–P ). Although there were variations in the areas with high polysaccharide content across different time periods under future climate scenarios, the total area of these regions showed minimal disparity.

In contrast, high heterophyllin B content was predicted for different regions depending on the climate change scenario and time period. High content was located mainly in Shandong, Guizhou, western Hunan and Hubei under the SSP1-2.6 scenario, while it was located mainly in Hubei and Shaanxi under the SSP2-4.5 scenario for the period from 2041 to 2060. In contrast, high content was located mainly in Henan and northern Hubei under the SSP5-8.5 scenario for the period from 2081 to 2100. This modeling suggests that the area covered by high heterophyllin B content will increase in the future, especially under the SSP5-8.5 scenario ( Figure 10 ).

In this study, There were large differences in suitable range shift locations of P. heterophylla under different climate scenarios. We found that the center of mass suitable P. heterophylla habitat and the center of mass levels of bioactive components are not in the same area. Under the future climate scenarios SSP1-2.6, SSP2-4.5, the migration direction of the center of mass in the suitable area is southwest, which is located in the north of Hubei Province. Under the future climate scenario SSP5-8.5, the direction of the center of mass migration in the suitable area is the north. Under the three Climate Scenarios, Migration with a migration distance of 295 ~ 669 km. The migration trend of polysaccharide centroid is consistent with that of suitable area, and the migration distances are 411 km, 601 km and 585 km. The center of mass of Heterophyllin B content was predicted to shift to the north from the current area in the future by 445 km under SSP1-2.6, 534 km under SSP2-4.5, and 912 km under SSP5-8.5. Under the future climate scenarios SSP1-2.6 and SSP2-4.5, the center of mass Heterophyllin B and the center of mass of the suitable area move in different directions. Under the future climate scenario SSP5-8.5, the migration direction of the center of mass of the suitable area aligns with that of polysaccharide centroid and Heterophyllin B centroid ( Figure 11 ).