Predicting the Potential Geographic Distribution of Invasive Freshwater Apple Snail Pomacea canaliculate (Lamarck, 1819) under Climate Change Based on Biomod2

Abstract

Pomacea canaliculata is widely distributed in the Chinese provinces south of the Yangtze River, causing serious damage to aquatic ecosystems, rice cultivation, and human health. Predicting the potential geographic distributions (PGDs) of P. canaliculata under current and future climate conditions in China is crucial for developing effective early warning measures and facilitating long-term monitoring. In this study, we screened various species distribution models (SDMs), including CTA, GBM, GAM, RF, and XGBOOST, to construct an ensemble model (EM) and then predict suitable habitats for P. canaliculata under current and future climate scenarios (SSP1-26, SSP2-45, SSP3-70, SSP5-85). The EM (AUC = 0.99, TSS = 0.96) yielded predictions that were more precise than those from the individual models. The Annual Mean Temperature (Bio1) and Precipitation of the Warmest Quarter (Bio18) are the most significant environmental variables affecting the PGDs of P. canaliculata. Under current climate conditions, the highly suitable habitats for P. canaliculata are primarily located south of the Yangtze River, collectively accounting for 17.66% of the nation’s total area. Unsuitable habitats predominate in higher-latitude regions, collectively covering 66.79% of China’s total land area. In future climate scenarios, the total number of suitable habitats for P. canaliculata is projected to expand into higher latitude regions, especially under SSP3-70 and SSP5-85 climate conditions. The 4.1 °C contour of Bio1 and the 366 mm contour of Bio18 determine the northernmost geographical distribution of P. canaliculata. Climate change is likely to increase the risk of P. canaliculata expanding into higher latitudes.

1. Introduction

Biological invasion refers to the process by which organisms encroach upon a new environment, either naturally or artificially, resulting in economic losses or ecological disasters for biodiversity, agriculture, forestry, animal husbandry, fishery production, and human health in the invaded area [1,2]. The majority of invasive species exhibit pronounced reproductive capacities and high tolerance, enabling them to surpass native species in the competition for resources and habitat [3]. Consequently, once an invasive alien species (IAS) has established a population in a new habitat, eradicating it becomes extremely challenging [4,5]. With the escalation of economic and trade interactions, as well as international exchanges between China and other nations, the number of IAS has been increasing in China at an unprecedented rate, representing a significant threat to China’s agricultural development, ecological environment, and public health [6]. As of 2016, China had confirmed 560 IAS, of which 125 species are pests, including 92 that adversely impact agricultural ecosystems. IAS are estimated to incur economic losses exceeding $18.9 billion annually [7].

Pomacea canaliculata (Lamarck, 1819), Ampullariidae, Architaenioglossa, and Gastropoda, originating from the Amazon Basin in South America, was listed by the International Union for Conservation of Nature (IUCN) in 2000 as one of the 100 most destructive invasive species in the world [8]. In the 1980s and early 1990s, P. canaliculata initiated their invasion of numerous countries in Southeast Asia and East Asia [9,10]. As of the latest data, in Southeast Asian countries like Thailand, the Philippines, Malaysia, and Vietnam, P. canaliculata has emerged as one of the most significant rice pests, resulting in over $1 million expended annually for their control [11,12,13,14]. Concurrently, P. canaliculata exerts significant detrimental effects on aquatic vegetation [15]. Its high population density resulted in the near-complete eradication of aquatic vegetation and led to the eutrophication of water bodies and increased phytoplankton biomass [15,16]. This culminates in water turbidity and the disturbance of ecosystem function [17]. Additionally, P. canaliculata serves as the intermediate host for Angiostrongylus cantonensis, with up to 6000 parasites detectable in a single snail [18]. Consumption of raw or undercooked snail meat can readily lead to eosinophilic meningitis, as well as myelitis and spinal meningitis [19]. Moreover, the empty shells of deceased snails represent a health hazard, as they may cause injuries to the feet of individuals involved in planting, harvesting, or managing crops [20,21].

Pomacea canaliculata demonstrates characteristics of amphibiousness, euryphagy, remarkable adaptability, and high fecundity, enabling it to rapidly establish substantial populations in invaded areas [22]. After its introduction into Zhongshan City, Guangdong Province, in 1981 as a valuable species for freshwater fisheries, P. canaliculata subsequently expanded its range to include Hubei, Zhejiang, Fujian, Jiangsu, Shanghai, Jiangxi, and other provinces through artificial introduction [23]. In 2003, P. canaliculata was included in the inaugural “List of invasive alien species in China”. It then disseminated over long distances in China primarily by attaching its eggs to ships or serving as prey for carnivorous fish (herring) or leeches [24,25]. Over the past two decades, Chinese government authorities have engaged in extensive efforts to prevent and control P. canaliculata on a large scale. Although there has been some reduction in the damage and spread of P. canaliculata, the overall impact remains relatively limited. By 2023, P. canaliculata had spread as far north as Chaoyang City, Liaoning Province, and as far west as Lanzhou City, Gansu Province [26]. Although low temperatures may limit P. canaliculata’s spread to higher latitudes, global climate warming aids its survival further north in the future [27,28,29]. Consequently, given the condition of global climate warming, predicting the potential geographic distributions (PGDs) of P. canaliculata could be a more effective strategy for preventing and controlling its invasion in non-endemic areas.

Species distribution models (SDMs) utilize known distribution locations of species and related environmental data to ascertain the ecological needs of species and then forecast the potential distribution of species through mathematical-statistical methods or machine learning theory [30]. SDMs are extensively employed in biodiversity conservation research, species invasion studies, evolutionary biology, and the prediction of species’ suitable habitats [31,32]. However, various SDMs exhibit differences in data requirements and algorithms, which leads to variations in the impact of environmental factors on species distribution, thus introducing uncertainty in the simulation results [33]. For instance, Classification Tree Analysis (CTA) [34], Random Forest (RF) [35], and eXtreme Gradient Boosting (XGBoost) [36] can process nonlinear and complex data, offering significant advantages like no need for data preprocessing, automatic feature value selection, and the ability to handle missing values. However, these models have complex structures, complicating the explanation of the relationship between variables and predictions, and are prone to overfitting. Conversely, the Generalized Linear Model (GLM) [37] and Generalized Additive Model (GAM) [38] excel in flexibly handling various data types, are well-suited for regression and classification problems, and offer more interpretable results. However, these models are sensitive to outliers and limited in processing high-dimensional data. MaxEnt, widely utilized in species distribution prediction, excels in modeling uncertainty and managing small sample data. However, it is highly sensitive to the quality of training data and requires a complex parameter adjustment process [39]. Therefore, leveraging the strengths of different models while mitigating their shortcomings has become a focal point in research. Ensemble models (EM) preserve the authentic data information captured in each individual model and subsequently minimize the errors inherent in the model, thus enhancing the reliability of the modeling results [40,41,42]. BIOMOD2 (BIOdiversity MODeling) has the capability to select superior SDMs for constructing EM, thereby enabling a more precise simulation of PGDs and permitting the deduction of the primary environmental factors influencing species distribution [43]. Due to BIOMOD2’s ability to compensate for the limitations of individual models, it has rapidly emerged as the predominant model for predicting species distribution in recent years [44,45,46].

Since the introduction into China of P. canaliculata in 1981, numerous researchers have employed various SDMs such as Maxent [47], CLIMEX [48], and GARP [49] to predict the suitable habitat of P. canaliculata in China. However, thus far, no researchers have employed Biomod2 to predict the suitable habitat of P. canaliculata in China under the conditions of future climate change. Therefore, in this study, we hypothesize that Biomod2 could provide more precise prediction of PGDs for P. canaliculata in China than traditional SDMs and that the PGDs of P. canaliculata may continue to spread toward higher latitudes in northern China under future warming climate conditions of P. canaliculata. In order to identify the important environmental variables that affect the spatial distribution of P. canaliculata and predict the distribution pattern and change law of P. canaliculata in China under future climate change conditions, We (1) selected five optimal SDMs to combine with an EM; (2) analyzed key environmental variables influencing PGDs of P. canaliculata; (3) predicted PGDs of P. canaliculata under current climate and climate change conditions in the 2021–2040, 2041–2060s, 2061–2080 and 2081–2100 by EM model; and (4) analyzed centroid shifts of P. canaliculata under current climate conditions versus climate change conditions.

2. Materials and Methods

2.1. Occurrence Records of P. canaliculata

The global distribution data of P. canaliculata are obtained from the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/, accessed on 1 December 2023), Mediterranean Plant Protection Organization Global database (EPPO; https://gd.eppo.int/, accessed on 1 December 2023), CABI Compendium Invasive Species (https://www.cabidigitallibrary.org/product/qi, accessed on 1 December 2023). The distribution data of snails in China were obtained by searching the China National Knowledge Infrastructure (CNKI; https://www.cnki.net/, accessed on 15 December 2023) database. A total of 107 articles, ranging from 1980 to the present, focusing on the birthplace of the P. canaliculata, were identified. A total of 5001 original records were obtained, of which 3001 records were selected through screening and deduplication using R for subsequent data processing and analysis (Supplementary S1). Due to the lack of data on actual absence points, we followed the methodologies outlined by Zhao (2022) [50] and Crego (2022) [51], selecting 1000 coordinate points at random from the global geographical scope to serve as pseudo-absence points [52] (Supplementary S2).

2.2. Environmental Variables Related to P. canaliculata

Given that P. canaliculata are aquatic organisms, environmental factors like soil conditions and wind speed do not influence their PGDs. Consequently, to enhance the efficiency of the model’s operation, 19 bioclimatic variables (downloaded from the WorldClim Database (www.worldclim.org/data/worldclim21.html, accessed on 18 December 2023)) [53], along with three terrain factors (elevation, aspect, and slope), downloaded from the Geospatial Data Cloud (www.gscloud.org, accessed on 18 December 2023), and water bodies coverage (WBG), downloaded from the ASTER Global Water Bodies Database (https://search.earthdata.nasa.gov/search, accessed on 1 December 2023), were identified as the primary environmental factors. A total of 19 future bioclimate variables were obtained from BCC-CSM2-MR, a climate model from the Coupled Model Intercomparison Project Phase 6 (CMIP6). These variables are derived from the latest shared socioeconomic pathways (SSPs) and representative concentration pathways (RCPs) to reflect various greenhouse gas emission scenarios. These scenarios encompass SSP1-26 (low forcing), SSP2-45 (moderate forcing), SSP3-70 (moderate to high forcing), and SSP5-85 (high forcing) for the periods of the 2021–2040, 2041–2060, 2061–2080, and 2081–2100. Correlation analysis was conducted using R to address the potential issue of multicollinearity among 19 bioclimate variables, which may result in model overfitting and affect prediction accuracy. Whenever the correlation coefficient between environmental variables exceeded 0.8, variables demonstrating a higher contribution rate to the distribution of P. canaliculata were prioritized for inclusion. Finally, seven bioclimate indicators (Bio1, Bio4, Bio7, Bio12, Bio17, Bio18, Bio19), one terrain indicator (elevation), and Water Bodies Coverage (WBC) were employed as the environmental variables. Additionally, the study assumed that elevation and WBC would remain constant over time to guarantee model comparability.

2.3. Construction and Evaluation of SDMs

Biomod2 is an integrated platform that establishes the relationship between species and environmental variables using statistical and machine learning algorithms. Biomod2 (version 4.2.1) includes 11 SDMs, such as the Generalized Linear Model (GLM, Generalized Boosting Model (GBM), Generalized Additive Model (GAM), Multivariate Adaptive Regression Splines (MARS), Classification Tree Analysis (CTA), Artificial Neural Networks (ANN), Surface Range Envelope (SRE, commonly known as BIOCLIM), Flexible Discriminant Analysis (FDA), Random Forest (RF), Maximum Entropy (MaxEnt) and eXtreme Gradient Boosting (XGBOOST). In each process, the distribution data of P. canaliculata were randomly divided into two subsets. One subset, comprising 80% of all distribution data, was used to train the model, while the other, containing 20% of the data, was used to evaluate the model’s performance. The Area under the ROC Curve (AUC) and True Skill Statistic (TSS) were utilized to assess the model’s accuracy. As the AUC value is unaffected by the diagnostic threshold and remains insensitive to species occurrence distribution rate, it enables a comprehensive comparison of the accuracy of diagnostic tests [54]. AUC is currently recognized as one of the most effective evaluation indicators for SDMs [55]. TSS was proposed by Allouche and Kadmon in 2006. TSS not only inherits the advantages of Kappa statistics but also addresses the limitation of Kappa statistics related to unimodal curve responses to species incidence [56]. In this study, PGM models with an AUC > 0.95 and a TSS > 0.85 were selected to integrate the model via the weighted average method. The weight of each model is determined by the ratio of its AUC value to the sum of the AUC values of all models. Ultimately, the EM model was applied to identify the PGDs for P. canaliculata under current and future climate scenarios.

2.4. Division of the PGDs and Calculation of Centroids

The suitable habitat of P. canaliculata is classified into four categories using ArcGIS 10.6: unsuitable habitats (0–0.25), low suitability habitats (0.25–0.50), moderate suitability habitats (0.50–0.75), and high suitability habitats (0.75–1.00). Simultaneously, the study compares the total suitable area of P. canaliculata between the current period and the 2021–2040, 2041–2060, 2061–2080 and 2081–2100, evaluating changes in the habitat (increases, decreases, or no change). Calculate the direction and distance of the spatial movement of the centroids based on the changes in the coordinates of the centroids, and then predict the overall changes in the future distribution trend of P. canaliculata. The direction and distance of the spatial movement of the centroids are calculated based on changes in their coordinates, followed by predictions of the overall changes in the future distribution trend of P. canaliculata.

3. Results

3.1. Evaluation of SDMs

Utilizing the global and China-specific distribution data of P. canaliculata, along with environmental variables, the TSS and AUC values of 11 models were obtained after conducting ten runs. XGBOOST (TSS = 1, AUC = 1), RF (TSS = 0.99, AUC = 1), and CTA (TSS = 0.915, AUC = 0.97) exhibited the highest average TSS and AUC values, demonstrating the best accuracy and stability. GAM (TSS = 0.86, AUC = 0.97) and GBM (TSS = 0.88, AUC = 0.98) ranked second in overall performance. However, SER, ANN, FDA, GLM, MAXENT, and MARS exhibited lower TSS and AUC values, indicating poorer accuracy and stability. Models with an AUC value > 0.95 and a TSS value > 0.85 (CTA, GBM, GAM, RF, and XGBOOST) were screened to derive the EM (AUC = 0.99, TSS = 0.96) (Figure 1).

3.2. Evaluation of Environmental Variables

The weights assigned to 9 environmental variables differ across various SDMs. Bio1 and Bio18 have the highest weight proportions and are the most significant environmental variables affecting the PGDs of P. canaliculata, followed by Water Bodies Coverage (WBC), Bio4, and Bio19. Bio7, Bio12, Bio17, and elevation each contribute less than 10 to the PGDs of P. canaliculata. In the XGBOOST, RF, and EM models, Bio1 has the highest weight, with respective values of 32.61, 29.21, and 25.13, while in the GBM, GAM, and CTA models, Bio18 has the highest weight, with values of 33.24, 33.23, and 32.96 respectively. (Figure 2). There is a non-linear relationship between the habitat adaptability of P. canaliculata and various environmental variables. The habitat adaptability of P. canaliculata increases with the rise in Bio1 and Bio18 values. As WBC, Bio4, and Bio19 have minor contributions to the PGDs, the probability of occurrence, as indicated by the response curves of these environmental variables, is relatively stable, leading to a lower impact on the suitable distribution of P. canaliculata. Concurrently, variations are observed across different models in how the changing intervals of environmental variables affect the suitability of P. canaliculata habitats. In models such as EM, CTA, GBM, GAM, XGBOOST, and RF, the critical Bio1 thresholds for P. canaliculata survival are 4.1 °C, 10 °C, 4.6 °C, 6.9 °C, 4.6 °C, and 8.6 °C, respectively (Figure 3). However, across the various models, there is a consensus regarding the critical Bio18 threshold for P. canaliculata, identified as 366 mm.

3.3. Occurrence Records and PGDs of P. canaliculata in China under the Current Climate Conditions

Since the initial discovery of P. canaliculata in China in 1981, its distribution has gradually extended from the southern regions towards the north. From 1981 to 1990, P. canaliculata was primarily found in the coastal areas of Guangdong Province and was introduced to various counties in Hubei, Sichuan, Anhui, Chongqing, Zhejiang, and Yunnan provinces through aquaculture. Between 1991 and 2000, P. canaliculata expanded its range from the coastal to inland areas of Guangdong Province. Concurrently, its presence steadily increased in Guangxi Province. The period from 2001 to 2010 witnessed P. canaliculata’s continued northward expansion, with notable invasions in Hunan, Chongqing, Zhejiang, and Sichuan. This expansion resulted in P. canaliculata occurring in most counties of Zhejiang and Chongqing. From 2011 to 2020, the distribution of P. canaliculata further expanded into Fujian, Hubei, and Sichuan provinces. Additionally, P. canaliculata invasions were recorded in southern regions, including Shanxi, Gansu, and Beijing. Since 2020, P. canaliculata has progressively spread to northern Jiangsu provinces, with invasions also reported in Henan and Shandong provinces (Figure 4A). Under the current climate, high-suitability habitats for P. canaliculata in China are mainly located in Guangdong, Guangxi, Hunan, Hubei, Jiangxi, Anhui, Zhejiang, Fujian, Chongqing, eastern Sichuan, and northern Yunnan provinces, together accounting for 17.66 of the nation’s total area. Moderate suitability habitats for P. canaliculata are mainly found in Yunnan, central Guizhou, central Henan, central Shandong, and parts of Beijing and Hebei, covering 8.73 of China’s total area. Low suitability habitats are primarily situated in select areas of western Yunnan, northern Henan and Hebei, and southern Inner Mongolia, encompassing 6.81 of the nation’s area. Unsuitable habitats are predominant in Inner Mongolia, Liaoning, Gansu, Shanxi, Shaanxi, Qinghai, Xinjiang, and other higher-latitude regions, collectively occupying 66.79 of China’s total land area (Figure 4B).

3.4. PGDs of P. canaliculata under Future Climate Conditions

Compared to its current climate suitability habitats, the total suitable habitats for P. canaliculata show an increasing trend under the four future climate scenarios (SSP1-26, SSP2-45, SSP3-70, SSP5-85) (Figure 5). Under the SSP1-26 climate scenario, there is a notable expansion of the total suitable area for P. canaliculata, particularly in 2021–2040. The proportion of highly suitable habitats for P. canaliculata significantly grows, reaching 28.39%. Regions including northern provinces like Shandong, southern Henan, eastern Hebei, and parts of Beijing emerge as highly suitable habitats for P. canaliculata. Subsequently, from 2041 to 2100, the proportion of highly suitable habitats gradually decreases while moderately suitable habitats increase. In the SSP2-45 climate scenario, the total suitable habitat for P. canaliculata increases, primarily due to the expansion of highly and moderately suitable habitats, with a decrease in low suitable habitats. Both highly and moderately suitable habitats extend northward, reaching the southern part of Liaoning Province. Under SSP3-70 climate conditions, the total suitable habitats for P. canaliculata gradually increase, with an expansion in both highly and moderately suitable areas. However, from 2021 to 2040, there will be a temporary decrease in the total suitable area for P. canaliculata, mainly due to a reduction in low-suitability habitats. Similarly, under SSP5-85 climate conditions, there is a gradual increase in the total, highly, moderately, and low suitable areas for P. canaliculata. The habitat for moderately suitable areas extends to the central region of Heilongjiang (Figure 5 and Figure 6).

3.5. Centroids Shifts of P. canaliculata under Future Climate Change Conditions

The centroid of total suitable habitats for P. canaliculata under various climatic scenarios is depicted in Figure 7. Under current climate conditions, the centroid of suitable habitats is located in the central region of Chongqing. In the SSP1-26 climate scenario, the centroid of suitable habitats initially shifts northward before gradually moving southward over time. Under SSP2-45 climate conditions, the centroid of suitable habitats exhibits oscillatory lateral movements. Under SSP3-70 conditions, the centroid of suitable habitats initially shifts southward, then gradually migrates northward, ultimately settling in Shaanxi province. Lastly, under SSP5-85 conditions, the centroid of suitable habitats steadily migrates northward over time.

4. Discussion

4.1. Model Evaluation

The precision of the PGDs for IAS depends on various factors, such as the availability of species distribution data, the choice of SDMs, and the types of environmental variables. Predictions incorporating both true presence and absence data generally outperform those based only on true presence data in forecasting PGDs [57]. Additionally, MacLeod (2008) [58] suggests that in cases where a species has a broad distribution range, and the true presence data is authentic, its potential distribution range can be theoretically predicted with accuracy, although sampling bias may affect the analysis results. Therefore, this study has employed a method of gathering global distribution data for P. canaliculata and incorporating data of P. canaliculata invasion in China to ensure the accuracy of predicting its PGDs. Due to algorithmic variations among different SDMs, various environmental factors have diverse impacts on species distribution, resulting in uncertainties in predictions from various SDMs. Elith (2006) [59] evaluated the predictive outcomes of 16 models for 226 species, classifying their performance into three categories: Group 1 consists of the most effective models, mainly using machine learning techniques such as MARS community and Boosted Regression Trees (BRT); Group 2 includes models with the second-best predictive performance, primarily based on standard regression methods like Generalized Additive Models (GAM) and Generalized Linear Models (GLM); Conversely, group 3 encompasses models with the lowest predictive performance, primarily using only presence data without inferred absences, such as BIOCLIM, LIVES, and DOMAIN. Similarly, our findings are largely in agreement with those of Elith (2006) [59]. In this study, machine learning methods such as XGBOOST, GTA, and RF were found to have the highest average TSS and AUC values, indicating superior accuracy and stability. Meanwhile, standard regression models such as GAM and GBM ranked second in overall performance. However, models such as SRE (commonly known as BIOCLIM) exhibited lower TSS and AUC values, placing them in the third group. While machine learning models are capable of achieving high predictive accuracy, Breiman (2001) [35], Loh (2011) [34], and Chen (2016) [36] note that this model often lacks explicit statistical principles and is prone to risks like overfitting. The Ensemble model, which combines multiple models, enhances predictive accuracy, especially for rare species with low occurrence rates. This approach effectively mitigates overfitting concerns while maintaining prediction accuracy [60]. Marmion (2009) [61] and Crimmins (2013) [62] used data proportional segmentation for verification evaluation and found that the Ensemble model outperformed individual models, achieving a notable prediction accuracy of 0.97. Predictions from the Ensemble model align with the accuracy of those from machine learning methods. Therefore, this study integrates five different SDMs (XGBOOST, RF, GBM, GAM, and CTA) to enhance the accuracy of the Ensemble model and overcome the limitations of single-model approaches, yielding a more accurate representation of the geographic distribution of P. canaliculata.

4.2. PGDs of P. canaliculata

Although certain regions in Yunnan, Guizhou, and Chongqing are moderately suitable habitats and low suitable habitats for P. canaliculata, it is clear that P. canaliculata can survive in all areas south of the Yangtze River (Figure 4B). Therefore, our study focuses more on the distribution range of P. canaliculata in suitable and unsuitable habitats north of the Yangtze River or at higher latitudes. Under current climate conditions, unsuitable habitats for P. canaliculata are primarily found in Shaanxi, Shanxi, northern Sichuan, Tibet, Qinghai, Xinjiang, and other northern regions. This finding contradicts Yin’s (2022) [47] Maxent predictions, which identified Yunnan, Guizhou, and Chongqing as unsuitable habitats and the western regions of Hunan and Hubei as moderately suitable for P. canaliculata. This discrepancy may stem from Yin’s (2022) [47] reliance solely on global distribution data (GBIF data) and the neglect of population distribution data of P. canaliculata in China. As mentioned in my discussion earlier, predictions of a species’ potential geographic distribution are influenced by data on its actual presence. In Figure 4A, P. canaliculata is found in Chongqing, Yunnan, and Guizhou, damaging rice fields across up to 50,000 hectares [63]. This also explains why our prediction model identifies Chongqing, Yunnan, and Guizhou as suitable areas for P. canaliculata. At the same time, the prediction results align closely with the current occurrence areas of P. canaliculata in China, further validating the accuracy of our integrated model. Furthermore, compared to Zhang (2016) [48] Maxent-based predictions, which indicated that suitable habitats for P. canaliculata were mainly within the Yangtze River Basin and its southern regions, our predictions show a continued northward expansion of suitable habitats for P. canaliculata, the northernmost reaches Jilin Province. This expansion is attributed to the rapid proliferation of P. canaliculata in recent years, particularly in areas of northern Hubei and northern Jiangsu, north of the Yangtze River (Figure 4A). Concurrently, sporadic damage caused by P. canaliculata has been reported in higher latitude areas such as Henan, Shandong, and Beijing within China. The distribution data of P. canaliculata in China significantly influences the prediction of their suitable habitat within the country.

4.3. Environmental Variables

Environmental variables Bio1 and Bio18 have the highest weight, influencing the water and heat requirements crucial for the northward expansion of P. canaliculata. These findings are consistent with the research of Lv (2011) [64], Yang (2022) [26], Yin (2023) [47], Seuffert (2021, 2024) [28,29], and Qin (2023) [65], all highlighting the primary environmental variables—water and heat—as influencers of the suitable habitat for P. canaliculata. P. canaliculata predominantly thrives in warm and humid environments, while dry conditions and cold winters act as limiting factors for its range expansion. In the EM, CTA GBM, GAM, XGBOOS, and RF, the survival thresholds of P. canaliculata on annual mean temperature (Bio1) are 4.1 °C, 10 °C, 4.6 °C, 6.9 °C, 4.6 °C, and 8.6 °C, respectively. This means that different models lead to differing estimates of suitable habitat areas for P. canaliculata in China. For example, the CTA model, which has the highest Bio1 threshold (10 °C), predicts a smaller suitable habitat area for P. canaliculata compared to the EM, GBM, GAM, and XGBOOST models. The tolerance of P. canaliculata to extremely low-temperature conditions in the invaded areas is a critical factor in determining whether it can successfully establish new populations there. Yoshida (2009) [66] has demonstrated that although P. canaliculata eggs cannot overwinter in extremely low-temperature conditions, they can endure cold winters and resume reproduction when temperatures rise in the following year. Furthermore, P. canaliculata can increase its energy supply and tolerance capacity, thereby enhancing its cold resistance [67]. For example, Wada and Matsukura (2007) [27] discovered that P. canaliculata in Ibaraki Prefecture, Japan, not only survived cold winters but also exhibited greater cold tolerance compared to P. canaliculata in southern Japan. This phenomenon helps explain why, despite P. canaliculata being introduced as food in Tianmen City, Hubei, and Nanjing City, Jiangsu, in 1988, a large-scale outbreak occurred in Hubei only after 2015 and in the Jiangsu region after 2020, potentially due to the gradual increase in P. canaliculata cold tolerance.

4.4. PGDs under Climate Change

Environmental variables Bio1 and Bio18 play a crucial role in influencing the geographic distribution of P. canaliculata. Consequently, changes in Bio1 and Bio18 under different future climate scenarios will dictate the PGDs of this species. To elucidate the relationship between changes in P. canaliculata’s suitable habitats and Bio1 and Bio18, we visualized raster data of Bio1 and Bio18 under current and future climates and marked the lowest thresholds of Bio1 (4.1 °C) and Bio18 (366 mm) for P. canaliculata’s survival, as well as the optimal thresholds of Bio1 (16.7 °C) and Bio18 (897 mm) (Figures S1–S3). Under SSP1-26 climate conditions, the 4.1 °C contour does not exhibit significant northward movement. However, under the SSP2-45, SSP3-70, and SSP5-85 climate scenarios, both the 4.1 °C and 16.7 °C contours exhibit significant northward movement over time (Figure S2). Similarly, under SSP1-26 and SSP2-45 climate conditions, the 366 mm contour has not extended northward. However, under SSP3-70 and SSP5-85 climate scenarios, the 366 mm contour for P. canaliculata survival has moved northward over time (Figure S3). These observations suggest that due to climate change, currently inhospitable dry and cold regions may become warm and humid, providing suitable environments for P. canaliculata’s survival. In our study, under SSP1-26 and SSP2-45 climate scenarios, the total suitable habitats for P. canaliculata remain stable. In contrast, under SSP3-70 and SSP5-85 climate scenarios, the total suitable habitats for P. canaliculata increase. The northward expansion of the PGDs for P. canaliculata under various climate conditions aligns closely with the northward shift of the 4.1 °C and 366 mm contours. Therefore, limiting high fossil fuel consumption not only protects the environment but also helps curtail the colonization and spread of invasive P. canaliculata.

5. Conclusions

In this study, we reconstructed the invasive history and analyzed the PGDs of P. canaliculata in China. We selected five SDMs—CTA, GBM, GAM, RF, and XGBOOST—to develop an ensemble model (EM) that predicts suitable habitats for P. canaliculata under current and future climate scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5). The EM, exhibiting an AUC of 0.99 and a TSS of 0.96, produced more precise predictions compared to the individual models. Under current climate conditions, the highly suitable habitats for P. canaliculata are primarily located south of the Yangtze River. The 4.1 °C contour of Bio1 and the 366 mm contour of Bio18 determine the northernmost geographical distribution of P. canaliculata. Sufficient water and heat conditions provide a favorable living environment for P. canaliculata. Under future climate scenarios, the total suitable habitats for P. canaliculata are expected to expand into higher latitude regions, particularly under SSP3-7.0 and SSP5-8.5 conditions. Future efforts should include field surveys in high-risk areas, such as Shandong, Henan, and Hebei provinces, and increased focus on preventing the further spread of P. canaliculata in northern China.