Multi-taxa colonisation along the foreland of a vanishing equatorial glacier

species associated to glacier-influenced habitats but also prevent cold-adapted and hygrophilous species from using these habitats as refuges in a warming world.


Introduction
One major concern in alpine ecology is to determine the time lag between climatic shifts and the upward migration of species to adapt to the effects of warming (Dullinger et al. 2012, Gottfried et al. 2012, Alexander et al. 2018, and to assess whether high-elevation organisms would have time to colonise and establish in higher zones before the arrival of species migrating from lowlands. While the area of suitable habitat available for colonization declines with increasing elevation, thereby restricting the possibility of upward migration for high-elevation species (Elsen andTingley 2015, Urban 2018), glacier retreat releases new potential habitats for both aquatic and terrestrial species (Breen and Levesque 2006, Milner et al. 2011, Lee et al. 2017. Worldwide, mountain glaciers have been shrinking since the end of the little ice age (LIA), about 150-250 years ago, and shrinkage rates have accelerated dramatically over the last five decades (Zemp et al. 2019). However, rates of glacier shrinkage and upward migration of high-elevation species might not be synchronised. Squeezed at both ends, between the glacier forefront and the front of invaders, high-elevation species may thus be particularly vulnerable to global warming. Studies on early succession after glacier retreat are therefore crucial.
Research on early succession in glacier forelands is mostly based on biological observations along a chronosequence of zones with known deglacierization dates. Succession in glacier forelands has a long research history, with earlier studies focusing on plants (Clements 1916, Robbins 1918, Cooper 1923, Matthews 1992, later on aquatic (Milner 1994, Milner et al. 2008) and terrestrial invertebrates (Kaufmann 2001, Gobbi et al. 2006, Vater and Matthews 2015, and more recently on microbial communities (e.g. fungi, bacteria), especially since the developments of high-throughput sequencing techniques (Schmidt et al. 2008, Brown and Jumpponen 2014, Schmidt et al. 2014. Previous studies described significant temporal variability in community composition, with pioneer species present at early successional stages and increasing diversity, biomass and trophic network complexity with increasing time since deglacierization (Matthews 1992, Raso et al. 2014, Ficetola et al. 2021b). Among the studied glacial successions, two different patterns have been observed: the replacement model (turnover) where early colonizers are replaced by more tolerant and/or competitive species (Kaufmann 2001, Gobbi et al. 2007) and the addition model (nestedness) where early colonizers persist in older zones (Tampucci et al. 2015, Vater andMatthews 2015). Successional patterns are driven by spatial processes (dispersal limitation), environmental filters and species interaction (facilitation versus competition; Ficetola et al. 2021b). However, the relative importance of these mechanisms is poorly understood and varies according to the geographical zone and the taxonomical group considered (Matthews 1992, Brown and Jumpponen 2014, Vater and Matthews 2015. Although quite rare, it is important to compare simultaneously, at the same site, patterns of succession across different taxa (Jumpponen et al. 2012, Lencioni andGobbi 2018) to disentangle the ecological processes from the geographical influences (e.g. geology, geomorphology, climate, altitude). In particular, in the tropics, colonization patterns may differ compared to higher latitudes, where most ecological studies on glacier forelands were conducted (Ficetola et al. 2021b), due to differences in environmental conditions and seasonality (e.g. higher ultraviolet radiation and lower oxygen concentration at high elevations, no permanent or seasonal snow cover; Madsen et al. 2015, Jacobsen 2020. We applied a multi-taxa approach to identify the processes influencing the successional patterns after glacier retreat in a tropical glacier foreland (Carihuairazo, Ecuador) along a gradient of time since deglacierization from the end of the LIA to the present. For this, we compared the successional patterns across various taxonomical groups: aquatic invertebrates, ground beetles, terrestrial plants and soil eukaryotes (algae, invertebrates, plants). We also examined the functional traits of species dominating pioneer communities. For all taxonomic groups, we expected a significant dissimilarity in taxa assemblage along the gradient of time since deglacierization, characterized by higher diversity values at the older stages of glacier foreland. Nevertheless, we hypothesised that the strength of the relationships between age of deglacierization and diversity, as well as the community patterns (nestedness or turnover) along the glacier foreland, would differ across the taxonomic groups, linked to their functional traits. In particular, we expected strongest age-diversity relationships for passive dispersers and turnover pattern for competitive communities. Finally, we examined whether our study highlighted remarkable successional patterns that could be specific to the tropics.

Study area
This study was conducted in a small Ecuadorian glacierized catchment of less than 1.5 km 2 at the Carihuairazo volcano (5018 m a.s.l., 01°24′25″S, 78°45′00″W, Fig. 1). The only remaining glacier on this mountain, with an area of about 0.015 km 2 in 2017, has lost around 95% of its surface area since 1956, and will probably disappear in the early 2020s, due to its very small size and relatively low elevation for a tropical glacier (Rabatel et al. 2013). The study catchment included one glacier-fed stream, originating at the glacier snout. The study was conducted in five different zones delimited by dated glacier snout position: these zones were deglacierized between 2015 and 2017, 2001 and 2005, 1991 and 2001, 1956 and 1991, and between the LIA maximum and 1956 (although not locally precisely dated, the LIA could likely be estimated from the early 18th to early 19th centuries (Jomelli et al. 2009)). Hereafter, zones were referred to zones 2015,2001,1991,1956 and LIA. Glacier outlines were delineated manually based on aerial photographs from 1956 and 2005, LANDSAT image from 1991, and directly measured in the field from topographical measurements (differential-GNSS) in Cáceres 2015). Study zones were located between 4770 and 4630 m a.s.l., at distances from the glacier front location in 2017 ranging from 25 to 900 m. Climatic data of the period 1981-2018, extracted from the reanalysis ERA5-Land (spatial resolution 0.1°, monthly temporal resolution) provided by the ECMWF (Hersbach and Dee 2016), showed, at the decadal scale, that mean annual precipitation is slightly decreasing, from 570 to 520 mm y −1 , while mean annual temperature is slightly increasing, from 5.3 to 5.6°C ( Fig. 1). At the Carihuairazo foreland, six distinct ecological communities were characterized: aquatic invertebrates, ground beetles, terrestrial plants at zones 2015, 2001, 1991 and 1956, and soil eukaryotes: soil algae, soil invertebrates and soil plants at zones 2001, 1991, 1956 and LIA.

Aquatic invertebrates
At zones 2015, 2001, 1991 and 1956, we obtained six temporal replicates of benthic macroinvertebrate samples in the glacier-fed stream. At each sampling zone and date (six times between May 2015 and May 2017 for zones 2001, 1991 and 1956, but only in May 2017 for zone 2015, which was previously covered by ice), we randomly collected five Surber samples (0.05 m 2 ; mesh size 200 µm) along a 25 m reach from pebble-cobble substratum. Samples were preserved in 70% ethanol and taken to the laboratory, where they were rinsed through a 200-μm sieve and sorted thoroughly by 4 hand without the use of magnification. No subsampling was applied. We identified macroinvertebrates under a stereoscope (at 7-90× magnification) mostly to family level and separated into morphospecies, according to (Domínguez and Fernández 2009). At each sampling zone and date, we measured temperature, conductivity (at 25°C), pH, oxygen concentration and saturation, with portable meters (model Cond 315i, pH 315i and YSI 55, WTW, Weilheim, Germany). We estimated food resources available to macroinvertebrates by quantifying chlorophyll a in epilithic algae on nine randomlycollected pebbles and by measuring benthic coarse particulate organic matter from Surber samples (Jacobsen et al. 2014 for details). At each zone, during the entire study period we recorded water temperature at a 10 min time-step, using Hobo temperature loggers (Onset Computer Corp, USA).

Ground beetles (Coleoptera: Carabidae)
At zones 2015, 2001, 1991 and 1956, we obtained three temporal replicates of ground beetle samples. At each sampling zone and date (three times between May 2015 and May 2017), we randomly deployed six pitfall traps about 10 m apart from each other during a period of 13 days (Moret et al. 2020 for more details). Traps consisted of plastic vessels (diameter 7 cm, height 10 cm) baited with a mixture of wine-vinegar and salt . Ground beetles were preserved in 70% ethanol and identified to species level, following Moret (2005). At each zone, we recorded soil temperature and moisture (at ~2 cm depth) at 10 min time-steps during the sampling periods using Hobo temperature and humidity loggers (Onset Computer Corp, USA). Although other families of arthropods were collected in pitfall traps, we only focused on ground beetles, because they represent a major component of the meso-and macro ground-dwelling fauna on recently deglacierized terrains in terms of abundance, and because their taxonomy and ecology are fairly well known (Moret et al. 2020).

Terrestrial plants
At zones 2015, 2001, 1991 and 1956, we performed ten spatial replicates of terrestrial plant observation in May 2017. Terrestrial plant coverage was recorded in ten squared-plots of 1 m 2 , a surface that proved to be representative of the plant diversity in this type of ecosystems (Zimmer et al. 2018). At each plot, we listed all observed vascular plants and visually estimated their relative cover (%), as well as the relative cover of lichens and mosses.

Soil eukaryotes
At zones 2001,1991,1956 and LIA, we obtained five spatial replicates of soil samples in May 2015. At each sampling zone, we collected soil samples (15 g) about 25 m apart from each other, using sterilized shovels. We preserved soil in sterilized tubes containing 30 g of silicagel, that allows long-term preservation of eDNA for metabarcoding studies (Guerrieri et al. 2021). In the laboratory, extra-cellular environmental DNA was extracted from soil samples following Bienert et al. (2012). For each sample, we mixed ~15 g of soil with 15 ml of saturated phosphate buffer (Na 2 HPO 4 ; 0.12 M; pH ≈ 8) during 15 min. Then, 2 ml of the mixture were centrifuged (10 min, 10 000 g) and 400 µl of the resulting supernatant were kept as starting material for DNA extraction, using the NucleoSpin® Soil kit (Macherey-Nagel; Taberlet et al. 2012). Eukaryote DNA was amplified with the Euka02 primer, adapted to eukaryote taxa (Guardiola et al. 2015, Taberlet et al. 2018. This marker is able to amplify most of eukaryotes with very limited bias, but with a limited taxonomic resolution, particularly for some taxa (e.g. within vascular plants; Taberlet et al. 2018, Ficetola et al. 2021a). All samples and controls underwent four PCR replicates (Ficetola et al. 2015). Sequencing was performed on Illumina MiSeq platform. We filtered DNA sequences using the OBITOOLS software (Boyer et al. 2016), as described in Pansu et al. (2015). Using the ecotag program, we assigned sequences to relevant taxa by comparing them with a global database, generated from EMBL (Boyer et al. 2016). At each sample, we recorded soil temperature (at ~2 cm depth) at 10 min time-steps during one month, using Hobo temperature loggers (Onset Computer Corp, USA). An aliquot of soil samples was dried and 2 mm-sieved. We measured soil pH following the ISO 10390:2005 norm using a pH-meter (pH7110, inoLab) in a 1/5 solution (1 volumic part soil sample and 4 volumic parts distilled water), and estimated soil organic matter (SOM) content by loss on ignition (after 4 h incubation in a muffle furnace at 550°C).

Data analyses
For each zone, community composition was characterised based on six temporal replicates for aquatic invertebrates, three temporal replicates for ground beetles, ten spatial replicates for terrestrial plants, and five spatial replicates for soil eukaryotes. Based on these replicates, we first calculated the average of chlorophyll a concentration (µg cm −2 ), coarse particulate organic matter quantity (g m −2 ), aquatic invertebrate density (ind m −2 ) in stream, ground beetle density (ind pitfall −1 day −1 ), terrestrial vegetation cover (%) and soil organic matters (%). We then calculated α-diversity of 1) aquatic invertebrates, as the total number of taxa collected in the five Surber samples at each sampling date, 2) ground beetles, as the total number of species collected in the six pitfall traps at each sampling date, 3) terrestrial plants as the number of taxa at each plot and 4) soil algae, invertebrates and plants as the number of operational taxonomic units (MOTUs) measured for each taxonomic group at each soil-eukaryote sample. For metabarcoding data, considering all detected MOTUs may overestimate actual diversity, but allows successful identification of relationships with environmental variation (Calderón-Sanou et al. 2020). Similarly, for each zone and taxonomic group, we calculated the relative abundance of dominant taxa. For metabarcoding data, the proportion of reads of dominant taxa calculated within the three taxonomical groups (soil algae, invertebrates and plants) was used as a proxy of the relative abundance even though it provides an imperfect measure (e.g. variability in the number of reads linked to match with the primers, length of the amplified sequence). To compare the local diversity of communities across zones, we performed a Kruskal-Wallis test, followed by a Dunn post-hoc test. To characterise spatial patterns in taxon assemblage along the gradient of deglacierization, we calculated the Sørensen-based multiple-site dissimilarity index (β-diversity based on occurrence matrices) as well as both turnover (species replacement) and nestedness (species gain) components following Baselga (2010). To examine patterns of dissimilarity in community composition among deglacierized zones, we performed non-metric multidimensional scaling (NMDS) analyses based on Bray-Curtis distances, calculated on relative abundance matrices. The NMDS goodness of fit is estimated with a stress function (which ranges from 0 to 1), with values close to zero indicating a good fit. Differences in community composition among zones were tested using an analysis of similarities (ANOSIM), which tested the null hypothesis that withingroups similarity was equal to between-groups similarity. The degree of separation between zones was estimated with a statistical parameter R (which ranges from 0 to 1), with a score of 1 indicating complete separation and a score of 0 indicating no separation. Significance of the R statistic was determined by comparing the observed R value with the null distribution of the R statistic obtained with 999 permutations of group membership. To determine whether the community structure differed according to environmental variables, we fitted environmental vectors onto the NMDS ordination space (Borcard et al. 2011). Significance of the fitted environmental vectors was assessed with randomisation tests (999 permutations). To reduce collinearity, prior to the analysis we removed strongly correlated environmental variables (Pearson r > 0.8, p < 0.05); for example, only the distance to glacier was retained (the main variable used in a recent meta-analysis on this topic, Cauvy-Fraunié and Dangles 2019) among altitude, distance to glacier front and time since deglacierization. All analyses and figures were performed in R (<www.r-project.org>, ver. 3.6.2) using packages ggplot2, funrar, vegan, cowplot, PMCMR.
Dissimilarity in taxon assemblage (β-diversity) was mainly attributed to the nestedness component for aquatic invertebrates, terrestrial plants and soil algae, and to the turnover component for soil invertebrates (Table 1). Both nestedness and turnover components contributed equally to soil plant β-diversity, while β-diversity was very low for ground beetles.
We found significant relationships between the successional stages and community composition for all studied taxonomic groups (ANOSIM, p < 0.05 in all cases, Fig. 4, Table 2). Aquatic invertebrate communities were exclusively composed by chironomids (mainly Podonominae) in the youngest zone (zone 2015). From zone 2001 on, more insect Fit of the environmental variable vectors onto the NMDS ordination space showed that community structure was mainly related to distance to the glacier (R 2 = 0.512, p = 0.05) and water temperature (R 2 = 0.68, p = 0.01) for aquatic invertebrates and to distance to the glacier (R 2 = 0.60, p = 0.001) for terrestrial plants. We did not find significant relationships between environmental variables and ground beetle community composition. For soil eukaryotes, community structure was mainly related to distance to the glacier (R 2 = 0.59, p = 0.002) and pH (R 2 = 0.44, p = 0.016) for algae, and to distance to the glacier for both invertebrates (R 2 = 0.92, p = 0.001) and plants (R 2 = 0.72, p = 0.001, Table 2).

Multiple-taxa succession patterns
We compared successional patterns of terrestrial and aquatic communities along a gradient of time since deglacierization. For each taxonomic group considered, we observed significant dissimilarities in community composition along the successional stages of deglacierization. As globally observed in most previous studied in glacier forelands (Cauvy-Fraunié and Dangles 2019), successional patterns of the majority of  ( , 1991( , 1956 taxonomic groups was characterised by higher abundances (density or biomass) and diversities at the oldest zones, indicating that colonisation occurred from lower altitudes. On the contrary, algal richness was higher close to the glacier, suggesting potential algae colonisation from the glacier ecosystem (Kaštovská et al. 2005, Stibal et al. 2020).
Our results showed slight variation in diversity pattern across taxonomic group. First, the strength of the relationship between taxon richness and time since deglacierization differed among the groups of organisms considered (Fig. 2). Indeed, increase in diversity was more pronounced for passive dispersers, i.e. terrestrial plants, suggesting that dispersal was a key mechanism structuring succession after glacier retreat (Zimmer et al. 2018). Second, spatial β-diversity was not necessarily dominated by the same component throughout the studied communities. In accordance with the addition and persistence model, spatial β-diversity in ground beetles resulted from a nestedness pattern in equatorial glacierized catchments (Moret et al. 2020). Similarly, spatial β-diversity in terrestrial plants and in aquatic invertebrates was mainly linked to the nestedness component, whereas turnover pattern usually prevailed in glacier-fed streams and have been often related to longitudinal environmental changes (Milner et al. 2001, Cauvy-Fraunié et al. 2014, Khamis et al. 2016. In particular, the longitudinal pattern in taxa assemblage along glacier-fed streams has been linked to increase in water temperature and channel stability (Milner et al. 2001, Jacobsen et al. 2010. This dominance of the nestedness patterns suggested that only a small fraction of the local species pool, from lower elevations, was able to colonise the first stages of succession. This also highlighted the small gradient of time since deglacierization within the zones from 2015 to 1956, and associated low environmental variability, preventing late-coloniser species to establish at the end of the   ( , 1991( , 1956. gradient (oldest zones). On the contrary, spatial β-diversity in soil invertebrates across successional stages was due to turnover processes suggesting competition exclusion at the oldest successional stage (Kaufmann et al. 2002).
Overall, species successions described along the Carihuairazo glacier foreland (Fig. 3) were consistent with previous observations in glacierized catchments worldwide ( Fig. 4 in Cauvy-Fraunié and Dangles 2019). Indeed, aquatic  (zones 2015(zones , 2001(zones , 1991(zones , 1956. Vectors of all tested environmental variables were overlaid onto the ordination space. The arrow indicates the direction of most rapid change in the variable, and its length is proportional to the correlation between ordination and the variable (Table 2 for statistics). invertebrate assemblages were composed by chironomids at the early stages of deglacierization, followed by other insect larvae, annelids, and, finally, nematodes, at older stages (Milner et al. 2001, Brown andMilner 2012). Microscopic invertebrates such as tardigrades, rotifers and nematodes prevailed in the soil invertebrate community close to the glacier, followed by larger arthropods (beetles such as carabids and spiders), and annelids at the oldest zones (Simmons et al. 2009). Finally, excluding both bacteria and fungi (not analysed in this study), the youngest soils were first colonised by algae, lichen and moss, while vascular plants appeared later (Fernández-Martínez et al. 2017, Hågvar et al. 2020.

Common functional traits of pioneer species
Close to the Carihuairazo glacier, observed pioneer communities were mostly dominated by species with high dispersal ability, mainly through aerial transport accentuated by wind. For example, chironomids, characterised by a winged adult stage with small body size (Milner et al. 2011), were the dominant aquatic larvae. Similarly, within the soil communities, microalgae and microfauna, such as jumping collembolans, are organisms easily transported by air currents (Hågvar 2012, Hågvar et al. 2020, also referred to as aeroplankton (Gobbi et al. 2010). Likewise, plants colonised the new deglacierized soils through seed and spore wind dispersal (Oehl et al. 2011), resulting in an overrepresentation of anemochorous species close to the glacier (Parolo andGraziano 2008, Anthelme et al. 2021). For example, the most frequent vascular species in the pioneer stages of the early succession at the Carihuairazo were the anemochorous Poa cucullata (Poaceae), Xenophyllum humile and Werneria pumila (Asteraceae), as well as Cerastium floccosum (Caryophyllaceae) dispersed by both wind and animals (Meneses et al. 2015).
On the contrary, observed ground beetles were all brachypterous (flightless), thus with low dispersal ability. Nevertheless, these species could be considered as good colonizers due to their foraging activity leading them to move as far as 10 m per night (PM pers. obs.). Among all environmental conditions, temperature was the variable most strongly related to the structure of the aquatic communities (Fig. 4, Table 2). Indeed, species occurring close to glaciers are usually well adapted to cold temperatures, with behavioural and/or physiological capability to survive freezing and grow under permanently cold temperatures (Cauvy-Fraunié and Dangles 2020). For example, certain chironomids synthesise cryoprotectants such as sugars to avoid freezing (Lencioni et al. 2015). In contrast, ground beetle assemblages were not affected by local environmental conditions (Moret 2005), and showed no specific micro habitat or microclimatic preferences, unlike their temperate (Tampucci et al. 2015) and polar counterparts (Vater and Matthews 2015). Concerning plants, the high relative cover of lichens and mosses was also probably enhanced by high levels of precipitation at Carihuairazo, partly due to the orographic control on atmospheric processes (Buytaert et al. 2006, Laraque et al. 2007).
Glacier forelands usually are oligotrophic systems with limited nutrient availability and reduced primary production (Cauvy-Fraunié and Dangles 2019, Khedim et al. 2020). Close to the glacier (zones 2015 and 2001), we observed very low organic matter content in soil (< 0.9%) and low percentage of terrestrial plant cover (< 3.8%), leading to low availability of allochthonous organic matter in the stream (< 1.67 g m −2 ). Likewise, autochthonous organic matter was extremely low in the stream (< 0.11 µg cm −2 ). This explained why recorded animals in glacier forelands (both aquatic and terrestrial) were not herbivorous but mainly detritivores and/or predators feeding on allochthonous matter such as wind-blown debris and organisms from adjacent ecosystems (Cauvy-Fraunié and Dangles 2020) or microorganisms originating from the glacier and transported by meltwater (Wilhelm et al. 2013). Observed ground beetles were indeed active predators, probably mainly feeding on collembolans (high densities observed in pitfall traps) and chironomids as observed for high-elevation Alpine carabids (Raso et al. 2014). In addition, a complementary study at the same glacier-fed stream recorded an omnivorous food diet, composed of microphytes, coarse and fine particulate matter for two chironomid taxa (Rivadeneira Barba 2019).
All these processes (dispersion, environmental evolution, soil development) are directly or indirectly related to time since glacier retreat (Ficetola et al. 2021b). Experimental studies with transplantation and/or modification of environmental and soil proprieties would allow distinguishing the relative contribution in dispersal processes, environmental filtering and biotic interactions in driving the glacier foreland succession processes. Nevertheless, our study design highlighted that dispersal limitation and nutrient lack were primary drivers.

Glacier retreat and species sensitivity in the tropics
Successional patterns observed at the Carihuairazo glacier foreland were similar to those observed at higher latitudes, i.e. globally characterised by a positive effect of time since deglacierization on diversity, but with differences among the taxonomical groups considered (Cauvy-Fraunié and Dangles 2019). All pioneer communities were mainly composed of environmental specialists (adapted to the harsh glacial habitat, especially to the permanently cold temperatures), feeder generalists (exhibiting flexible feeding strategies, Niedrist andFüreder 2018, Crespo-Pérez et al. 2020), with high dispersal capacity (Oehl et al. 2011, Hågvar 2012. However, although not tested in this study, colonisation rates might be higher in the tropics because dispersal and establishment could occur all year round compared to the limited windows of opportunities at higher latitudes, mainly due to both seasonal snow cover and seasonality in the life cycles of temperate and polar organisms (Uehlinger et al. 2010).
The rapid shrinkage of tropical glaciers induces an accelerated reduction of glacial influence in adjacent ecosystems, and may have considerable consequences on biodiversity. The decrease in meltwater supply reduces aquatic habitat availability for aquatic species (Milner et al. 2017), and increases the risk of desiccation for terrestrial habitats (Breen andLevesque 2006, Anthelme et al. 2021). Indeed, as precipitation gradient along elevations is often inverted in the tropics (Leuschner 2000, Anthelme andDangles 2012), tropical alpine plants rely significantly on ground surface water sourcing from glacier melting. In particular, plant aggregation has been frequently observed along glacier-fed stream corridors. Similarly, two hygrophilous ground beetle species (Paratrechus boussingaulti and Dyscolus oopteroides) were observed close to the perennial glacier-fed stream at Carihuairazo but were not detected along the intermittent glacier-fed stream at the rain-shadow south-west slope of the neighbouring Chimborazo volcano (unpublished hand searching survey over the same study period). In addition, the modification in environmental conditions (e.g. turbidity in the stream) linked to glacial influence will discriminate specialised species adapted to these specific habitats (Milner et al. 2008). For example, in agreement with Espinosa et al. (2020) at the Antisana volcano, certain aquatic larvae from the beetle families Dytiscidae and Staphylinidae were exclusively observed in the studied glacier-fed stream, and not in a neighbouring groundwater-fed stream (unpublished observations from the same sampling dates).
In the tropics, accelerated glacier shrinkage exposes new aquatic and terrestrial habitats, quickly colonised by high-elevation species. Nevertheless, the imminent extinction of many tropical glaciers (Granados et al. 2007, Rabatel et al. 2018, Permana et al. 2019) will first provoke the loss of glacierinfluenced habitats and associated species (Cauvy-Fraunié and Dangles 2020, Stibal et al. 2020), while also preventing cold-adapted and/or hygrophilous species from using these habitats as refuges in a warmer world. These effects can be particularly strong in the tropics where the lack of persistent snowfields precludes cold habitats as well as a supplementary seasonal water supply after the complete disappearance of glaciers contrary to high-latitude regions (Muhlfeld et al. 2020).
Acknowledgments -We warmly thank Leonardo Punina for his technical support in the field. Part of the glaciological data was provided by the Service National d'Observation GLACIOCLIM (<http://glacioclim.osug.fr/>) funded by CNRS-INSU, IRD, Université Grenoble Alpes and GREAT ICE. We thank the two anonymous reviewers for their constructive suggestions on a previous version of the manuscript. Funding -Part of this work was funded and conducted in the framework of the International Joint Laboratory GREAT-ICE, a joint initiative of the IRD, universities and institutions in Bolivia, Peru, Ecuador and Colombia. Part of the fieldwork was funded by Pontificia Universidad del Ecuador, Grant Agreement no. M13434 (Efecto del rápido retroceso glaciar sobre la biodiversidad en ecosistemas tropicales de altura). GFF and LG were funded by the European Research Council under the European Community's Horizon 2020 Programme, Grant Agreement no. 772284 (IceCommunities). Co-authors from University of Grenoble acknowledge the support of LabEx OSUG@2020 (Investissements d'avenir -ANR10 LABX56). This study is part of the project Life Without Ice, funded by the BNP Paribas Foundation. The funder had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. Conflict of interest -The authors declare to have no conflict of interest. Permits -Organisms were collected and analysed under permits: MAE-DNB-CM-2015-0030 and 002-16-IC-FAU-DNB/MA, granted by the Ecuadorian Ministry of the Environment.