Gwendolyn Peyre
email: gf.peyre@uniandes.edu.co

Jonathan Lenoir

Dirk N Karger

Monica Gomez

Alexander Gonzalez

Antoine Guisan

The fate of páramo plant assemblages in the sky islands of the northern Andes

Keywords: Andes, community assemblages, climate change, maximum dispersal distance, páramo, plant species, Species Distribution Models, Spatially Explicit Species Assemblage Modelling, sky islands

Aims: Assessing climate change impacts on biodiversity is a main scientific challenge, especially in the tropics, therefore, we predicted the future of plant species and communities on the unique páramo sky islands. We implemented the Spatially Explicit Species Assemblage Modelling framework, by i) calculating species' maximum dispersal distance, ii) modelling species distributions at present up to 2100, iii) assembling models into communities. Finally, we assessed the vulnerability of sky islands based on richness and composition changes.

Location: Ecuadorian super-páramo (>4200 m)

Methods: Using species trait data, the maximum dispersal distance of 435 species was calculated. Species distribution models (SDM) were fitted to obtain current and future distribution predictions based on dispersal and bioclimatic factors. The final assemblages for present and 2100 were achieved by stacking all probabilistic SDMs and applying the probability ranking rule. The vulnerability of each sky island was evaluated by quantifying richness and composition changes.

Results: Maximum dispersal distances ranged between 0.008-6027 m/year, and across all scenarios, 70% of models showed a net loss in species distribution while 9% of all species were predicted to undergo extinction by 2100. Local richness was estimated to decrease by 56.63% on average, and composition changes in each sky island suggested a mean loss of 64.74% of their original species pool against a 12.97% gain. Finally, 5% of the sky island floras reconverted from high-elevation to low-elevation species. These numbers were usually more important for high-elevation species and the mountains Pichincha, Ilinizas and Antisana.

Conclusions: Our study is methodologically pioneer and provides novel insight on the future of páramo biodiversity. Significant losses in species distribution and changes in community richness and composition suggest drastic impacts and call for further study considering additional factors, such as land-use. Finally, we recommend focusing monitoring and conservation strategies on the northern sky islands in priority.

Introduction

The climatically-driven redistribution of life on Earth is one of the most challenging environmental threat humanity faces today (Pecl et al. 2017;[START_REF] Diaz | Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services[END_REF]. It is widely accepted that the current climate change is human-induced [START_REF] Team | climate change 2014: synthesis report[END_REF] and that its main impacts on biodiversity include: shifts in phenology [START_REF] Visser | Shifts in phenology due to global climate change: the need for a yardstick[END_REF]; changes in population densities [START_REF] Graae | Stay or go-how topographic complexity influences alpine plant population and community responses to climate change[END_REF], local extinctions [START_REF] Panetta | Climate warming drives local extinction: Evidence from observation and experimentation[END_REF], species range shifts [START_REF] Lenoir | Climate-related range shifts-a global multidimensional synthesis and new research directions[END_REF], and accelerated changes in species richness and assemblages [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF][START_REF] Steinbauer | Accelerated increase in plant species richness on mountain summits is linked to warming[END_REF]. Because biodiversity redistribution implies important changes in ecosystem functioning, human well-being and climate change itself (Pecl et al. 2017), there is an urgent need to increase our knowledge on the complex interplay between climate change and biodiversity and implement integrative scientifically-supported strategies to adapt and mitigate climate change [START_REF] Diaz | Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services[END_REF].

The high-elevation ecosystems in the northern Andes belong to the páramo biogeographical province, a tropical biodiversity hotspot critically threatened by climate change [START_REF] Tovar | Diverging responses of tropical Andean biomes under future climate conditions[END_REF][START_REF] Cuesta | New land in the Neotropics: a review of biotic community, ecosystem, and landscape transformations in the face of climate and glacier change[END_REF]). The páramo is the youngest biodiversity hotspot and the richest tropical alpine province on Earth, totaling 15 different phytogeographical units and around 5000 plant species, of which 60-80% are endemic [START_REF] Sklenář | Island biogeography of tropical alpine floras[END_REF][START_REF] Peyre | Phytoregionalisation of the Andean páramo[END_REF]. It can be altitudinally divided at present into three main elevation belts: i) the sub-páramo (⁓ 3000-3500 m) or lower tree or shrub-dominated ecotone with Andean forests; ii) the mid-páramo (⁓ 3500-4200 m), or páramo proper, including vast grasslands and giant rosette communities; and iii) the super-páramo (> 4200 m), with gradually scattering vegetation [START_REF] Cuatrecasas | Aspectos de la vegetación natural de Colombia[END_REF]. Contrarily to the sub-páramo and mid-páramo, which show certain spatial connectivity and are strongly shaped by human activities, the super-páramo is, hitherto, geographically isolated and above the current agriculture frontier (⁓ 4200m). The super-páramo hence represents the iconic sky islands of the northern Andes and hosts a very rare and endemic biodiversity with little species migration capacity on mountain tops.

Recent estimations predicted a temperature rise of 3ºC (± 1.5) across the páramo by the end of the century, accompanied by a general increase up to + 300 mm in vertical precipitation (rain), however highly variable depending on geography and topography [START_REF] Urrutia | Climate change projections for the tropical Andes using a regional climate model: Temperature and precipitation simulations for the end of the 21st century[END_REF][START_REF] Anderson | Consequences of climate change for ecosystems and ecosystem services in the tropical Andes[END_REF][START_REF] Buytaert | Potential impacts of climate change on the environmental services of humid tropical alpine regions[END_REF].

Although there is a general consensus on the need to increase our knowledge on the climate change impacts on the páramo biota and ecosystems, studies remain very scarce to date [START_REF] Tovar | Diverging responses of tropical Andean biomes under future climate conditions[END_REF][START_REF] Feeley | Most 'global'reviews of species' responses to climate change are not truly global[END_REF][START_REF] Anthelme | Alpine and arctic plant communities: a worldwide perspective[END_REF][START_REF] Anthelme | Biogeography of South American Highlands. Reference Module in Earth Systems and Environmental Sciences, Encyclopedia of the World's biomes[END_REF]. For instance, [START_REF] Anderson | Consequences of climate change for ecosystems and ecosystem services in the tropical Andes[END_REF] predicted a general upslope migration of 600 m for species to track the isotherms throughout the Tropical Andes by 2100. Other studies have also predicted up to a 50-60% increase in threats and local extinction risks by the second half of the century in the region [START_REF] Ramírez-Villegas | Using species distributions models for designing conservation strategies of Tropical Andean biodiversity under climate change[END_REF]. At the ecosystem level, changes in structure and functions are suggested for the Tropical Andes (e.g. [START_REF] Tovar | Diverging responses of tropical Andean biomes under future climate conditions[END_REF][START_REF] Cuesta | New land in the Neotropics: a review of biotic community, ecosystem, and landscape transformations in the face of climate and glacier change[END_REF]) but lack proper quantification for the páramo alone, and the consequent impacts on species assemblages and ecosystem services remain unclear [START_REF] Buytaert | Potential impacts of climate change on the environmental services of humid tropical alpine regions[END_REF]. In order to fill the important knowledge gap on species redistribution under anthropogenic climate change and the subsequent structural re-organisation of plant assemblages in the páramo, new endeavors such as the Institute von Humboldt's Biomodels (http://biomodelos.humboldt.org.co) or the Pontificia Catolica Universidad del Ecuador's Bioweb (https://bioweb.bio), as well as individual research initiatives (e.g. [START_REF] Ramírez-Villegas | Using species distributions models for designing conservation strategies of Tropical Andean biodiversity under climate change[END_REF], looking to model the potential distribution of páramo species under climate change using species distribution models (SDMs) fitted on bioclimatic data are increasing.

Dispersal is key in enabling species to track climate change, essentially through the mechanisms of seed production, migration over certain dispersal distance and along suitable habitat routes, germination and establishment [START_REF] Kammer | Increasing species richness on mountain summits: Upward migrationdue to anthropogenic climate change or re-colonisation[END_REF][START_REF] Vittoz | Seed dispersal distances: A typology based on dispersal modes and plant traits[END_REF][START_REF] Lenoir | Climate-related range shifts-a global multidimensional synthesis and new research directions[END_REF]. This factor prevails in alpine areas where plants face many challenges to disperse upslope or across mountains, this due to i) often reduceddispersal capacity [START_REF] Vittoz | Diaspore traits discriminate good from weak colonisers on high-elevation summits[END_REF][START_REF] Morgan | Alpine plant species have limited capacity for longdistance seed dispersal[END_REF], ii) the presence of biogeographic barriers such as forests and valleys, and iii) steep gradients favouring certain dispersal paths [START_REF] Dirnböck | Habitat distribution models, spatial autocorrelation, functional traits and dispersal capacity of alpine plant species[END_REF][START_REF] Engler | The MIGCLIM R package-seamless integration of dispersal constraints into projections of species distribution models[END_REF]. Accounting for dispersal when studying changes in species distribution in tropical mountains over time is crucial, because in comparation to higher latitudes, plants' dispersal capacities tend to be greater and more diverse [START_REF] Chen | Seeds tend to disperse further in the tropics[END_REF] while spatial climatic variations are less pronounced. At the community level, the resulting persistence vs. migration mismatch between local species, competitive regional pool and newcomers, such as invasive species, translates in important restructuration and composition changes [START_REF] Singer | Community dynamics under environmental change: How can next generation mechanistic models improve projections of species distributions?[END_REF][START_REF] Alexander | Lags in the response of mountain plant communities to climate change[END_REF]. Despite the general agreement on the need to include dispersal capacity in SDMs to improve current and future species distribution predictions [START_REF] Araújo | Five (or so) challenges for species distribution modelling[END_REF][START_REF] Thuiller | Predicting global change impacts on plant species' distributions: future challenges[END_REF][START_REF] Dullinger | Extinction debt of high-mountain plants under twenty-first-century climate change[END_REF][START_REF] Hattab | A unified framework to model the potential and realized distributions of invasive species within the invaded range[END_REF], this factor remains overlooked today. Moreover, when accounted for, it usually takes the form of an "all or nothing" parameter, i.e. no dispersal or full dispersal [START_REF] Guisan | Predicting species distribution: offering more than simple habitat models[END_REF]Araújo et al. 2006;[START_REF] Holloway | Incorporating movement in species distribution models: how do simulations of dispersal affect the accuracy and uncertainty of projections?[END_REF], albeit partial dispersal scenarios based on mean and/or maximum dispersal distance often are more reliable [START_REF] Dullinger | Extinction debt of high-mountain plants under twenty-first-century climate change[END_REF][START_REF] Engler | The MIGCLIM R package-seamless integration of dispersal constraints into projections of species distribution models[END_REF][START_REF] Bateman | Appropriateness of full-, partial-and no-dispersal scenarios in climate change impact modelling[END_REF].

Nevertheless, partial dispersal approaches are increasing thanks to new techniques allowing to calculate significant approximations of dispersal distances based on a suite of plant trait data [START_REF] Vittoz | Seed dispersal distances: A typology based on dispersal modes and plant traits[END_REF][START_REF] Thomson | Seed dispersal distance is more strongly correlated with plant height than with seed mass[END_REF][START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF][START_REF] Bullock | A synthesis of empirical plant dispersal kernels[END_REF]). Among the most dependable traits used to date is plant height, assuming that a tall plant disperses further than a short one [START_REF] Muller-Landau | Interspecific variation in primary seed dispersal in a tropical forest[END_REF][START_REF] Thomson | Seed dispersal distance is more strongly correlated with plant height than with seed mass[END_REF]. Seed mass has also been evaluated as a significant trait, generally suggesting that light seeds disperse further than heavy ones [START_REF] Parolo | Upward migration of vascular plants following a climate warming trend in the Alps[END_REF], even though this pattern can be significantly influenced by the correlation between plant height and seed mass [START_REF] Moles | Small-seeded species produce more seeds per square metre of canopy per year, but not per individual per lifetime[END_REF][START_REF] Thomson | Seed dispersal distance is more strongly correlated with plant height than with seed mass[END_REF]. Finally, dispersal syndrome is usually considered important as, for instance flightenhancing structures may increase anemochorous dispersal, and fruit characteristics can influence zoochorous dispersal [START_REF] Vittoz | Seed dispersal distances: A typology based on dispersal modes and plant traits[END_REF][START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF][START_REF] Thomson | Can dispersal investment explain why tall plant species achieve longer dispersal distances than short plant species?[END_REF]).

Yet, to go beyond individual species predictions from SDMs, recent scientific advances in ecology propose to model entire biotic communities, relying on species distributions and adding a frame of abiotic and biotic factors such as historical-dispersal constraints, ecosystem richness capacity and biotic interactions among species [START_REF] Guisan | SESAM-a new framework integrating macroecological and species distribution models for predicting spatio-temporal patterns of species assemblages[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF][START_REF] Mod | Biotic interactions boost spatial models of species richness[END_REF]. In this context, the Spatially Explicit Species Assemblage Modelling framework (SESAM; [START_REF] Guisan | SESAM-a new framework integrating macroecological and species distribution models for predicting spatio-temporal patterns of species assemblages[END_REF] has been extensively used to predict species assemblages on a wide range of taxonomic groups, from insects to plants, and for several habitats, from forests to grasslands (e.g. D'Amen, Pradervand, & Guisan, 2015a;[START_REF] D'amen | Predicting richness and composition in mountain insect communities at high resolution: a new test of the SESAM framework[END_REF][START_REF] Mod | Biotic interactions boost spatial models of species richness[END_REF][START_REF] D'amen | Spatial predictions at the community level: from current approaches to future frameworks[END_REF][START_REF] Mateo | Biodiversity models: what if unsaturation is the rule?[END_REF]. A wealth of research efforts has focused on improving the different aspects of SESAM, either regarding the macroecological constraints that represent the carrying capacity of communities (D'Amen et al. 2015a;[START_REF] Mateo | Biodiversity models: what if unsaturation is the rule?[END_REF], or the assembly rules that prioritize certain species co-occurrences [START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF][START_REF] D'amen | Predicting richness and composition in mountain insect communities at high resolution: a new test of the SESAM framework[END_REF]). However, the historical-dispersal constraints remain overlooked to date and the need persists to improve techniques and launch more integrative models to refine community modelling at present and under climate change scenarios [START_REF] D'amen | Spatial predictions at the community level: from current approaches to future frameworks[END_REF][START_REF] Guisan | Scaling the linkage between environmental niches and functional traits for improved spatial predictions of biological communities[END_REF].

The present study aims at implementing the SESAM framework, accounting for dispersal constraints for the first time, to predict changes in plant distributions and assemblages of the páramo sky islands during the 21 st century. First, maximum dispersal distances were calculated and compared for mid-páramo and super-páramo species. It was expected that tall species such as shrubs that present either light seeds and dispersalenhancing structures, e.g. wings, or heavy seeds with fruits dispersed by efficient zoochory, e.g. by birds, would be the best dispersing species. Second, the current and future distributions of each species up to the year 2100 were modelled based on bioclimatic variables while accounting for dispersal constraints to explore migration and extinction rates in the páramo. It was assumed that mid-páramo species would migrate further and faster than super-páramo species, which would in turn present greater extinction risks. Third, current and future plant communities were assembled on the sky islands by stacking all previous probabilistic models and applying richness constraints and assembly rules. The vulnerability of the network of sky islands was assessed based on the magnitude of changes in species richness and species composition in terms of species numbers and ratio of midpáramo and super-páramo species, expecting that locally diverse areas today [START_REF] Peyre | Fine-scale plant richness mapping of the Andean páramo according to macroclimate[END_REF] would suffer faster and more drastic changes.

Methods

Study area

Ecuador was set as study case as a model representant of páramo countries due to its abundance of páramos and equatorial geolocation. The most representative mountains carrying the super-páramo belt, or sky islands, in Ecuador are from North to South: Cayambe, Antisana, Pichincha, Ilinizas, Cotopaxi, Chimborazo, Tungurahua, Sangay and Cajas (Fig. 1A). The super-páramo is usually found between 4200-4800 m a.s.l, although its lower and upper limits as well as main vegetation strongly depend on local topography and topoclimate. For example, low shrublands and cushion communities can be found at lower elevations under a dense upper cloud elevation zone, while meadows, deserts with shrubs and periglacial deserts are most likely to occur in the higher areas shaped by volcanic activity and/or glacier melting (Fig. 1B; [START_REF] Sklenář | Diversity of zonal páramo plant communities in Ecuador[END_REF]. To conduct the following modelling analyses, the 4200 m isohypse and mountain tops were assumed as the lower and upper limits of the actual super-páramo belt respectively. Similarly, the 3000 m isohypse was set to delimit the broader páramo province (mid-páramo plus super-páramo), in order to consider upward migration of plant species from the mid-páramo to the super-páramo as climate will warm under future scenarios.

Vegetation data

All vegetation plots corresponding to the Ecuadorian páramo were downloaded from

VegParamo, an open access database for floristic and vegetation páramo data, compiling information from 40 data sources (www.vegparamo.com;[START_REF] Peyre | VegPáramo, a flora and vegetation database for the Andean páramo[END_REF]. Because VegParamo comprises the original authors' classification of the plots into sub-páramo, midpáramo, super-páramo and azonal vegetation, i.e. Polylepis forests, bogs and marshes and rupicolous vegetation, only the plots belonging to the zonal mid-páramo and super-páramo elevational belts were retained. Plots located outside the study area, e.g. on Amazonian mountains, were also removed, as well as plots with a coarse georeferencing precision superior to 1km in the UTM system. The final plots, all sampled according to the phytosociological method [START_REF] Braun-Blanquet | Pflanzensoziologie: Grundzüge der Vegetationskunde[END_REF], were rescaled into presence/absence data.

The floristic contents of the plots were checked, and the following groups were eliminated from the dataset: pteridophytes, bryophytes, lichens and vascular plants determined at the genus level or higher. When a taxon was determined at the infra-specific level, such as varieties and subspecies, it was aggregated at the species level. Taxonomic synonymy was checked using the Plant List (www.plantlist.com) and Tropicos (www.tropicos.org). Because the consequent vegetation dataset showed an important imbalance of species presences versus absences, additional occurrence data was complemented for these species from biological and herbarium databases consulted online, including Tropicos, Aarhus University Herbarium (www.aubot.dk) and GBIF (www.gbif.org). The additional occurrence data was revised to remove outlier points and check duplicates between occurrence data and vegetation plots (based on UTM coordinates).

Decimal spatial coordinates were then obtained for all plots, using the centroid of the corresponding UTM cell, as well as for the new added occurrence data, relying on the original georeferencing. Species that still presented less than 6 occurrences were removed and the complete dataset includes 642 vegetation plots from VegParamo (157 UTM) and 2095 additional occurrences for 435 vascular plant species (Fig. S1; Table S1, S2; Peyre, 2020). Finally, species were classified into super-páramo and mid-páramo species, when at least 50% of their occurrence data occurred above and below the 4200 m isohypse, respectively, resulting in 84 super-páramo species and 351 mid-páramo species.

Bioclimatic data 19 bioclimatic variables were downloaded from the CHELSA V1.2 database (www.chelsaclimate-org; [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF]) and adjusted to the study area. All variables were obtained for the period 1979-2013 and averaged across this time interval, considered as the current conditions (year 2000), which fitted the occurrence and vegetation data .

Multicollinearity between variables was evaluated using a variance inflation factor correlation analysis (vif function, usdm R package; [START_REF] Naimi | Where is positional uncertainty a problem for species distribution modelling?[END_REF], and all variables below a threshold value of 0.7 for the Pearson correlation coefficient were retained (Dormann et al. 2012), resulting in the following selection: mean diurnal temperature range (bio2); temperature seasonality (bio4); mean temperature of the wettest quarter (bio8); precipitation seasonality (bio15); and precipitation of the coldest quarter (bio19). This set of predictor was not fitted per species (D'Amen et al. 2015a;[START_REF] Araújo | Standards for distribution models in biodiversity assessments[END_REF]) but was considered appropriate for the study area as a whole because of its certain focus on precipitation-related factors and seasonality, which usually prevail as drivers of plant diversity variation in tropical mountain areas [START_REF] Peyre | Fine-scale plant richness mapping of the Andean páramo according to macroclimate[END_REF]. To represent future scenarios and in order to encompass sufficient variance while reducing uncertainty in predictions [START_REF] Thuiller | Uncertainty in ensembles of global biodiversity scenarios[END_REF], a total of 8 climate change scenarios were selected, based on two representative concentration pathways (RCPs): CIMC5-RCP45 and CIMC5-RCP85, as well as four global circulation models (GCMs): bbc-csm1-1, CESM1-BGC, HadGEM2-AO and MRI-CGCM3. These GCMs were privileged based on their dissimilarity, according to the distance matrix calculated for all available GCMs presented in [START_REF] Sanderson | A Representative Democracy to Reduce Interdependency in a Multimodel Ensemble[END_REF]. Future bioclimatic predictions for each variable of interest were then downloaded from CHELSA V1.2 for every scenario at two dates, 2050 and 2070, and later cropped to fit the study area. Because dispersal is considered here as a key factor to constrain future species distributions, the bioclimatic data was interpolated by decade to provide frequent steps to include dispersal limitations [START_REF] Engler | The MIGCLIM R package-seamless integration of dispersal constraints into projections of species distribution models[END_REF]). To do so, simple linear regressions were used to obtain for each bioclimatic variables the 2010, 2020, 2030 and 2040 values using a coefficient calculated for the 2000-2050 period following a similar procedure to the one employed by [START_REF] Adhikari | Potential impact of climate change on the species richness of subalpine plant species in the mountain national parks of South Korea[END_REF]. The same procedure was applied for 2060 based on the 2050-2070 period, and for decades 2080, 2090 and 2100, values by extrapolating the 2050-2070 regression coefficient into the future. All statistical analyses were conductet in R 4.3.1 (R [START_REF] Core | R: A language and environment for statistical computing[END_REF].

Dispersal capacity

The widely-used dispeRsal function developed by [START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF] was chosen to estimate dispersal capacity in form of maximum dispersal distance, which is usually more suitable to include in SDMs than mean distance [START_REF] Vittoz | Seed dispersal distances: A typology based on dispersal modes and plant traits[END_REF][START_REF] Thomson | Seed dispersal distance is more strongly correlated with plant height than with seed mass[END_REF]. Beforehand, a basic trait database including plant height, vegetation stratum, main dispersion mode and seed mass was compiled for the 435 species. For plant height and vegetation stratum, information at the species level was retrieved from online herbarium material, such as the Field Museum (https://plantidtools.fieldmuseum.org), COL herbarium (www.biovirtual.unal.edu.co) and Plant JSTOR (https://plants.jstor.org), directly measured and averaged using at least four specimens per species when available. Dispersion mode and seed mass data at the species level were usually unavailable and therefore, information provided at the genus level from the Kew Botanical Garden seed collection initiative (http://data.kew.org/sid), and specific páramo literature were used [START_REF] Frantzen | Dispersal and growth form patterns of some zonal páramo vegetation types[END_REF][START_REF] Melcher | Seed dispersal in páramo plants: epizoochorous and hydrochorous taxa[END_REF], 2004). Assuming for simplicity reasons that species would disperse only once yearly, the dispeRsal function was run as a kernel-shape probability density function of dispersal distances to calculate maximum dispersal distances (in m/yr) [START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF][START_REF] Bullock | A synthesis of empirical plant dispersal kernels[END_REF]. The raw results were then upscaled to km per decade to fit the future bioclimatic predictions. Finally, dispersal capacity of midpáramo and super-páramo species were sommely compared by means of a t-Student test.

Species distribution models

The current distribution of each species was modelled based on actual environmental conditions, represented by the 5 bioclimatic variables selected, and dispersal capacity. The latter factor was introduced as a new variable, set by log-transforming the species maximum dispersal distance previously computed from plant traits and scaling the obtained values between 0 and 1 to rank species from bad to good dispersers (function iForce, iSDM package; [START_REF] Hattab | A unified framework to model the potential and realized distributions of invasive species within the invaded range[END_REF]. SDMs were then run using two different families of algorithms fitted for presence-absence data: i) generalized linear models (GLMs), which explore the linear relationship between explanatory and explained variables through a flexible Gaussianidentity distribution-link approach and were fitted using second degree polynomial curves and ii) random forests (RFs) that integrate the data through an ensemble learning, coupling decision trees with classification and regression settings/approaches (biomod2 package; [START_REF] Thuiller | Package 'biomod2'. Species distribution modeling within an ensemble forecasting framework[END_REF]. For each species, the data was divided into a 75% training dataset and a 25% testing dataset, and run with 50 models per algorithm. The true skill statistic (TSS) value was computed for each model to evaluate SDMs' performances and all models with a TSS value greater than 0.6, either obtained with the GLM or RF algorithms, were retained [START_REF] Allouche | Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS)[END_REF]. The final models were then ensembled to obtain the final probabilistic preditions [START_REF] Araújo | Ensemble forecasting of species distributions[END_REF][START_REF] Marmion | Evaluation of consensus methods in predictive species distribution modelling[END_REF], which were finally transformed into binary values using a threshold approach that equalizes the sensitivity and specificity metrics (optimal.threshold function; ecospat package; [START_REF] Liu | Selecting thresholds of occurrence in the 551 prediction of species distributions[END_REF][START_REF] Di Cola | ecospat: a R package to support spatial analyses and modeling of species niches and distributions[END_REF].

The calibrated current models were then used to project species' distributions into the future per decade and under each of the 8 different climate-change scenarios. To do so, the previoulsy-used dispersal covariate was set to 0 and future distributions were projected against bioclimatic predictors only, in order to alleviate the static species-specific dispersal constraint. As formerly done, the probabilistic projections for each decade were then transformed into binary values and constrained in temporal order by the species maximum dispersal distance [START_REF] Engler | The MIGCLIM R package-seamless integration of dispersal constraints into projections of species distribution models[END_REF], hence allowing to track progressive changes in species distributions depending on available bioclimatic conditions and capacity to disperse from the previous distribution. Additional biogeographical barriers represented by the Andean forest and anthropogenic activities from the lowlands were also considered by using a raster mask at each time step that assumed that all pixels below the 3000 m isohypse were unsuitable. Thus, a species could only cross the matrix of unsuitable conditions in a timestep if its maximum dispersal distance was greater than the distance separating two suitable pixels. The set of decadal binary maps obtained for each species separately offered time series of species redistribution that can be used to assess the species-specific proportion of areas that were lost, gained or that remained occupied by the focal species between present and 2100 (functions migclim.distance and migclim.plot, migclim package; [START_REF] Engler | The MIGCLIM R package-seamless integration of dispersal constraints into projections of species distribution models[END_REF]. Finally, the 2100 binary prediction was used to fit the 2100 probabilistic prediction for richness capacity (see below).

Changes in species distribution for the 21 st century were evaluated as a binary result, net loss or net gain, by comparing and quantifying the number of pixels belonging to the 2000 and 2100 distributions. Finally, extinction risk was considered when a species had at least one scenario predicting its complete extinction (i.e. the species completely disappeared from the studied area) by 2100. The severity of the risk was evaluated based on the number of scenarios predicting extinction by 2100, i) low risk (1-2 scenarios), ii) intermediate risk (3-4 scenarios), high risk (5-6 scenarios) and extinct (7-8 scenarios).

Model richness and assemblages

The SESAM framework relies on the predict first assemble later principle [START_REF] Guisan | SESAM-a new framework integrating macroecological and species distribution models for predicting spatio-temporal patterns of species assemblages[END_REF], meaning that final community predictions are achieved after applying an assembly procedure that consists in selecting species from the predicted pool until reaching the carrying capacity of a community defined as potential richness [START_REF] Guisan | SESAM-a new framework integrating macroecological and species distribution models for predicting spatio-temporal patterns of species assemblages[END_REF][START_REF] Mateo | Biodiversity models: what if unsaturation is the rule?[END_REF]. This local carrying capacity (at pixel level) is often estimated by either stacking the probabilities predicted by the SDMs or by using macroecological models (D'Amen et al. 2015a;[START_REF] Mateo | Biodiversity models: what if unsaturation is the rule?[END_REF] but the former has the advantage of being embedded within a same S-SDM framework and allowing for a varying timeline, so that we used it here. All SDMs in their probabilistic form were stacked in 2000 and 2100 for each scenario to evaluate richness. Net richness change at the pixel level was assessed as a percentage of the current (2000) richness, with positive and negative values meaning a net gain and loss, respectively, in potential species richness.

Finally, species composition at the pixel level was predicted through the probability ranking rule (PRR; D'Amen et al. 2015b), which consists in selecting species from the predicted pool in decreasing order of their SDM probability of occurrence until reaching the previously predicted potential richness. The PRR was applied to the stacks of SDM predictions for 2000 and 2100 using the prr.ecospat function of the ecospat package [START_REF] Di Cola | ecospat: a R package to support spatial analyses and modeling of species niches and distributions[END_REF]. Changes in composition, expressed as numbers of species gained and lost, as well as ratio of mid-páramo and super-páramo species to emulate upward migration and high-elevation competition, were evaluated for the 10 main páramo sky islands. Lastly, based on the estimated changes in richness and composition, local vulnerability to climate change was assessed.

Results

Species varied greatly in terms of their maximum dispersal distance, ranging from 0.008 to 6027 m/yr in a negative exponential manner with a mean of 193 m/yr (± 647) (Fig. S2; Table S2). Mid-páramo and super-páramo species showed significantly different maximum dispersal distances (t-test: 4.2972, df: 423.26, p-value: 2.148e-05), with mid-páramo species dispersing relatively well in average at 227.14 m/yr (± 714), while super-páramo species showed shorter dispersal distances at 51.39 m/yr (± 136). Long dispersal distances were usually associated with shrub genera with relatively heavy seeds (10-100 mg) and using zoochory or assisted anemochory, such as Monnina (Polygalaceae) and Gaiadendron (Loranthaceae), whereas short dispersal distances was common in small herbaceous plants with lightweight seeds (0.0001-0.001 mg) and using unassisted anemochory or autochory as principal dispersal mode, as seen in Aa (Orchidaceae) and Ourisia (Plantaginaceae).

The performance of SDMs varied between species, and the average values were considered acceptable overall (across models: mean TSS = 0.614, Sensitivity = 79.835, Specificity = 78.952; Table 1). In general, the RF algorithm performed better than the GLM algorithm, and more of its models were retained, approximately 29% across all species. About 70% of the SDMs across all scenarios predicted a net loss in species' spatial distribution and only 30% predicted a net gain (Fig. 2). A total of 65 species (15% of all species) presented a degree of extinction risk under at least one scenario by the end of the century, 18 of which as a low risk, eight as an intermediate risk, 15 as a high risk and finally 39 were predicted to become totally extinct from the studied area by 2100. The most endangered genera encountered in numbers were Cerastium (3 species), Draba (3 species), Gentianella (4 species) and Viola (3 species). In addition, 22 of the 65 species were superpáramo species, 5 of which are known endemics at the national scale: Aetheolaena involucrata, Bartsia pumila, Festuca chimborazensis, Loricaria antisanensis and Viola polycephala. Finally, of the 39 species with the highest extinction risk, 10 were superpáramo species known as páramo endemics, for example Cerastium candicans, Draba depressa, Lachemilla tanacetifolia, Lupinus alopecuroides and Xenophyllum crassum.

Regarding net richness changes in the páramo sky islands, the general trend at pixel level (alpha diversity) was a loss of 56.63 % (± 20.83) of the original species richness between the present and 2100 (Fig. 3). The most affected areas were the northern mountains of Pichincha and Ilinizas with an average loss of 82.64 (± 2.15) and 77.24% (± 3.95) respectively. Contrarily, the least affected areas were the central mountains of Tungurahua and Altar with an average loss of 13.52 % (± 11.12) and 31.90 % (± 6.52) respectively, and punctual pixels with positive values hence gaining species, up to 18.73 % on the Chimborazo.

At the sky island level, important significant composition changes in the final plant communities were observed (Fig. 3; Fig S3; Table S3). Sky islands lost an average 64.74% (± 14.12) of their original species pool (gamma diversity), Pichincha and Ilinizas being the most significantly affected with losses of 85.68% (± 5.70) and 74.73% (± 5.12) respectively. By contrast, Altar and Tungurahua were the least affected, losing 46.49 (± 10.76) and 50.00% (± 9.74) respectively. Gains at this scale did not compensate losses, with the average gain reaching 12.97 % (± 11.41), some páramos gaining more species such as Tungurahua and Altar with 39.21% (± 5.03) and 20.86% (± 6.31), while others gained less, such as Cayambe and Antisana with 7.17 % (± 3.02) and 7.47 % (± 3.05). The most drastic changes were therefore observed for the northern mountains of Cayambe and Pichincha.

On average in 2000, the ratio between mid-páramo and super-páramo species in the páramo sky islands reached 68.79 / 31.21 (± 6.09), and changed across all scenarios in 2100 to 73.61 / 26.39 (± 4.46), meaning that approximately 4.82% of super-páramo species were replaced with mid-páramo ones (Table 2). The most drastic tendencies were observed for Ilinizas and Antisana, losing 9.76 and 8.01% of their super-páramo species respectively, while the least affected sky islands were Cayambe, which gained 0.78% of super-paramo species, and Sangay that lost 0.48% of super-paramo species. Based on richness at pixel and sky island level and composition changes, Ilinizas, Pichincha and in lesser means Antisana were therefore considered as most threatened by climate change.

Discussion

Mountains are particularly sensitive to climate change because they are often prone to accelerated changes and more rapid responses of biota in comparison with lowlands [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF][START_REF] Smith | Near-term acceleration in the rate of temperature change[END_REF][START_REF] Steinbauer | Accelerated increase in plant species richness on mountain summits is linked to warming[END_REF]. Neotropical mountains in particular are classified as highly vulnerable to climate change [START_REF] Team | climate change 2014: synthesis report[END_REF] and are expected to suffer drastic changes in biodiversity and its associated ecosystem services that will affect millions of people on the short to mid-term [START_REF] Buytaert | Potential impacts of climate change on the environmental services of humid tropical alpine regions[END_REF][START_REF] Tovar | Diverging responses of tropical Andean biomes under future climate conditions[END_REF][START_REF] Anthelme | Alpine and arctic plant communities: a worldwide perspective[END_REF]. In addition, the fast and important changes in rural occupation and land uses that are occurring due to improving socio-economic conditions and the subsequent population growth, contribute greatly to the extent and intensity of environmental degradation [START_REF] Hofstede | Los páramos del mundo: Proyecto Atlas Mundial de los Páramos[END_REF].

It is also in this region that predicting the impacts of future global changes becomes particularly challenging, due to data availability and the lack of knowledge on biodiversity, climatic variations, anthropogenic dynamics and the responses of biodiversity to these drivers [START_REF] Anderson | Consequences of climate change for ecosystems and ecosystem services in the tropical Andes[END_REF][START_REF] Ramírez-Villegas | Using species distributions models for designing conservation strategies of Tropical Andean biodiversity under climate change[END_REF][START_REF] Lenoir | Climate-related range shifts-a global multidimensional synthesis and new research directions[END_REF][START_REF] Feeley | Most 'global'reviews of species' responses to climate change are not truly global[END_REF].

Dispersal capacity

To date, trait data and functional analyses for the páramo remain scarce, especially regarding non-leaf characteristics (e.g. [START_REF] Sánchez | Dynamic cloud regimes, incident sunlight, and leaf temperatures in Espeletia grandiflora and Chusquea tessellata, two representative species of the Andean Páramo, Colombia. Arctic[END_REF], and so far few works have focused on seed-related traits and characterized the dispersal modes of páramo species [START_REF] Frantzen | Dispersal and growth form patterns of some zonal páramo vegetation types[END_REF][START_REF] Melcher | Seed dispersal in páramo plants: epizoochorous and hydrochorous taxa[END_REF][START_REF] Melcher | Seed atlas of the monocotyledonous genera of the páramo[END_REF]. Our estimates of species maximum dispersal distances based on plant dispersal traits highlighted that páramo species disperse on the lower end of the terrestrial angiosperm dispersal spectrum [START_REF] Kinlan | Propagule dispersal in marine and terrestrial environments: a community perspective[END_REF] and in general similarly to other alpine floras [START_REF] Vittoz | Seed dispersal distances: A typology based on dispersal modes and plant traits[END_REF][START_REF] Morgan | Alpine plant species have limited capacity for longdistance seed dispersal[END_REF]. Our result showed some support for the hypothesis that low-elevation shrubs, characterized by heavy seeds and usually using either zoochory (e.g. Ericaceae) or assisted anemochory (e.g. Asteraceae), disperse best. Indeed, the correlation favouring plant height over seed mass seemed particularly advantageous for these shrub species, at least at the landscape level [START_REF] Parolo | Upward migration of vascular plants following a climate warming trend in the Alps[END_REF]Thomson et al. 2011, 17), whereas within a same dispersal syndrome or vegetation type, light-seed species that disperse anemochorously are usually considered as more efficient [START_REF] Muller-Landau | Interspecific variation in primary seed dispersal in a tropical forest[END_REF]. Our best dispersers were often zoochorous and presenting relatively large fruits visible to vertebrate animals, such as birds, which undoubtley enhanced local dispersal [START_REF] Thomson | Seed dispersal distance is more strongly correlated with plant height than with seed mass[END_REF][START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF] and also increased frequency of long-distance dispersal events, for example between mountain tops [START_REF] Vittoz | Diaspore traits discriminate good from weak colonisers on high-elevation summits[END_REF]). Consequently, an upcoming challenge for these plant species in the future is facing their own responses to climate change but also those of their animal disperser [START_REF] Travis | Dispersal and species' responses to climate change[END_REF]. With increasing elevation and distance from the treeline, Poaceae become progressively dominant in the páramo ecosystems, ensuing a larger proportion of anemochorous species that often present disperal-enhancing structures, i.e. wings or pappi. Such structures are also crucial in helping plants overcome the long-distance dispersal limitations they face in alpine areas [START_REF] Vittoz | Diaspore traits discriminate good from weak colonisers on high-elevation summits[END_REF][START_REF] Morgan | Alpine plant species have limited capacity for longdistance seed dispersal[END_REF][START_REF] Thomson | Can dispersal investment explain why tall plant species achieve longer dispersal distances than short plant species?[END_REF]). In the super-páramo belt, small characteristic plants such as Nototriche, Draba and Ourisia [START_REF] Peyre | Phytoregionalisation of the Andean páramo[END_REF]) are often featuring perreniality, small light-seeds, high autogamy rate and little seed productivity, which might explain the shortdistance dispersal tendency observed. In addition, the yet unhospitable environment and unavailability of efficient animal dispersers on these sky-islands often forces plant species to recur to vegetative (clonal) reproduction in priority and limit their sexual reproduction [START_REF] Vittoz | Diaspore traits discriminate good from weak colonisers on high-elevation summits[END_REF]. Accounting for establishment capacity, the availability of microsites, good microclimatic and edaphic conditions, and strong species interactions, such as facilitation of assisted dispersal, is therefore key to finally assess the migration capacity of these species [START_REF] Scherrer | Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming[END_REF][START_REF] Anthelme | Facilitation among plants in alpine environments in the face of climate change[END_REF][START_REF] Hupp | Alpine cushion plants have species-specific effects on microhabitat and community structure in the tropical Andes[END_REF][START_REF] Alexander | Lags in the response of mountain plant communities to climate change[END_REF]. Finally, considering the competitiveness of super-páramo species against advantaged shrub and herbaceous species from lower elevations, including opportunistic and exotic species, is necessary to evaluate the intensity of the threat (Pauchard et al. 2009;[START_REF] Llambí | Reciprocal interactions between a facilitator, natives, and exotics in tropical alpine plant communities[END_REF].

The importance of dispersal capacity in shaping a species distribution has long been acknowledged and increasingly included to improve SDM performances in the last decade [START_REF] Thuiller | Predicting global change impacts on plant species' distributions: future challenges[END_REF][START_REF] Dullinger | Extinction debt of high-mountain plants under twenty-first-century climate change[END_REF][START_REF] Bateman | Appropriateness of full-, partial-and no-dispersal scenarios in climate change impact modelling[END_REF][START_REF] Hattab | A unified framework to model the potential and realized distributions of invasive species within the invaded range[END_REF].

Community modelling has been considerably improving in parallel, with considerable efforts being put on better accounting for community assembly processes, including the carrying capacity of communities and assembly rules [START_REF] D'amen | Spatial predictions at the community level: from current approaches to future frameworks[END_REF][START_REF] Mateo | Biodiversity models: what if unsaturation is the rule?[END_REF], nonetheless, dispersal constraints had remained mostly overlooked in community modelling [START_REF] D'amen | Spatial predictions at the community level: from current approaches to future frameworks[END_REF][START_REF] Guisan | Scaling the linkage between environmental niches and functional traits for improved spatial predictions of biological communities[END_REF]). This study supports that dispersal is a key factor in defining the future distribution of species and their assemblages, as illustrated in Fig. 2 by the overall overprediction of potential versus realized distributions quantified as the suitable but unoccupied area. Our finding is therefore in agreement with previous works maintaining that dispersal capacity is particularly important in tropical mountains, where climate has had an important but not exclusively prevailing influence on biota distribution (e.g. [START_REF] Flantua | The flickering connectivity system of the north Andean páramos[END_REF]. Our results also suggest that super-páramo plants might not be able to migrate altitudinally along the same mountain but also be prevented in their attempts to colonize new mountain ranges, sustaining the key role of topographic distance and isolation [START_REF] Sklenář | Distribution patterns of páramo plants in Ecuador[END_REF][START_REF] Flantua | The flickering connectivity system of the north Andean páramos[END_REF].

Distribution responses of plant species

Of all species, 70% were predicted to show a significant distribution loss by 2100, which complements previous findings in which 50-60% of tropical Andean species reduced their distribution in the second half of the 21 st century [START_REF] Ramírez-Villegas | Using species distributions models for designing conservation strategies of Tropical Andean biodiversity under climate change[END_REF]. Even more concerning is that 15% of all páramo species showed a significant extinction risk by 2100, and 25% of all super-páramo species, supporting our hypothesis that super-páramo species are generally more at risk than mid-páramo ones. Among the vulnerable super-páramo species are several diagnostic species of super-páramo phytogeographical units, for example Ourisia muscosa of the Upper humid super-páramo and Draba depressa of the Upper dry super-páramo [START_REF] Peyre | Phytoregionalisation of the Andean páramo[END_REF]. Furthermore, our results posed as threatened many superpáramo species that are páramo endemics such as Cerastium candicans, Draba depressa and Xenophyllum crassum or even national páramo endemics, for example Bartsia pumila, Loricaria antisanensis and Viola polycephala. Our study alerts of maximum extinction risk for 9% of all studied species, some of which were already classified following the IUCN system as threatened, for example Bartsia pumila (VU) and Gentianella hirculus (EN), but also some considered as not threatened, such as Gentianella limoselloides (LC) and Halenia taruga-gasso (NT) (León-Yánez, 2012). Because SDMs can be a useful complement to further conduct UICN evaluations [START_REF] Breiner | Including environmental niche information to improve IUCN Red List assessments[END_REF]; but see [START_REF] Akcakaya | Use and misuse of the IUCN Red List Criteria in projecting climate change impacts on biodiversity[END_REF], we therefore urge, based on the SDM results to focus further study on these predicted-extinct species, especially regarding population health and individual fitness, and encourage ex-situ conservation efforts such as seed banks to preserve this unique biodiversity.

Changes in richness and plant assemblages

Local richness at the pixel level has been found to significantly decrease in most sky islands, even accounting for potential colonization from lower elevations. The northern sky islands of Ilinizas and Pichincha were the most affected, which contradicts previous findings stating richness stability in this region [START_REF] Ramírez-Villegas | Using species distributions models for designing conservation strategies of Tropical Andean biodiversity under climate change[END_REF], and exceptions to the rule included few central mountains such as Tungurahua and Altar where richnesst would remain stable or even increase slightly. As a result, these last sky islands might be experiencing the worldwide-observed accelerated colonization process on mountains [START_REF] Steinbauer | Accelerated increase in plant species richness on mountain summits is linked to warming[END_REF]), but buffered by specific environmental conditions and a particularly vulnerable species pool. Overall, we found no support for the negative relationship hypothesis between actual richness and future richness changes, as the more diverse (at the local scale, or plot scale) eastern mountains suffered less drastic changes than the less diverse northern mountains [START_REF] Peyre | Fine-scale plant richness mapping of the Andean páramo according to macroclimate[END_REF], Because this finding plausibly supposes that the communities' richness capacity is not reached on these sky islands, we recommend further scientific focus on that particular macroecological hypothesis.

Composition changes were also predicted as very abrupt in the region, estimating important species losses and little gains at the sky island level (gamma diversity). Northern mountains were predicted to lose substantial amounts of species as much as 85% in Pichincha and gain very few species overall. By contrast, central mountains saw moderate losses and in the case of Tungurahua, a potential even loss and gain, depending on the climate change scenario considered. Several complementary hypotheses to climate change might be advanced to explain such a pattern, for example the geographic isolation and elevation of the massif, which could either facilitate or constrain dispersal and migration [START_REF] Sklenář | Distribution patterns of páramo plants in Ecuador[END_REF]. Finally, the ratio of super-páramo vs. mid-páramo species decreased by 2100, especially in northern mountains, sky islands generally losing 5% of their super-páramo species to mid-páramo ones. Our results therefore support previous findings that mid-páramo species are moving upslope and competing with super-páramo species, although the role of adaptation, microenvironmental refugia and biotic interactions remain to be considered for more refined plant responses.

Our results evidenced that sky islands from northern Ecuador, including Pichincha and Ilinizas, were evaluated as the most vulnerable to future climate change. These páramos are located along the inter-Andean valley where an important population resides, i.e. Quito and Latacunga, making the anthropogenic change expected for the end of the century in the form of fast and intense land-use changes an even bigger threat than climate change [START_REF] Hofstede | Los páramos del mundo: Proyecto Atlas Mundial de los Páramos[END_REF]. Because the agricultural frontier could move upslope too when soil develops sufficiently at high altitude, entire super-páramo ecosystems could suffer drastic changes in their structure and functions. Therefore, the predictions provided here, based on climate-change solely, should be considered as conservative scenarios and even more forceful impacts should be expected under a global change scenario. In addition, exotic and efficiently dispersed opportunistic species could take advantage of the anthropogenic change to migrate and invade the sky islands, competing directly with local species and forming noanalog communities with little ecological value [START_REF] Le Roux | Rapid range expansion and community reorganization in response to warming[END_REF][START_REF] Anthelme | Biogeography of South American Highlands. Reference Module in Earth Systems and Environmental Sciences, Encyclopedia of the World's biomes[END_REF]. Precise monitoring should therefore be tracking the anthropogenic changes associated to climate change on these sky islands in priority so to prevent a critical threat to the sky islands of the northern Andes currently holding the reservoir or museum of superpáramo species.

Study limitations and future recommendations

Regarding dispersal estimates, we consider that even though trait-based dispersal distance methods such as [START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF] are very useful and give reasonable predictions, the fact that they have been shown to underestimate distances for wind-dispersed species, an important proportion of our species pool [START_REF] Herrmann | Connectivity from a different perspective: comparing seed dispersal kernels in connected vs. unfragmented landscapes[END_REF][START_REF] Morgan | Alpine plant species have limited capacity for longdistance seed dispersal[END_REF], would support the use of more detailed techniques in the future, for example mechanistic models that acknowledge wind velocities and dispersal pathways [START_REF] Skarpaas | How far can a hawk's beard fly? Measuring and modelling the dispersal of Crepis praemorsa[END_REF][START_REF] Holloway | Incorporating movement in species distribution models: how do simulations of dispersal affect the accuracy and uncertainty of projections?[END_REF]. Moreover, we did not account for certain crucial variation such as i) potential changes in biogeographical barriers with climate change such as glacier retreat or forest expansion, which might affect the dispersal reach of species [START_REF] Caplat | Looking beyond the mountain: dispersal barriers in a changing world[END_REF], nor did we consider ii) that specific dispersal distances might change over time [START_REF] Travis | Dispersal and species' responses to climate change[END_REF] or iii) the importance of long dispersal events promoted by climate change consequences such as the frecuency and intensity of storms [START_REF] Hellmann | Five potential consequences of climate change for invasive species[END_REF]. Another limitation regards data availability and coverage, a commonly encountered issue in tropical research. In fact, the species list was recovered from VegParamo, which is a pioneer database with important representativity of the páramo taxa, although it still accounts for certain geographic and taxonomic knowledge gaps [START_REF] Peyre | VegPáramo, a flora and vegetation database for the Andean páramo[END_REF][START_REF] Bottin | Phytosociological data and herbarium collections show congruent large-scale patterns but differ in their local descriptions of community composition[END_REF]. Finally, climatic interpolations in tropical mountain areas that are topographically and environmentally complex are particularly challenging to obtain, especially for a broad study area and future climatic scenarios, and more fitted climatic data with a finer resolution could have improved the models performance and predictions.

Future recommendations include complementing the trait data with observed data and calculate distance on site to provide better estimates of dispersal capacity of plants in the region. To improve the models' accuracy and ecological meaning, it would be useful to account for different dynamic factors unavailable at the time [START_REF] Anthelme | Facilitation among plants in alpine environments in the face of climate change[END_REF][START_REF] Graae | Stay or go-how topographic complexity influences alpine plant population and community responses to climate change[END_REF], for example: i) additional abiotic variables at the macroscale, such as horizontal precipitation that represents an important water intake for plants, but also meso and microscale, for example local night freezing, that could balance macroeffects [START_REF] Scherrer | Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming[END_REF][START_REF] Mod | Biotic interactions boost spatial models of species richness[END_REF], ii) fine biotic interactions to incorporate into the assembly rules, competition but also facilitation to dispersal or establishment [START_REF] Llambí | Reciprocal interactions between a facilitator, natives, and exotics in tropical alpine plant communities[END_REF], and iii) land-use changes, whose intensity and expansion might increase under future climate change [START_REF] Hofstede | Los páramos del mundo: Proyecto Atlas Mundial de los Páramos[END_REF][START_REF] Anderson | Consequences of climate change for ecosystems and ecosystem services in the tropical Andes[END_REF].

Tables Table 1: Evaluation metrics and their standard deviation (in parenthesis) showing the overall performance of the GLM and RF models for 435 plant species as well as their comparison between algorithm according to a paired t-student test.

Figure captions Figure 1 :

 captions1 Figure captionsFigure 1: Current potential distribution of the páramo in Ecuador (A): in dark gray, midpáramo (3000-4200 m) and in red, super-páramo (> 4200 m); and characteristic superpáramo vegetation (B), from bottom to top, increasing with elevation.

Figure 2 :

 2 Figure 2: Examples of distribution responses of páramo species to climate change by 2100, according to the CIMC5-RCP45-CESM1-BGC scenario: A) a mid-páramo species colonizing higher elevations (Vicia andicola) and B) a super-páramo species reducing its distribution (Astragalus geminiflorus).

Figure 3 :

 3 Figure 3: Net richness changes at the pixel level (1 km 2 ) between 2000 and 2100 in the 10 páramo sky islands above 4200 m (named in bold). Between parenthesis, balance of gained and lost species between the present and 2100 (in percentage of the original richness at the sky island level).

Figure S1 :

 S1 FiguresFigure1

Figure S2 :

 S2 Figure S2: Density plot of the dispersal capacity, as maximum dispersal distance, of the páramo species (m/yr), differentiating the mid-páramo and super-páramo species.

Table 2 :

 2 Composition changes (in %) in terms of proportions of mid-páramo and superpáramo species at the sky island level between 2000 and 2100. The final super-páramo balance shows net gain (in green) or loss (in red) of super-páramo species, expressed as a percentage.

		2000		2100			Final super-
		Mid-páramo	Super-paramo	Mid-páramo	Super-páramo	Standard	páramo balance
		species	species	species	species	deviation
	Altar	66.91	33.09	74.54	25.46	4.15	-7.63
	Antisana	67.30	32.70	75.30	27.70	6.67	-8.01
	Cajas	81.82	18.18	84.21	15.79	6.10	-2.39
	Cayambe	69.33	30.67	68.56	31.44	4.16	+0.78
	Chimborazo 64.47	35.53	71.15	28.85	2.40	-6.68
	Cotopaxi	67.11	32.89	71.16	28.84	1.72	-4.04
	Ilinizas	59.78	40.22	69.55	30.45	7.35	-9.76
	Pichincha	69.90	30.10	73.77	26.23	13.44	-3.87
	Sangay	75.51	24.49	75.99	24.01	5.33	-0.48
	Tungurahua	65.75	34.25	71.92	28.08	6.99	-6.17

Acknowledgements

We thank Dr. Fabien Anthelme for revising the final form of the paper.

Funding information

GP was financed through the Fondo de Apoyo para Profesores Asistentes -FAPA of the University of the Andes [PR.3.2016.3708].

Author contributions

G.P. and A.G. conceived of the research idea; G.P., M.G. and A.G. processed the data; G.P. performed the statistical analyses and wrote the paper; all authors discussed the methods, results and commented on the manuscript.

Data accessibility

The vegetation data is freely available from VegParamo (www.vegparamo.com)

Supplementary materials

Table S1: vegetation and occurrence data used in this study (Peyre, 2020) Table S2: List of species included in the analysis, including their representation in the final dataset, mid-páramo versus super-páramo status and calculated maximum dispersal distance (in m/yr).

Tabla S3: Species composition of the Ecuadorian sky islands for 2000 and 2100 according to the modelling analyses and under the different climate change scenarios.