Tropical mountains are renowned for their high biodiversity both on land and in water. However, the increasing threat of invasive species, exacerbated by human activities and climate change, poses a major threat to this rich diversity (He et al., 2023), particularly at higher elevations (Ruiz et al., 2008). The growth of human populations at high altitudes has caused increased habitat fragmentation (He et al., 2023), enabling the establishment and expansion of invasive species (Fuentes-Lillo et al., 2023). Additionally, elevated temperatures facilitate the upward propagation of both native and invasive species (Iseli et al., 2023), which subsequently intensifies competition at the summit. Although these processes could lead to the degradation of native biodiversity, a comprehensive assessment of the potential threats posed by the climate-driven expansion of invasive species is currently lacking in tropical mountainous regions, such as the northern Andes.

The spread of invasive species can negatively impact ecosystems and human well-being. In tropical mountains, impacts range from displacement of native and cultivated plant species to soil degradation by altering natural water and nutrient cycles (Cárdenas López et al., 2017; Costa et al., 2019; Le Maitre et al., 2011; Valdez-Ramírez et al., 2020). Currently, invasive species such as Genista monspessulana ((L.) L.A.S. Johnson, 1962), Melinis minutiflora (P. Beauv., 1812), and Ulex europaeus (L., 1753) also increase the risk of fire in rural and urban landscapes by increasing fuel loads and fire corridors (Anderson and Anderson, 2009; D’Antonio et al., 2011; Pauchard et al., 2008), whereas others, such as Ricinus communis (L., 1753), corrode the integrity of roads and homes (Sandoval et al., 2022). Combined, these impacts could cost as much as USD 10 billion to governments in Andean countries (Heringer et al., 2021) and cause incalculable losses in terms of Nature Contributions to People (Roy et al., 2024). Thus, the future of biodiversity in the high mountains and the well-being of the people living in these areas (more than 150 people per square kilometer, Williams, 2011) are inextricably linked to the effective containment of invasive species.

Climate change could increase the invasion risk in high-elevation ecosystems such as the páramos above the Andean tree line. Projections indicate that changes in climate patterns will cause significant shifts in the distribution of native biodiversity at higher elevations (Peyre, 2022). Additionally, climate change may weaken historical environmental constraints, making it easier for non-native plants to establish themselves in the region (Valencia et al., 2013). Although both processes may create 'open spaces' for invasive species, it is currently unknown which species will use these newly accessible habitats and to what extent Páramo may be invaded. To address these knowledge gaps, we utilized species distribution modeling to identify the present-day and future zones of higher suitability for the establishment of the most harmful invasive plants currently found in the high Andean Mountains. In particular, we asked: (1) How extensive is the range of conditions currently suitable for invasive species in the tropical Andes? (2) How will range change under different climate change scenarios? (3) How much of the Páramo ecosystem is likely suitable for the establishment of invasive species? We investigated these questions in relation to three global circulation models and three climate change scenarios, each with low, mid, and high levels of carbon emission. As such, we aimed to capture the uncertainties in habitat suitability projections that may arise from different climate change projections.

Our modelling approach consists of four stages (Fig. 1): (A) preparing data of species occurrences and environmental predictors, (B) training three algorithms to predict the distribution of each species in its native region, (C) validating the usefulness of the models to project the distribution of invasive species in the tropical Andes, and (D) employing the trained models to project the potential distribution of invasive species for the mid-century. All modelling stages were conducted using R ver. 4.2.1 and implemented in accordance with a standardized and replicable methodological pipeline, which has been successfully employed with other biological groups in situations of limited data (e.g., Gouvêa et al., 2024). R scripts to reproduce this pipeline are available in a permanent repository on GitHub (https://github.com/jorgeassis/speciesDistributionModelling).

Fig. 1.

Flowchart showing the step-by-step procedure for projecting the current and future distribution of invasive plant species in the tropical Andes. (A) Data preparation included (1) the compilation and curation of species occurrences in their native regions; (2) the implementation of a Principal Component Analysis (PCA) performed on the environmental variables in the native region; (3) the application of a kernel-based filter, which splits the environmental space into two subspaces associated with either the environmental conditions more suitable for the species (in yellow) or those associated with less or unsuitable environmental conditions; and (4) uniform sampling of pseudo-across a sampling grid of a chosen resolution superimposed on the 2-dimensional environmental space. (B) Model training and evaluation included calibration of three algorithms and evaluation of their individual and ensemble performance using a set of metrics. (C) Model validation included evaluation of model transferability (i.e., performance outside the training region) using occurrences registered in the tropical Andes. (D) Ensemble prediction involved using the ensemble of three algorithms to project the distribution of each invasive plant species under three scenarios of future climate change.

(0.58MB).

Below, we provide a detailed description of the steps and data involved in each stage of the modelling approach.

• A

  Data preparation

Occurrence dataset. Occurrence data for eleven plant invasive species (Table 1) in the Tropical Andes were compiled from various databases, including the Global Biodiversity Information Facility (GBIF, retrieved in November 2023), the Catalog of Plants of Colombia, Tropicos (www.tropicos.org), and scientific publications (e.g., Sandoval et al., 2022). To ensure data quality, we retained only records from the Tropical Andes and regions of origin for each species (Table 1) with taxonomic verification and GPS errors below 1 km. Furthermore, the retained records were uniformly sampled within a 2-dimensional environmental space created using principal component analysis (PCA) with all environmental predictors (Fig. 1A). The 'uniformSampling' function from the USE R package (Da Re et al., 2023) was used for this purpose. This environmentally driven thinning approach addresses the biases that stem from uneven sampling across the entire environmental range, as highlighted by Castellanos et al., 2019 and Da Re et al., 2023. By employing PCA to capture over 75% of the total environmental variation and thinning occurrences within this space, we were able to achieve a balanced representation across the environmental gradient, thereby reducing the oversampling of accessible areas and undersampling of remote regions. As such, this approach enhances the reliability of our distribution models, particularly regarding the projection of species distributions in future climate change scenarios (Varela et al., 2014). The final number of records used for calibrating and evaluating species distribution models is described in Table S1.

Pseudo-absence dataset. Given that the available species distribution data contain presence-only records, we adopted the approach of Da Re and colleagues (2023) to generate a sample of pseudo-absences. This involved three steps: recovering the core area of the species' bioclimatic niche within the environmental space using records within species native ranges; using a kernel-based filter to split the environmental space into subspaces associated with suitable and less suitable conditions; and uniformly sampling pseudo-absences across subspaces representing less suitable conditions. All three steps were implemented using the 'paSampling' function from the USE R package (Da Re et al., 2023). The native ranges of the species were delineated based on the extent of the phylogeographic region in which the evolutionary lineage of each species originated. This approach allows for the consideration of all potentially accessible areas from the lineage's center of origin and the species' dispersal potential and associated barriers, as proposed by Rojas‐Soto et al. (2024). The boundaries of each phylogeographic region were extracted from a study on floristic regions of the world by Liu et al. (2023).

Environmental Predictors. Fifteen environmental predictors were chosen based on established ecological relationships (see Table S1). These predictors were derived from various sources and encompassed climate (Brun et al., 2022), soil characteristics (Lembrechts et al., 2022; Poggio et al., 2021), and human-mediated effects (e.g., livestock density, Gilbert et al., 2018). All environmental layers were averaged to a grid scale of 1 × 1 km grid cells and reprojected to the ESPG:4326 coordinate system using the functions ‘resample’ and ‘project’ from the terra R package (Hijmans, 2024). The Tropical Andes Mountain range, spanning from central Peru to northern Colombia and Venezuela, was utilized for the projections. The Huancabamba depression, located at 6 °S–79 °W, serves as the southern limit as it denotes the transition between tropical and subtropical climates (Weigend, 2002).

• B

  Model training and evaluation

For each species under consideration, we employed three distinct machine learning algorithms to develop species distribution models (SDMs) (Fig. 1B). The algorithms utilized were Adaptive Boosting (AdaBoost) (Krause-Jensen et al., 2020), Boosted Regression Trees (BRT) (Elith et al., 2008), and Extreme Gradient Boosting (XGBoost) (Chen and Guestrin, 2016). These algorithms are recognized for their effectiveness in the presence of limited data (Barbet‐Massin et al., 2012) as well as their capacity to capture the intricate interplay between predictors and response variables and their adaptability to altered hyperparameters (Elith et al., 2008). Each algorithm has a distinct advantage. AdaBoost is an effective method for handling imbalanced datasets, BRT can capture intricate species-environment interactions, and XGBoost enhances computational efficiency while reducing overfitting through regularization. Accordingly, within an ensemble forecasting framework (Araujo and New, 2007), these models utilize their complementary strengths, thereby enhancing their transferability and facilitating reliable species distribution predictions within and outside their native regions (Elith et al., 2008). AdaBoost models were implemented using the 'mboost' function from the mboost package (Hofner et al., 2014); BRT models were implemented using the 'gbm' function from the gbm package (Greenwell et al., 2019); and XGBoost models were implemented using the 'xgboost' function from the xgboost package (Chen and Guestrin, 2016).

Following Gouvêa et al. (2024), we used a 10-fold cross-validation framework and the grid search method to identify the optimal combination of hyperparameter values for the three algorithms. The grid for BRT changes with respect to the learning rate (0.1, 0.01, and 0.001), tree complexity (1–4), and the number of trees (50–1000, step 50). The grid for AdaBoost changes with respect to the number of interactions (50–250, step 50), degrees of freedom (1–12), and shrinkage (0.25–1, step 0.25). The grid for XGBoost changes with respect to gamma (0–5, step 1), interaction depth (1–4), shrinkage (0.1−0.5, step 0.1), and the number of rounds (10–100, step 10). We also ensured that the models had monotonic responses to positive or negative constraints (Table S2). This helped to reduce overfitting and align with the expected relationships between species occurrences and environmental predictors (Hofner et al., 2011). Prior to modeling, we examined collinearity among predictors using the variance inflation factor (VIF) metric and selected a reduced set of uncorrelated predictors through a stepwise procedure to avoid multicollinearity (Dormann et al., 2013). The threshold for selecting uncorrelated predictors was set at a VIF < 3.

The performance of the individual models and their weighted average ensemble (Araújo et al., 2006) in predicting the distribution of each species within each species' native range (see Table 1) was evaluated using the Boyce index, area under the receiver operating characteristic curve (AUC), and total score (TSS) (Boyce et al., 2002; Fielding and Bell, 1997). The values of these metrics range from −1 to +1 for the Boyce index and from 0 to 1 for the AUC and TSS. A positive Boyce index above 0, or AUC and TSS values above 0.5, indicate that the model's predictions outperform random expectations, with values approaching 1 indicating strong agreement between the predicted and observed distribution of each species (Hirzel et al., 2006).

• C

  Model validation

As the models were trained (and evaluated) with data from regions where the species are native, it is essential to ascertain their transferability to regions where the species are considered invasive (Fig. 1C). In other words, the models should be validated by determining whether the projected distribution in the Andes chain aligns with documented records. We validated the transferability of the ensemble model using the Boyce index (Boyce et al., 2002). One advantage of the Boyce index is that it is not dependent on threshold-based metrics such as TSS and AUC. The Boyce index assesses the performance of a model across the entire suitability gradient, thereby obviating the need for a fixed cutoff. This makes it particularly well suited for situations in which reliable absence data are lacking (Hirzel et al., 2006). Moreover, because it captures the relationship between predicted suitability and observed occurrences across all suitability levels, the Boyce index provides a more comprehensive assessment of model performance and reliability (Leroy et al., 2018).

The presence data of each species were used to validate the transferability of the model to the tropical Andes. As with the occurrence dataset (see above), occurrence records across the Tropical Andes were uniformly sampled using the 'uniformsampling' function of the USE R package (Da Re et al., 2023). This was performed to prevent the inclusion of an excessive number of occurrences derived from sites with similar environmental characteristics, which would otherwise inflate metric values. The numbers of occurrences used to validate the distribution models for each species are presented in Table S1.

• D

  Ensemble forecasting

For each species, we combined the outcome of the different models within an ensemble forecasting framework (Araujo and New, 2007) to produce maps of future (2040–2070) habitat suitability (Fig. 1D). To assess the uncertainties in projections of habitat suitability that may result from various scenarios of future climate change, we generated projections based on three shared socio-economic pathways (SSP) scenarios: SSP126, SSP370, and SSP585. These scenarios correspond to the lowest (best), mid, and highest (worst) greenhouse gas emission scenarios (Riahi et al., 2017). Each SSP scenario was produced using three different global circulation models (GFDL-ESM [Dunne et al., 2020], MPI-ESM [Gutjahr et al., 2019], IPSL-CM6A [Boucher et al., 2020]). Distribution maps are available in the Zenodo repository (DOI: 10.5281/zenodo.10907376).

Habitat suitability maps for the present and future were transformed into binary maps using a threshold that optimizes TSS (sensitivity and specificity) (Liu et al., 2016). To measure the distribution area of each species, we summed all cells where its presence was predicted. Similarly, we determined the elevational ranges of each species based on the altitude of the cells in which they were predicted, obtained from a digital elevation model (DEM). To determine whether the species' elevation range will increase or decrease in the future, we used Cliff's delta. This effect size measures the probability that a selected value from one group will be greater than a selected value from the other group, ranging from −1 to +1, with values closer to zero indicating no overlap between the two groups. Positive values indicated an increase in the elevational range of the species, while negative values indicated a decrease. The extent of elevational change was determined using the thresholds outlined in Romano et al. (2006): D < 0.147 is considered 'negligible,' D < 0.33 is considered 'small,' D < 0.474 is considered 'medium,' and anything beyond that is considered 'large’.

Tropical mountains provide a favorable environment for the potential establishment of invasive plant species. Currently, 18–60% of the region is estimated to have suitable conditions for at least one of the studied species. Our models indicated that suitability is higher in areas with higher dry season rainfall, stable temperatures throughout the year, specific soil properties (such as high pH and organic carbon), and/or dense livestock populations. However, as the climate changes, the spatial distribution of these conditions may change. By the mid-21st century, under the most benign scenario, invasive species could lose between 2,000 and 28,000 km2 of suitable environmental areas and gain between 1,000 and 10,000 km2 of more suitable areas. Although this result may seem counterintuitive, it is consistent with the expected effects of rising temperature on mountain plant species (Dirnböck et al., 2011). Specifically, rising temperatures are expected to affect the probability of species to establish new populations at lower elevations, particularly for cold-adapted species, resulting in upward shifts in their range (Alexander et al., 2018). Such upward shifts and range contractions have already been documented for U. europaeus in Colombia (Ángel-Vallejo et al., 2024), with a limited seed germination range of 0–15 °C (Udo et al., 2016) identified as the limiting factor, confining the species to the cooler highlands of the tropical Andes.

Our projections can be used to guide future conservation and management plans. Areas predicted to be highly suitable for invasive species could be targeted for conservation or management actions, such as population removal, blockage of invasion pathways, and active/passive restoration of ecosystems at a high risk of invasion (Roy et al., 2024). The implementation of such measures may prove pivotal in areas exhibiting a higher degree of suitability for more than one species (illustrated by the red cells in Figs. 4b and c), and notably in areas encompassing páramo ecosystems. The establishment and propagation of invasive species into páramos has the potential to pose a substantial risk to the region's water regulation services. Grasses such as A. donax and M. minutiflora have the potential to alter the soil structure and hinder the regeneration of native species (Zenni et al., 2022), which could ultimately impact the system's ability to retain water. Furthermore, the establishment of species, such as G. monspessulana and U. europaeus, combined with the projected longer droughts in the area (González‐Trujillo et al., 2024), may increase soil impoverishment and fire hazards (Galappaththi et al., 2023), ultimately compromising the paramos' capacity for water regulation and provisioning.

It is important to emphasize that given the life histories of the eight species studied, both modeling and management efforts should be species-specific and integrate multiple strategies. The suitability of the habitat for each species is determined not only by distinct abiotic factors (Supp. Mat. 1) but also by varying human activities that contribute to the establishment and spread of these species (Roy et al., 2024). Consequently, effective invasion control must go beyond managing and protecting the most suitable areas; it must also address key factors, such as invasion pathways. For several species, human-facilitated spread is driven by economic and cultural factors. For example, R. communis and M. minutiflora have been introduced for agricultural purposes, while T. Alata, A. donax, and H. coronarium have been spread for ornamental use. Thus, managing invasion pathways alongside the designation of priority protection areas is critical for effective control.

A more comprehensive approach to invasion control can be found in the restoration of native ecosystems. Restoration can function as a filter, impeding the establishment of invasive species, while facilitating the colonization of native species (Bakker and Wilson, 2004). Besides, local efforts to restore high Andean forests and páramos have shown promising results in controlling the spread of exotic grasses, gorses, and shrubs when combined with proactive management actions, such as the mechanical removal of species, prescribed burns, and artificial shading (Castiblanco-Álvarez, 2012; Gómez-Ruiz et al., 2013). The challenge ahead is to raise awareness of the pressing issue among stakeholders to increase the resources available for testing and scaling up the most feasible combinations of restoration and active management efforts.

Although our projections identify potential management priority areas until the mid-century, it is important to acknowledge the inherent limitations of our modelling approach. First, it should be noted that our models are conservative in nature, and that habitat suitability is predicted based on the assumption of no niche shift, which may not be applicable to all species (Atwater et al., 2017). For example, recent evidence indicates that U. europaeus has expanded its climate niche and colonized areas with climates distinct from its native range (Christina et al., 2020). This indicates that areas deemed suitable for invasion may be more expansive than our estimates suggest. Second, despite the careful selection of predictors, the analysis did not include important factors such as groundwater levels, spatial distribution of man-made corridors, and incidence of biotic interactions. The inclusion of data on these predictors could enhance the precision of our projections; however, they are currently unavailable at a global level or are not projected for future climate scenarios. For instance, the distribution of hydrophyte species such as A. donax is constrained to areas with shallow groundwater (Quinn and Holt, 2008), whereas species such as R. communis are primarily dispersed along man-made corridors such as roads and urban areas (Sandoval et al., 2022).

In the absence of more detailed knowledge, our projections provide the first comprehensive assessment of the suitability range of invasive species in the tropical Andes under different climatic scenarios. This information provides an overview of the magnitude of the problem and can help prioritize areas for implementing ecosystem-based management actions, such as restoration. In addition, these findings can inform future theoretical and empirical research on the environmental factors contributing to the spread of invasive species in the tropics, which has historically received less attention when studying biological invasions (Chiu et al., 2023).

JDGT conceived the study. Species data collection and analyses were performed by JDGT. The first draft of the manuscript was written by JDGT, and all authors commented on subsequent versions of the manuscript. All authors read and approved the final manuscript.

The authors have no relevant financial or non-financial interests to disclose.