The study area encompassed the tropicalpine páramo region in the northern Andes, ranging from the depression Huancabamba in northern Peru, often described as a main biogeographic barrier between the northern and central Andes, to the mountains of Ecuador, Colombia and Venezuela (∼Lat. 11°N-5°S; Long: 65°–82°W; Luteyn et al., 1999; Weigend, 2002). Small páramo areas out the main Andean cordilleras, such as Amazonian volcanoes or central American páramos were excluded (Kappelle and Horn, 2016; Lozano et al., 2020). To predict species’ distributions in the future and account for potential downward migrations, a small elevation buffer was added to the actual ∼3,500 m treeline (Peyre et al., 2021) and therefore the 3,000 m isohypse was used to delimit the study area (Figure 1). No upper páramo limit was considered since the natural páramo frontier is represented by mountain summits or glaciers, which currently encompass about 1% of the northern Andes and will undergo further reductions by 2070 (Dussaillant et al., 2019; Peyre et al., 2021).

The substantial species dataset used was identical to Peyre et al. (in review) to ensure comparability of the predictive results for future conditions with those for the LGM. All analyses were performed in R 4.0.2 (R Core Team, 2020). The species list was retrieved from the VegParamo – the flora and vegetation database for the Andean páramo (Peyre et al., 2015)¹ and included 664 vascular plant species spread over the main páramo plant families. The VegParamo vegetation plots were then downloaded and complemented for rare species (<50 occurences) with additional occurrence records (timeperiod 1980–2019) obtained from 21 open-access herbarium databases. Then, the data cleaning procedure included: (i) merging sub-species at the species level; (ii) removing non-vascular species and taxa identified at the genus level or above; (iii) removing plots and occurrence points either poorly georeferenced (>1 km) or outside the study area; and (iv) checking and removing observation duplicates between vegetation plots and species occurrences by applying a 500 m radius buffer to all VegParamo plot coordinates (data source and year or observation were final selection criteria). To prevent overfitting issues in the subsequent statistical models, very rare species with less than 10 occurences in the cleaned dataset were removed (Merow et al., 2014). This threshold was lower than the usually recommended minimum 15–25 observations to reduce model uncertainty and was defined to account for the overall sampling gap in the Andes (van Proosdij et al., 2016; Santini et al., 2021). As a result, the final dataset consisted of 3,439 georeferrenced presence-absence entries for the 664 species: (i) 1,986 observations from VegParamo with presence and absence information registered in vegetation plots; and (ii) 1,453 additional occurrences from open-access herbarium databases (Peyre, 2022).

Endemics -- defined as those species almost exclusively present in the tropicalpine belt of the northern Andes -- were then drawn from the 664 species list. This list of páramo endemic species is the first attempt to systematically assess endemicity in the region and was achieved through a two-step procedure. For each species, all its occurrence points were downloaded from GBIF² between 1980 and 2021 to fit the study’s timescale. Duplicate observations as well as obvious poorly georeferenced points, i.e., in habitats unfit of the species’ ecology, and data from collections, i.e., on other continents and in major cities, were discarded. To classify as endemic, a species first needed to present at least 90% of its occurrence points in the longitude-latitude spatial frame defined by the study area (a 10% of outlier points was allowed). Second, at least 90% of its datapoints had to be distributed above the elevation mask set at the 2,500 m isohypse, to include points with unprecise georeferencing (based on the 1 km UTM scale). As a result, 247 of the 664 species in the dataset constituted the endemic species pool (Supplementary Table 1).

Because the following models required a predictor able to emulate species’ dispersal capacity (DP), I retrieved the maximum dispersal distance (DD), and associated dispersal factor (DF) calculated for the 664 species in Peyre et al. (in review). DD were previously computed using the disperSal function for the R environment (Tamme et al., 2014) that relies on trait data (plant height, seed mass, and dispersal mode) and predicts kernel-shaped dispersal at a given distance (in km/year). Because dispersal frequency remains largely unknown for this species pool, periodic annual dispersal was assumed for all species for simplicity reasons. The DF was calculated by rescaling the DD onto a 0–1 scale, which allowed to rank species from worst to best disperser. Morever, I extrapolated the species’ dispersal range to the study area in the form of a raster map for current conditions, or dispersal predictor – DP. This procedure meant applying the iForce function to the study area, creating a buffer of currently accessible raster cells around the presence points, with a λ parameter defined as the DF (iSDM package, Hattab et al., 2017; Peyre et al., 2020).

Similarly to the previous step, the bioclimatic variables used in this study followed Peyre et al. (in review) to allow comparibility between future predictions (2070) and predicted climatic refugia since the LGM. The predictors set included (i) bio1 – mean annual temperature; (ii) bio2 – mean diurnal temperature range; (iii) bio4 – temperature seasonality; (iv) bio12 – total annual precipitation; (v) bio15 – precipitation seasonality; (vi) bio17 – precipitation of the driest quarter; and (vii) bio18 – precipitation of the warmest quarter. These variables were previously checked for multicollinearity using a variance inflation factor analysis (VIF < 10) and considered fit for the study area (vif function, usdm package; Naimi et al., 2014). Selecting the same variable set for all species may not fit the data as well as species-specific predictors, however, this methodological choice was made based on the following points: (i) little is known about the ecological correlations between each species’ distribution and environmental factors (Smith and Santos, 2020) and these correlations may be significantly affected by the large spatial study scale; (ii) keeping variable amounts relatively low helps maintain intermediate parameterization complexity in models (Brun et al., 2020); and (iii) the selected variables have been deemed important to explain current plant species richness patterns in the páramo (Peyre et al., 2019).

The seven bioclimatic variables were downloaded at a 30 arcsec resolution (∼1 km at the equator) from the CHELSA project v. 1.1³ (Karger et al., 2017) as mean values computed for the 1979–2013 period. The variable set was considered adequate to represent current conditions (set at 2010 to enable the tracking of species over time) and fitted the species distribution data spread over the 1980–2019 timeframe. The same bioclimatic predictors were downloaded for 2050 (mean values over 2040–2060) and 2070 (mean values over 2050–2070). My intent was to assess and compare the effects of plausible climate change scenarios on the páramo flora, and therefore two Representative Concentration Pathways (RCP), CMIP5-RCP45 (moderate scenario, hereafter referred to as RCP45), and CMIP5-RCP85 (severe scenario, hereafter referred to as RCP85), were chosen. Furthermore, to encompass sufficient climatic variance and reduce uncertainty in future predictions (Thuiller et al., 2019), 5 General Circulation Models (GCM) were selected per RCP: CNRM-CM5 (CNRM); FGOALS-g2 (FGOALS); IPSL-CMA-LR (IPSL); MPI-ESM-P (MPI); and MRI-CGCM (MRI). The 5 GCMs were selected on the bases of (i) their dissimilarity in the GCM correlation matrix presented in Sanderson et al. (2015), and (ii) because they matched the scenarios used to represent LGM climatic conditions with PMIP3 data (Peyre et al., in review).

In order to follow the progressive shift in species’ distribution between today and 2070, the future bioclimatic data were interpolated per decade so to apply step-by-step dispersal constraints. Simple linear regressions were used on the future variables cropped to the study area and regression coefficients were extracted per variable (one for the 2000–2050 period and one for 2050–2070 period). Then, the previously computed coefficient was successively applied to interpolate each variable for the years 2020, 2030, 2040, and 2060, while values for the years 2010, 2050, and 2070 were directly extracted from the original downloaded variables (Adhikari et al., 2018).

Species distribution models allow to correlate species observations (presence or presence-absence) to environmental data and extrapolate these relationships onto different spatio-temporal contexts (Araújo et al., 2019). I followed the procedure detailed in Peyre et al. (2020) that (i) uses SDMs to establish species distributions for current conditions and project them into the future per decade until 2070 for the 10 climate change scenarios; and (ii) successively applies dispersal limitations to the decadal projections.

First, to represent the species’ current (realized) distribution, SDMs were fitted with eight predictors: the seven bioclimatic variables and the previously computed species-specific predictor – DP, emulating the species’ actual dispersal range over the study area. The species data was retrieved from the final 3,439 observations dataset, from which were extracted both presence points (VegParamo data and curated additional occurrences) and absence points (VegParamo data). In order to enhance the models’ fit and predictive success, four algorithms across different families were included to perform the SDMs: Generalized Linear Model, Random Forest, Multivariate Adaptive Regression Splines, and Artificial Neural Network. For a given species, the modeling procedure consisted of (i) partitioning the dataset into training and testing sets through random splitting of the data in 75/25 proportions (Araújo et al., 2019); (ii) feeding the training set to an algorithm; and (iii) validating the resulting model with the testing set (biomod2 package; Thuiller et al., 2016). Out of the 100 models per species (four algorithms, 25 runs each), only those with a performance of 0.6 and above on the true skill statistic (TSS) metric were selected for the final ensembling of their predictive outputs (Allouche et al., 2006; Fernandes et al., 2019). Ensemble modeling was performed to better visualize the true correlation signal from the statistical noise intrinsic to each algorithm (Dormann et al., 2018; Hao et al., 2019) and then projected onto the study area to forecast species’ probabilistic distribution at present. Then, these predictions were binarized by applying a species-specific threshold that maximizes the sum of the sensitivity and specificity metrics in the retained models (optimal.thresholds function; PresenceAbsence package; Liu et al., 2005; Freeman and Moisen, 2008).

Second, to represent the future potential distribution of each species, SDMs were fitted with the seven bioclimatic variables only. An additional filter was applied to the species’ absence points to remove those absences that not necessarily represent environmental filtering of the species but instead dispersal limitations. This was achieved by applying the pDLA function, with a α parameter set at 0.5 that considered a 50% probability to classify an absence point as a dispersal-limitation-absence (iSDM package; Hattab et al., 2017). The rest of the modeling procedure remained identical, and every model ensemble was projected onto the variable set at each decade (2020, 2030, 2040, 2050, 2060, and 2070) for each RCP and GCM, and finally binarized.

Third for each species, its current realized distribution and future potential distributions as binary maps were ordered temporally and constrained by dispersal at each step (per decade) using the MigClim approach (Engler et al., 2012). At each step, the before and after binary predictions was compared and the maximum dispersal distance (DD) was applied to limit the later distribution to the area available and accessible to the species at that time (functions migclim.distance and migclim.plot, migclim package; Engler et al., 2012). In addition, a mask was used to emulate biogeographical barriers below the 3,000 m isohypse, corresponding to Andean forests and human-occupied areas that a species could cross if its DD was greater than the distance separating two suitable raster cells (exceptional long distance dispersal was not considered). The set of decadal binary maps per GCM and RCP offered time series of species redistribution from which can be assessed the species-specific proportion of areas that were colonized (gained), decolonized (lost), available but unoccupied (i.e., potential), and remained occupied (stable) between the present time and 2070. Therefore, these values were calculated and averaged per GCM at the species level (as raster cell numbers and as a percentage in distribution change by 2070), and then summed for the entire species pool and the endemic pool only for each RCP and overall.

The current (2010), future (2070, per RCP and averaged over the 5 GCM), and past (LGM, as a mean scenario) distributions were assembled for each species of the dataset to assess how much of their distribution remained in their climatic refugia between these three different time points. The previously computed SDM for current conditions was used to represent the present, and the LGM distributions were obtained from Peyre et al. (in review). The future distribution was obtained from the final map following the MigClim procedure, therefore accounting for (i) presences – summing total remaining and colonized raster cells; and for (ii) absences – including non-occupied cells, total decolonized cells and potential cells (available but unaccessible in 2070). The maps were stacked and their coinciding (climatic refugium) and non-coincinding pixels (distribution out of the refugia) were counted to express a mean percentage of coherence between present-future (PF) and past-present-future (PPF) timelines. A linear regression was finally performed for each timeline and altogether to reflect any relationship between the size of the climatic refugium per species (in Log form) and the distribution percentage occuring out of refugia, accross RCPs.

In terms of biodiversity loss, my results identified a high extirpation risk for 55 species, including 22 endemics by 2070. Many of these species are common from the sub-páramo transition with montane forests (∼3,000–3,500 m) and known in few locations, for example Brachyotum confertum in Ecuador, Monticalia abietina and Pentacalia pulchella in Colombia, or Diplostephium obtusum and Ternstroemia meridionalis in Venezuela. These species should be primary candidates for conservation measures, since they often are also threatened by anthropogenic activities and habitat degradation near the treeline (Hofstede et al., 2003; Peyre et al., 2018). The future upheaval of the treeline due to increased thermophilization of montane forests (Fadrique et al., 2018) might in addition procure additional competition with rising tree species which might cause a further reduction of the shrub species’ distribution in the sub-páramo, but also understorey species such as the endemic Epidendrum frutex and Geranium velutinum.

Models established that páramos species will tend to increase their actual distribution by 15.52% in 2070, despite showing important species-specific variations. This finding contradicts my first hypothesis H1, based on Valencia et al. (2020) for the Espeletiinae group and Peyre et al. (2020) for plant communities in the Ecuadorian super-páramo, both of which suggested overall reductions in species’ distributions. Our contrasted result is most probably due to the substantial diversity of life-forms and strategies represented within the 664 species pool. In fact, Valencia et al. (2020) and Peyre et al. (2020) focused on plant species with poor-dispersing strategies, namely the unassisted anemochorous Espeletiinae and mostly small plants adapted to severe conditions, hence investing in survival more than growth and reproduction. In contrast, the species pool modeled here also included many species representative of (i) the sub-páramo transition with Andean forests dominated by dwarf forests and shrublands, and (ii) the mid-páramo (approx. 3,500–4,200 m) dominated by mixed grasslands (Peyre et al., 2018). Such species emerging from lower elevations could be advantaged in their colonization endeavors since they often rely either on assisted anemochory (e.g., with wings) or zoochory, which means that they could expand their distribution, provided that they are able to maintain an important part of their current distribution range under environmental changes. The net gain in species distribution across scenarios is attributed to important area gains and losses rather than strong stable areas, which suggests that severe recompositions at the plant community level will occur during the 60 years timespam. We should keep in mind that the models are based on niche requirements only and therefore, they may underestimate species’ persistence in suitable–turned-unsuitable environments, while at the same time slightly overestimate migrations because they do not account for certain constrains that are relevant on site. Indeed, it is probable that species’ (i) acclimation capacity to novel environmental conditions (Leon-Garcia and Lasso, 2019), (ii) species interactions (Llambí et al., 2018; Zimmer et al., 2018), and (iii) establishment and extinction lags, for example due to soil conditions (Anthelme et al., in press; Alexander et al., 2018), will come into play and buffer distribution changes.

My results show that the richness gains on mountain tops already observed in European mountains (Steinbauer et al., 2018) and the Himalayas (Telwala et al., 2013) will likely also take place in the páramo in the next decades, in accordance with hypothesis H2. Models predicted that distribution gains will principally occur at high elevation (central part of cordilleras) while losses will mostly occur at low elevation (edges of the páramo distribution). As previously established, the richness gains expected at high elevations might be alarming, since they are associated to low elevation species rising quickly and competing with local high-elevation species for microrefugia in a more thermophilic setting. Indeed, more than 200 invasive species are present in alpine areas worldwide (Alexander et al., 2015) and such species migrate efficiently thanks to their important seed production and dispersal, as well as their ability to use local nurse plants to aid their progress (Llambí et al., 2018). These species are already taking advantage of higher temperature to establish in un-saturated super-páramo areas and after glacial retreat creates newly available niches (Dussaillant et al., 2019; Rosero et al., 2021) where no-analog temporary plant communities arise (Anthelme et al., in press; Peyre et al., 2020). Nevertheless, endemic species (through not necessarily super-páramo species) are predicted to follow a similar uprising trend, which suggests that the accelerated colonization process will include many native species and, which calls for further vulnerability assessment of the páramo flora with conservation priorization outcomes. One unexpected result was that the northern páramos of the Sierra de Santa Marta and Cordillera de Merida were predicted as important areas for species gains instead of losses that could be attributed to future complex climatic trends in these páramos (Anderson et al., 2011; Helmer et al., 2019; Pabón-Caicedo et al., 2020). In this case, an in-depth analysis of the composition shifts predicted to occur is necessary to establish what species are driving these distribution gains.

In parallel, distribution losses were mostly predicted for low elevations (the edges of the páramo distribution) as well as the eastern slope of the páramo, facing the Amazonia and Orinoquia regions, both for the 664 species pool and endemic pool only. These losses will likely be compensated by the progressive upheaval of montane species at the sub-páramo ecotone in the region (Fadrique et al., 2018) related to regional thermophilization. Regarding the eastern slope specifically, we can expect that the precipitation decrease and important warming of the western Amazonia (Wang et al., 2018) will likely push tree species up the piemount and Andean flank in an accelerated manner, as páramo species will too migrate upslope causing important ecosystem restructuration and changes in functions.

The hability of a species to track future climate change can be partly deduced based on its capacity to occupy its climatic refugium through range-filling since the last great climatic shift on Earth, i.e., the LGM 20,000 years ago (Svenning and Skov, 2004; Dullinger et al., 2012b; Cannone and Pignatti, 2014). It was previously established that vast majority of the 664 species pool showed a climatic refugium between the present and the LGM (84%; Peyre et al., in review), a number that reduced further in the future since only 77% of the species pool potentially had a LGM refugium by 2070 under the novel climates that create new conditions requiring fast responses from species. If a species is shifting its distribution out of its refugium, it likely means that it has good migration capacity. If it remains within its refugium, it means that the species relies on its adaptation capacity to new environments or it is potentially challenged and delaying its extirpation by other factors such as biotic interactions. For example, the endemic Eryngium humboldtii, which despite being predicted with a small future distribution of 10 km² should remain at 95% in its climatic refugium since the LGM. Across all studied species, the results suggested that the chances to find potential refugia are higher in central Ecuador and the Amotape-Huancabamba zone, an area coinciding with important future losses in species distribution according to the models. The refugium pattern encountered here repeats previous findings from a previous study focusing on refugia between the LGM and present time (Peyre et al., in review) and was interpreted as (i) the easier access to low elevation from out-of-the-páramo species, as well as (ii) climatic stability due to low elevations, and (iii) entry points for species colonizing from the central Andes (Weigend, 2002). Nevertheless, the spatial patterns observed here differed from previous findings by dismissing two refugia-rich areas located on the eastern slope of the cordillera and in Central Colombia, according to Peyre et al. (in review). These contrasting results could be, respectively, due to (i) the definition of endemics at the species level instead of the genus level, weighted by the Espeletiinae group diversified in Colombia; and (ii) quick climatic changes predicted to occur between 2010 and 2070 such as the progressive precipitation shifts of the western Amazonia, which could reduce long-term suitable areas for páramo species.

On average, 21.41% of a species’ future distribution in 2070 coincided with a LGM refugium, and an additional 47.59% to a PF refugium, with little distinction for endemics alone. It means that species have a strong tendency to shift out of their climatic refugium over time, and that this shift is supposed to occur quicker in the near future as 30% of a species distribution in 2070 will not coincide with any refugium. I supposed in hypothesis H3 that the smaller the refugium, the higher the probability that a species would climb out of it, however, there was no support for this hypothesis given the lack of PPF correlation, this probably due to the diversity of responses to environmental change on such a broad timescale. However, there was a slight significant positive correlation between refugium size and future distribution for the PF relationship. This finding means that a species with a large distribution in 2070 will have likely remained close to its current distribution, while a species with a small distribution will have likely climbed out of its current refugia.

The models also identified that on average, species will occupy 77.28% or their potential area by 2070, which means that about 20% of their potential distribution will be available but inaccessible due to dispersal limitations at that time. This finding strengthens my hypothesis H1 that dispersal is a significant constraint for páramo species, a pattern that is repeated for endemics who on average occupy 84.28% of their potential area. This insufficient dispersal at play can be related either in the form of (i) short DD, as páramo species disperse on the lower end of the terrestrial angiosperm dispersal spectrum (Peyre et al., 2020; Sklenář et al., 2021) and similarly to other alpine floras (Vittoz et al., 2009; Morgan and Venn, 2017); or (ii) as habitat natural fragmentation impeding dispersal, due to the sky island setting of the northern Andes (Anthelme et al., 2014). In addition, it suggests that narrow-range species with net distribution losses and little potential area could have issues finding niche area in the future, and therefore be at risk of extinction. For instance, Aa argyrolepis and Aa maderoi, both small terrestrial orchids that are predicted to occupy less than 10 km² in the future and experience a net loss of 28.82 and 79.59% of their distribution, respectively, by 2070 also occupy less than 1% of their potential niche due to dispersal limitations.

This study provides a first outlook of páramo plant responses to climate change in the Andean páramo and identified that, even though species distributions might increase overall through uplsope migration and lead to the accelerated colonization process of high mountain areas, there are several concerns. The speed of climate change induces little refugium stability and important losses in species distributions, urging them to rely on their DP to track climate change, which is an issue for many páramo plants. While 9% of the species pool was deemed at risk of extirpation, 10% of the endemics pool could undergo extinction by 2070, which calls for urgent monitoring and conservation measures targeting these species. In addition, I suspect that this proportion would increase accounting for a larger endemics species pool, since many were discarded from theses analyses due to the scarcity of datapoints, hence requiring rare-species model approaches (Breiner et al., 2017). More fitted models to predict these few occurrences species such as Joint SDMs would therefore help predict these species responses to climate change (e.g., Wilkinson et al., 2021). Climate change is already inflincting serious impacts on biodiversity worldwide, but its velocity and local directions in the páramo will likely lead to different levels of exposure and vulnerability for plant species (Dobrowski and Parks, 2016; Panetta et al., 2018). Focusing on the adaptative capacity of species is germane today to complement and fully assess the true vulnerability state of páramo species. In fact, topographically heterogenous landscapes create local microrefugia, which can respond independently or at least slower to environmental change than their surrounding macroclimate (Scherrer and Körner, 2011; Patsiou et al., 2014). This behavior could in turn imply that local páramo plants could be able to persist and show resilience in light of future climate change and delay their extirpation, as observed in several other alpine flora (Frei et al., 2014; Graae et al., 2018). Therefore, future research focusing on climate change responses from páramo plants should incorporate a meso to microscale outlook (Cannone and Pignatti, 2014; e.g., Lenoir et al., 2017). Finally, it would be crucial to account for evolving land-systems in future predictions, since the human response to climate change, in the form of population increase, land-use expansion and intensification, is bound to pose more threats to the páramo biodiversity than climate change itself (Anthelme and Peyre, 2020).