The 20 occurrences of A. adenophora were obtained from the Global Biodiversity Information Facility (GBIF, http://www.gbif.org/), the 56 occurrences of A. adenophora were obtained from Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn), the 80 occurrences of A. adenophora were obtained from National Specimen Information Infrastructure (NSII, http://www.nsii.org.cn), and the 2620 occurrences of A. adenophora were obtained from our field survey data. Finally, the total of 2776 occurrences of A. adenophora in China were obtained from online databases and our field survey data in China. We thereafter selected the occurrence data in the database since 1940, andremoved the occurrence data that did not have detailed geographic information. ENMTools was used to screen the occurrence data so that only one distribution point was retained within each 5 km × 5 km raster (Warren et al., 2010). Finally, we obtained 2,484 occurrences of A. adenophora for modeling the PGD in China using SDMs (Figure 1).

We obtained the current climate data from the PaleoClim database (http://www.paleoclim.org) with a resolution of 2.5’(Brown et al., 2018). The dataset of parameters based on high-resolution monthly temperature and precipitation values for 1979-2013 is broadly used in ecological applications worldwide. Future climate data, and altitude data were downloaded from the WorldClim database version 1.4 (http://www.worldclim.org) (Hijmans et al., 2005). The current and future climate data included 19 bioclimatic variables (Table S1). The future climate data with the Community Climate System Model (CCSM) 4.0 for the 2050s and the 2070s included two representative concentration pathways (RCP2.6 and RCP8.5). Multicollinearity between the bioclimatic and altitude variables can result in model overfitting. Therefore, we used ENMTools to calculate the correlation coefficients of bioclimatic and altitude variables (Warren et al., 2010). For two environmental variables with a correlation coefficient greater than 0.8 (| r| > 0.8), the more meaningful one was retained (Supplementary Figure 1). Finally, eight environmental variables were retained: mean diurnal range (bio2), isothermality (bio3), minimum temperature of the coldest month (bio6), temperature annual range (bio7), mean temperature of the warmest quarter (bio10), precipitation seasonality (coefficient of variation) (bio15), precipitation of the coldest quarter (bio19), and altitude.

In the present study, based on the occurrence data and environmental variables, we modeled the PGD of A. adenophora under climate change using eight different algorithms (GLM, GBM, CTA, ANN, SRE, FDA, RF, and MaxEnt) in the biomod2 package (Thuiller et al., 2016) in R Studio 4.1.3. For each modeling process, 75% of the distribution records were randomly selected as the training dataset while the remaining 25% were used as the testing dataset. The division of the training and testing data was repeated five times. In addition, 1,000 pseudoabsence points were randomly selected; this process was repeated three times. Finally, 120 models were obtained. After assessing the accuracy of the eight models, we selected four individual models with average values exceeding 0.9 to construct an EM. We then predicted the PGD of A. adenophora under various climate change conditions using the constructed EM.

In general, the true skill statistic (TSS), area under the receiver operating characteristic (ROC) curve (AUC), and Cohen’s Kappa (KAPPA) values are used to evaluate model accuracy. TSS considers the missing average error and is unaffected by the size of the verification dataset (Allouche et al., 2006). The ROC curve is an acceptance curve, with the horizontal coordinate indicating the false-positive rate (1-specificity) and the vertical coordinate indicating the true-positive rate (1-omission rate) (Peterson et al., 2008). KAPPA is affected by the incidence of occurrence data and thresholds, and determines the accuracy of predictions relative to random occurrences (Monserud and Leemans, 1992). Generally, the higher the values of these three indicators, the higher is the accuracy of the species distribution model results.

From the results of the EM, the ASCII raster layer was generated based on the probability (P) of the presence of A. adenophora, which ranged from 0 to 1,000. Finally, the PGD of A. adenophora was classified into four types using ArcGIS: high-suitability habitat (500 ≤ P ≤ 1000), moderate-suitability habitat (300 ≤ P < 500), low-suitability habitat (100 ≤ P < 300), and unsuitable habitat (0 ≤ P < 100).

The shifts in the distributional centroids of highly suitable habitat were used the species distribution model tool (SDM tool, http://www.sdmtoolbox.org) to detect between current and future climate. SDM tool can be integrated into ArcGIS software for calculating, which can analyze the shifts in the distributional centroids under different climate scenarios (Brown, 2014). The threshold of distributional centroids was the threshold of high-suitability habitat.

In the present study, we calculated the evaluation accuracy of various species distribution models, including GLM, GBM, CTA, ANN, SRE, FDA, RF, MaxEnt, and EM. For A. adenophora, our results indicated that GLM, GBM, FDA, and RF were more accurate than CTA, ANN, SRE, and MaxEnt (Figure 2). The mean TSS values for GLM, GBM, FDA, and RF were 0.891, 0.892, 0.886, and 0.895, respectively. The mean AUC values for GLM, GBM, FDA, and RF were 0.961, 0.963, 0.950, and 0.961, respectively. The mean KAPPA values for GLM, GBM, FDA, and RF were 0.926, 0.925, 0.921, and 0.927, respectively. Therefore, we selected these four models to construct the ensemble model (EM). The mean TSS, AUC, and KAPPA values for the EM were 0.925, 0.993, and 0.936, respectively. Our results indicated that the EM improved fitting accuracy and reduced fitting uncertainty. In other words, the PGD of A. adenophora predicted using EM was more reliable than that predicted using the single models.

We estimated the contribution values of each environmental variable using GLM, GBM, FDA, and RF, and analyzed the variables that significantly affected the PGD of A. adenophora. Our results indicated that the top three environmental variables with the highest mean contribution values were the temperature annual range (bio7, 0.449), mean diurnal range (bio2, 0.302), and minimum temperature of the coldest month (bio6, 0.297) in the GLM, Fugure 3A; the temperature annual range (bio7, 0.522), minimum temperature of the coldest month (bio6, 0.061), and precipitation of the coldest quarter (bio19, 0.042) in the GBM, Fugure 3B; the minimum temperature of the coldest month (bio6, 0.501), temperature annual range (bio7, 0.321), and precipitation of the coldest quarter (bio19, 0.072) in the FDA, Fugure 3C; and the temperature annual range (bio7, 0.320), minimum temperature of the coldest month (bio6, 0.138), and mean diurnal range (bio2, 0.034) in the RF (Fugure 3D). In summary, the significant environmental variables affecting the PGD of A. adenophora were the mean diurnal range (bio2), minimum temperature of the coldest month (bio6), temperature annual range (bio7), and precipitation of the coldest quarter (bio19).

We predicted the PGD of A. adenophora in China under the current climate based on nine SDMs in biomod2 (Figure 4). Our results indicated that the PGD of A. adenophora predicted using the GLM, GBM, FDA, RF, and EM was more consistent with the actual geographic distribution than that predicted using the other SDMs. The accuracy of the EM based on TSS, AUC, and KAPPA was higher than that of GLM, GBM, FDA, and RF. Therefore, we selected the PGD of A. adenophora predicted using the EM to analyze its PGD in China under the current climate (Figure 4I). The PGD of A. adenophora under the current climate was mainly distributed in southwestern and southeastern China under different climate zones, such as the tropical edge, north subtropical, middle subtropical, and south subtropical zones. The total, high, and moderate suitability habitat areas of A. adenophora in China were 153.82 × 10⁴, 92.13 × 10⁴, and 21.04 × 10⁴ km², respectively, accounting for 16.02, 9.60, and 2.19% of the Chinese mainland area, respectively.

The PGD of A. adenophora under RCP2.6 and RCP8.5 in the 2050s and the 2070s are presented in Figure 5. The gain, loss, and stable areas of highly suitable habitat were shown in Supplementary Figure 2. The PGD of A. adenophora under the future climate change scenarios was also distributed in Southwest and Southeast China. The range of high-suitability habitats for A. adenophora under the current climate was quite different from that under the future climate. The highly suitable habitat in Southwest China, which was distributed at the edge of the tropical and middle subtropical zones, shifted to moderately suitable habitat under the future climate change.

The low-suitability habitat areas of A. adenophora for the 2050s and 2070s were 36.57 × 10⁴ and 34.24 × 10⁴ km², respectively, under RCP2.6 (accounting for 3.81 and 3.57% of the Chinese mainland area, respectively), and 35.41 × 10⁴ and 40.59 × 10⁴ km², respectively, under RCP8.5 (accounting for 3.69, and 4.23% of the Chinese mainland area, respectively). The low-suitability habitat areas showed a decreasing trend from the current to the 2050s and the 2070s, and decreased significantly under RCP2.6 in the 2070s (Figure 6).

The moderate-suitability habitat areas of A. adenophora for the 2050s and 2070s were 49.94 × 10⁴ and 52.11 × 10⁴ km², respectively, under RCP2.6 (accounting for 5.20 and 5.43% of the Chinese mainland area, respectively), and 55.87 × 10⁴ and 55.6 × 10⁴ km², respectively, under RCP8.5 (accounting for 5.82 and 5.79% of the Chinese mainland area, respectively). The moderate-suitability habitat areas showed an increasing trend from the current to the 2050s and the 2070s, and increased significantly under RCP8.5 in the 2050s (Figure 6).

The high-suitability habitat areas of A. adenophora for the 2050s and 2070s were 31.42 × 10⁴ and 31.92 × 10⁴ km², respectively, under RCP2.6 (accounting for 3.27 and 3.33% of the Chinese mainland area, respectively), and 50.4 × 10⁴ and 61.12 × 10⁴ km², respectively, under RCP8.5 (accounting for 5.25 and 6.37% of the Chinese mainland area, respectively). The high-suitability habitat areas showed a decreasing trend from the current to the 2050s and the 2070s, and decreased significantly under RCP2.6 in the 2050s (Figure 6).

The total-suitability habitat areas of A. adenophora for the 2050s and 2070s were 117.93 × 10⁴ and 118.27 × 10⁴ km², respectively, under RCP2.6 (accounting for 12.28 and 12.32% of the Chinese mainland area, respectively), and 141.68 × 10⁴ and 157.31 × 10⁴ km², respectively, under RCP8.5 (accounting for 14.76 and 16.39% of the Chinese mainland area, respectively). Compared with the current climate, the total-suitable habitat areas in the 2050s under RCP2.6 and RCP8.5, and the 2070s under RCP2.6, showed a decreasing trend, and decreased significantly in the 2050s under RCP2.6. Under RCP8.5, the total suitability habitat area in the 2070s increased compared to that in the current climate.

The centroid of the high-suitability habitat was determined to be located at the point with coordinates 104.35°E, 26.77°N under the current climate. In the 2050s, under RCP2.6 and RCP8.5, the centroid of the high suitability habitat shifted from its current position to 105.83°E, 26.95°N, and 105.85°E, 27.17°N, respectively, shifting 144.33 and 24.86 km, respectively. In the 2070s, under RCP2.6 and RCP8.5, the centroid of the high suitability habitat shifted from its current position to the points with coordinates 106.81°E, 27.23°N and 107.08°E, 27.13°N, respectively, shifting 249.44 and 30.15 km, respectively. Therefore, the centroid of the highly suitable habitat of A. adenophora tended to shift to the northeast and southeast and to higher latitudes under future climate change (Figure 7).

SDMs are frequently being used to predict the PGD of IAPs; however, prediction results based on non-equilibrium occurrence data can lead to inaccurate model outputs (Jiménez-Valverde et al., 2008; Briscoe Runquist et al., 2019). The quantity and quality of input data are still crucial issues in terms of the reliability of predictive distributional modeling efforts (Welk, 2004). The non-equilibrium occurrence data of IAPs do not include data representative of all climatic tolerances. When IAPs are still in the process of biological invasion, the background points cannot be used in any SDM to classify locations beyond the existing range of the IAPs as (pseudo) absence data (Václavík and Meentemeyer, 2012). This is because such regions do not represent the absence of IAPs; rather, the IAPs have yet to occupy such regions. Therefore, we should use the equilibrium occurrence data to construct reliable SDMs for predicting the PGD of IAPs. For example, a study on the PGD of Lythrum salicaria in North America based on the currently known occurrence data and the occurrence data recorded before 1988 indicated that it was unreasonable to predict the PGD of this species based on the occurrence data recorded before 1988 (Welk, 2004). In the present study, we used the equilibrium occurrence data of A. adenophora (the occurrence data from online databases and our field survey data until 2021) and the EM constructed using GLM, GBM, and FDA to predict the PGD of A. adenophora under climate change in China. The mean TSS, AUC, and KAPPA values for the EM were 0.925, 0.993, and 0.936, respectively, indicating that the prediction results of this model were accurate. Our findings provide a framework for preventing and controlling the spread of A. adenophora under climate change conditions in China.

In our study, the PGD of A. adenophora was mainly distributed in southwestern and southeastern China, which is also located in the tropical and subtropical zones. The minimum temperature of the coldest month in the PGD of A. adenophora is above −3°C and the annual precipitation is usually between 800 and 1500 mm (Zhen-Feng et al., 2013). Our results indicated that temperature variables, such as the mean diurnal range, minimum temperature of the coldest month, and annual temperature range, were the most important variables affecting the PGD of A. adenophora in China, followed by precipitation (precipitation of the coldest quarter) and altitude. The successful invasion of A. adenophora can be mainly attributed to the synergy of temperature, precipitation, and altitude variables. Suitable environmental variables, including temperature, precipitation, and light, played key roles in the survival and growth of A. adenophora in the invasive zones, and the individuals from invasive populations showed greater plant height, stem diameter, leaf length, leaf breadth, and leaf area than those from native populations (Feng et al., 2009). The epigenetic modifications can lead A. adenophora to spread into the cooler zones of China (Zhu et al., 2007). For instance, A. adenophora preferred lower temperatures in the invasive areas compared with in the native areas (Tererai and Wood, 2014).

Previous studies have indicated that the winter temperature and minimum temperature of the coldest month have a significant impact on the distribution of A. adenophora in China and Nepal (Poudel et al., 2020; Fang et al., 2021). In general, A. adenophora prefers to grow in warm temperature conditions, such as between 10–25°C (Sun et al., 2004; Tererai and Wood, 2014). Compared with the native surroundings, experimental warming enhanced the biomass output and canopy of A. adenophora and decreased its mortality (He et al., 2012). In invasive zones, epigenetic modifications can facilitate the spread of A. adenophora to cooler regions, such that it prefers lower temperatures than normal (Huang et al., 2009). Therefore, temperature changes have a significant effect on the distribution patterns of A. adenophora. In brief, climate warming promotes environmental tolerance of A. adenophora and accelerates its invasion. The findings of these studies further verify the accuracy of our findings. Our findings showed that the precipitation of the coldest quarter also affected the PGD of A. adenophora. Climate change has altered precipitation patterns globally, while also profoundly affecting the invasion process of IAPs, such as increasing the invasibility of the community (Kreyling et al., 2008). Some evidence has shown that shifts in precipitation patterns attributed to climate change have accelerated the invasion process of IAPs such as Bromus hordeaceus in the arid habitats of North America (Hobbs et al., 2007). In the present study, we selected only the bioclimatic and altitude variables to predict the PGD of A. adenophora in China because A. adenophora has strong adaptability and low requirements for soil conditions. As long as the temperature and precipitation conditions are suitable, they flourish (Cao et al., 2004).

Some studies have indicated that the distribution areas of the IAPs have expanded or contracted under climate change. For example, the PGD of 11 IAPs is likely to expand in China under different climate change scenarios (Guan et al., 2020). The PGD of 108 IAPs is expected to decrease in Brazil in the future (Fulgêncio-Lima et al., 2021). In the 2050s and the 2070s, the total suitable habitat of A. adenophora in China is projected to decrease under RCP 2.6. The total suitable habitat areas of A. adenophora will increase under RCP 8.5 in the 2070s compared with that under the current climate. Some studies have indicated that the PGD of A. adenophora will expand by 45.3% and 2% under RCP 8.5 in India and Nepal, respectively (Poudel et al., 2020; Chaudhary et al., 2021). Our findings were consistent with this pattern.

Under climate warming, some IAPs will shift toward higher altitudes and latitudes (Petitpierre et al., 2016; Shrestha and Shrestha, 2019). The spread of IAPs is a result of synergy between abiotic and biotic variables. The interplay of multiple factors contributing to the rapid spread of A. adenophora in China has also been mainly attributed to the interaction of multiple factors. For instance, owing to the small size and light weight of its seeds, A. adenophora disperses over long distances through water and wind (Wang et al., 2011). Previous studies have indicated that flooding facilitates the invasion of A. adenophora into China. In particular, the great flood in 1998 along the Jinshajiang-Yangtze River significantly increased the spread rate of A. adenophora in China (Wang et al., 2011; Horvitz et al., 2014). The mechanism of the synergy of multiple factors increased the risk of invasion and spread in China. The findings of the present study indicate that the PGD of A. adenophora will spread to southeastern China under climate change conditions. The degree of suitability and invasion risk in southwestern China will decrease under climate change. This seems to be due to the saturated geographical distribution of A. adenophora in southwestern China. The high-suitability habitat of A. adenophora indicated a trend of spreading to southeastern China as a result of climate warming. Therefore, more attention should be paid to prevent its further spread in southeastern China. Our findings indicate that early warning and monitoring measures for A. adenophora in southeastern China should be formulated. In the present study, we estimated the distribution pattern and invasion risk areas under future climate scenarios. The results strongly support the need for developing appropriate and efficient strategies to prevent or minimize the invasion and spread potential in China.