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

Pedro Rosero, Verónica Crespo-pérez, Rodrigo Espinosa, Patricio Andino, Álvaro Barragán, Pierre Moret, Mauro Gobbi, Gentile Francesco Ficetola, Ricardo Jaramillo, Priscilla Muriel, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Pedro Rosero, Verónica Crespo-pérez, Rodrigo Espinosa, Patricio Andino, Álvaro Barragán, et al. Multi-taxa colonisation along the foreland of a vanishing equatorial glacier. Ecography, 2021, 44 (7), pp.1010-1021. 10.1111/ecog.05478. hal-03220245

HAL Id: hal-03220245

https://hal.inrae.fr/hal-03220245

Submitted on 7 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research

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

Pedro Rosero, Verónica Crespo-Pérez, Rodrigo Espinosa, Patricio Andino, Álvaro Barragán, Pierre Moret, Mauro Gobbi, Gentile Francesco Ficetola, Ricardo Jaramillo, Priscilla Muriel, Fabien Anthelme, Dean Jacobsen, Olivier Dangles, Thomas Condom, Ludovic Gielly, Jérôme Poulenard, Antoine Rabatel, Rubén Basantes, Bolívar Cáceres Correa and Sophie Cauvy-Fraunié

P. Rosero and S. Cauvy-Fraunié (sophie.cauvy-fraunie@inrae.fr), Inrae, UR RIVERLY, Centre de Lyon-Villeurbanne, Villeurbanne Cedex, France. – V. Crespo-Pérez, R. Espinosa and P. Andino, Laboratorio de Limnología, Museo de Zoología QCAZ I, Escuela de Ciencias Biológicas, Pontificia Universidad Católica del Ecuador, Quito, Ecuador. RE also at: Grupo de Biogeografía y Ecología Espacial – BioGeoE2, Facultad de Ciencias de la Vida, Universidad Regional Amazónica Ikiam, Tena, Ecuador. – Á. Barragán, Laboratorio de Entomología, Museo de Zoología QCAZ I, Escuela de Ciencias Biológicas, Pontificia Universidad Católica del Ecuador, Quito, Ecuador. – P. Moret, Laboratoire TRACES UMR 5608, CNRS, Université de Toulouse, France. – M. Gobbi, Section of Invertebrate Zoology and Hydrobiology, MUSE-Science Museum, Corso del Lavoro e della Scienza 3, Trento, Italy. – G. F. Ficetola, Dipartimento di Scienze e Politiche Ambientali, Università degli Studi di Milano, Milano, Italy. – GFF and L. Gielly, Laboratoire d'Écologie Alpine (LECA), UMR 5553, Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, Grenoble. – R. Jaramillo and P. Muriel, Laboratorio de Ecosifisiología, Escuela de Ciencias Biológicas, Pontificia Universidad Católica del Ecuador. – F. Anthelme, AMAP, U. Montpellier, IRD, CIRAD, CNRS, INRA, Montpellier, France. – D. Jacobsen, Freshwater Biological Section, Dept of Biology, Univ. of Copenhagen, Copenhagen, Denmark. – O. Dangles, CEFE, Univ Montpellier, CNRS, EPHE, IRD, Univ Paul Valéry Montpellier 3, Montpellier, France. – T. Condom and A. Rabatel, Univ. Grenoble Alpes, IRD, CNRS, Grenoble INP, Institut des Géosciences de l'Environnement (IGE, UMR 5001), Grenoble, France. – J. Poulenard, EDYTEM, Université Savoie Mont-Blanc, Université Grenoble Alpes, CNRS, Le Bourget du Lac, France. – R. Basantes, Earth and Climate Science Research Group, Universidad Regional Amazónica Ikiam, Tena, Ecuador. – B. Cáceres Correa, Instituto Nacional de Meteorología e Hidrología (INAMHI), Quito – Ecuador.

## Ecography

44: 1–12, 2021

doi: 10.1111/ecog.05478

Subject Editor: Kenneth Feeley

Editor-in-Chief: Miguel Araújo

Accepted 24 March 2021

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

www.ecography.org

Retreating glaciers, icons of climate change, release new potential habitats for both aquatic and terrestrial organisms. High-elevation species are threatened by temperature increases and the upward migration of lowlands species. Improving our understanding of successional processes after glacier retreat becomes urgent, especially in the tropics, where glacier shrinkage is particularly fast. We examined the successional patterns of aquatic invertebrates, ground beetles, terrestrial plants, soil eukaryotes (algae, invertebrates, plants) in an equatorial glacier foreland (Carihuairazo, Ecuador). Based on both taxonomical identification and eDNA metabarcoding, we analysed the effects of both environmental conditions and age of deglacierization on community composition. Except for algae, diversity increased with time since deglacierization, especially among passive dispersers, suggesting that dispersal was a key driver structuring the glacier foreland succession. Spatial β-diversity was mainly attributed to nestedness for aquatic invertebrates, terrestrial plants and soil algae, likely linked to low environmental variability within the studied glacier foreland; and to turnover for soil invertebrates, suggesting competition exclusion at the oldest successional stage. Pioneer communities were dominated by species exhibiting flexible feeding strategies and high dispersal ability (mainly transported by wind), probably colonising from lower altitudes, or from the glacier in the case of algae. Overall, glacier foreland colonisation in the tropics exhibit common characteristics to higher latitudes. High-elevation species are nevertheless threatened, as the imminent extinction of many tropical glaciers will affect

© 2021 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos
This is an open access article under the terms of the Creative Commons
Attribution License, which permits use, distribution and reproduction in any
medium, provided the original work is properly cited.
species associated to glacier-influenced habitats but also prevent cold-adapted and hygrophilous species from using these
habitats as refuges in a warming world.

Keywords: early succession, equatorial glacier foreland, glacier retreat, multiple-taxa

**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 migra-
tion for high-elevation species (Elsen and Tingley 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 syn-
chronised. Squeezed at both ends, between the glacier fore-
front 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 gla-
cier 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, bacte-
ria), especially since the developments of high-throughput
sequencing techniques (Schmidt et al. 2008, Brown and
Jumpponen 2014, Schmidt et al. 2014). Previous stud-
ies described significant temporal variability in community
composition, with pioneer species present at early succes-
sional stages and increasing diversity, biomass and trophic
network complexity with increasing time since deglacieriza-
tion (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 com-
petitive species (Kaufmann 2001, Gobbi et al. 2007) and the
addition model (nestedness) where early colonizers persist
in older zones (Tampucci et al. 2015, Vater and Matthews
2015). Successional patterns are driven by spatial processes
(dispersal limitation), environmental filters and species inter-
action (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 simultane-
ously, at the same site, patterns of succession across different
taxa (Jumpponen et al. 2012, Lencioni and Gobbi 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 seasonal-
ity (e.g. higher ultraviolet radiation and lower oxygen con-
centration 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 inverte-
brates, ground beetles, terrestrial plants and soil eukaryotes
(algae, invertebrates, plants). We also examined the func-
tional traits of species dominating pioneer communities. For
all taxonomic groups, we expected a significant dissimilarity
in taxa assemblage along the gradient of time since deglacier-
ization, characterized by higher diversity values at the older
stages of glacier foreland. Nevertheless, we hypothesised that
the strength of the relationships between age of deglacieriza-
tion and diversity, as well as the community patterns (nest-
edness or turnover) along the glacier foreland, would differ
across the taxonomic groups, linked to their functional traits.
In particular, we expected strongest age-diversity relation-
ships for passive dispersers and turnover pattern for competi-
tive communities. Finally, we examined whether our study
highlighted remarkable successional patterns that could be
specific to the tropics.

**Material and methods**

**Study area**

This study was conducted in a small Ecuadorian glacierized
catchment of less than 1.5 km² 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² 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
<figure><img src="image_5.png" /><figcaption>Figure 1. Map of the study area showing the sampling sites for each ecological community: aquatic invertebrates, ground beetles, terrestrial plants and soil eukaryotes represented in blue, red, green and brown, respectively. The location of the Carihuairazo volcano is indicated on the map of Ecuador by a red triangle. The plot on the lower left corner illustrates decadal means of annual precipitations and temperatures during the period 1981–2018 for the study zone (Methods for details). The aerial orthophoto was taken in 2005.</figcaption></figure>

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 deglaciated 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 2001, 2015 and 2017 (Fig. 1; 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⁻¹, 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²; 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
hand without the use of magnification. No subsampling was applied. We identified macroinvertebrates under a stereo-scope (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 randomly-collected 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 (Gobbi et al. 2018). 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², 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₂HPO₄; 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⁻²), coarse particulate organic matter quantity (g m⁻²), aquatic invertebrate density (ind m⁻²) in stream, ground beetle density (ind pitfall⁻¹ day⁻¹), 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 similarity 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 within-groups 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*.

## Results

Mean (± SD) daily water temperatures between May 2015 and 2017 increased from 0.36 (± 0.21) in zone 2015 to 2.44 (± 0.50) °C in zone 1956 (based on in-stream logger data, Table 1). From zone 2015 to 1956, mean chlorophyll a

concentrations increased from 0.01 (± 0.01) to 0.52 (± 0.53) µg cm⁻² and CPOM from 0.04 (± 0.33) to 6.77 (± 2.28) g m⁻². Mean (± SD) daily soil temperatures between May and June 2015 increased from 0.38 (± 0.51) in zone 2001 to 1.12 (± 0.83) in zone LIA (based on logger data in soil). Mean (± SD) SOM measured in soil sampled varied from 0.90 (± 0.01) in zone 2 to 2.70 (± 0.59) % in zone LIA. Mean vegetation cover increased from 0.01 (± 0.02) in zone 2015 to 48.29 (± 24.88) % in zone 1956. Mean density in aquatic invertebrates varied from 56 in zone 2015 to 901 (± 557) ind m⁻² in zone 1956 and ground beetles from 0.34 (± 0.24) in zone 2015 to 3.67 (± 3.05) ind pitfall⁻¹ day⁻¹ in zone 1956.

In the deglacierized catchment, we observed 36 taxa of aquatic invertebrates, 24 of terrestrial plants and four species of ground beetles. Based on eDNA samples, we detected 36 MOTUs of algae, 162 of soil invertebrates (including annelids, arthropods, nematodes, flatworms, tardigrades and rotifers) and 46 of terrestrial plants (including 28 mosses, and 18 vascular plants). Even though not significant, we observed an increasing trend in aquatic invertebrate richness with time since deglacierization (X²=6.88, df=3, p=0.076, Kruskal–Wallis (KW) test, Fig. 2), with higher values in zone 1956. We found no effect of the time since deglacierization on ground beetle diversity (X²=5.16, df=3, p=0.161, KW), with the four species (*Bembidion andinum*, *Dyscolus oopteroides*, *Dyscolus oreas* and *Paratrechus boussingaulti*) in all zones (Fig. 3). The genus *Bembidion* was also detected in the eDNA samples, but exclusively in zone 1991. Terrestrial plant richness significantly increased with time since deglacierization (X²=28.2, df=3, p=0.001, KW; Fig. 2), with only bryophytes (and lichens) occurring close to the glacier (zone 2015), and 83% of the recorded species occurring in zone 1956. eDNA analyses showed a significant positive relationship between the time since deglacierization and MOTU richness of both invertebrates (X²=9.01, df=3, p=0.029, KW) and plants (X²=12.52, df=3, p=0.006, KW), whereas for algae we detected higher diversity at early successional stages (X²=9.66, df=3, p=0.022, KW).

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

<table><caption>Table 1. Multiple-site Sørensen β-diversity, with both turnover and nestedness components calculated for all communities.</caption><thead><tr><th>β-Diversity</th><th>Aquatic inv.</th><th>Ground beetles</th><th>Terrestrial plants</th><th>Soil algae</th><th>Soil inv.</th><th>Soil plants</th></tr></thead><tbody><tr><td>Turnover</td><td>0.12</td><td>0</td><td>0.21</td><td>0.11</td><td>0.44</td><td>0.23</td></tr><tr><td>Nestedness</td><td>0.43</td><td>0.12</td><td>0.41</td><td>0.23</td><td>0.08</td><td>0.18</td></tr><tr><td>Sørensen</td><td>0.55</td><td>0.12</td><td>0.62</td><td>0.34</td><td>0.53</td><td>0.40</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 2. α-Diversity of (a) aquatic invertebrates, (b) ground beetles, (c) terrestrial plants, (d) soil algae, (e) soil invertebrates and (f) soil plants along the gradient of deglacierized zones (2015, 2001, 1991, 1956 and LIA).</figcaption></figure>

families (e.g. Muscidae) as well as nematodes, annelids (e.g. Lumbriculidae) and water mites appeared, but chironomids remained dominant. For ground beetles, *Dyscolus oreas* was dominant in the youngest zone (zone 2015) while *Bembidion andinum* dominated from zones 2001 to the oldest zones (zones 2001, 1991 and 1956). Mosses and, to a lesser extent, lichens dominated terrestrial plant cover in all zones. Vascular plants started to recruit in zone 2001 and they reached their highest relative cover at the end of the chronosequence (LIA, approx. 25% of the total cover). Among them, the first colonisers were Asteraceae (*Xenophyllum humile*, *Senecio nivalis*, *Werneria pumila*), Poa cucullata (Poaceae) and *Cerastium floccosum* (Caryophyllaceae). Note that the genus *Cerastium* was also detected in the eDNA samples in zones 1991, 1956 and LIA. In eDNA samples, plant communities were also globally dominated by mosses, while Zygnemophyceae were the dominant algae. Invertebrate communities in soils were dominated by nematodes in all zones, mostly Chromadorea in zone 2001 and Enoplea in older zones (zones 1991, 1956 and LIA). After nematodes, tardigrades (in particular Hypsibiidae) and arthropods (mainly the ground beetle *Bembidion*) prevailed in zones 2001 and 1991, respectively, while annelids (mainly Enchytraeidae) were the second most abundant invertebrates in zones 1956 and LIA.

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).

# Discussion

## 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
<figure><img src="image_7.png" /><figcaption>Figure 3. Relative abundance of the dominant taxa of (a) aquatic invertebrates, (b) ground beetles, (c) terrestrial plants, (d) soil algae, (e) soil invertebrates and (f) soil plants along the gradient of deglacierized zones (2015, 2001, 1991, 1956 and LIA).</figcaption></figure>

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 eco-system (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
<figure><img src="image_8.png" /><figcaption>Figure 4. Non-metric multidimensional scaling ordination of (a) aquatic invertebrates, (b) ground beetles, (c) terrestrial plants, (d) soil algae, (e) soil invertebrates and (f) soil plants along the gradient of deglacierized zones (zones 2015, 2001, 1991, 1956 and LIA). 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).</figcaption></figure>

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

<table><caption>Table 2. Results of NMDS analyses, including the goodness of the NMDS fit (stress function), ANOSIM statistics and significance of the fitted environmental vectors for each community, empty cells correspond to environmental parameters excluded from the models due to high correlation with other variables, or unmeasured variables (Methods for details).</caption><thead><tr><th></th><th></th><th>Aquatic inv.</th><th>Ground beetles</th><th>Terrestrial plants</th><th>Soil algae</th><th>Soil inv.</th><th>Soil plants</th></tr></thead><tbody><tr><td>NMDS</td><td>Stress</td><td>0.203</td><td>0.066</td><td>0.147</td><td>0.134</td><td>0.107</td><td>0.139</td></tr><tr><td rowspan="2">ANOSIM</td><td>R</td><td>0.186</td><td>0.284</td><td>0.305</td><td>0.267</td><td>0.670</td><td>0.383</td></tr><tr><td>P</td><td>0.027</td><td>0.022</td><td>0.001</td><td>0.003</td><td>0.001</td><td>0.001</td></tr><tr><td rowspan="2">Distance</td><td>R<sup>2</sup></td><td>0.512</td><td>0.067</td><td>0.597</td><td>0.586</td><td>0.920</td><td>0.721</td></tr><tr><td>P</td><td>0.050</td><td>0.950</td><td>0.001</td><td>0.002</td><td>0.001</td><td>0.001</td></tr><tr><td rowspan="2">Temp.</td><td>R<sup>2</sup></td><td>0.676</td><td>0.763</td><td></td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.011</td><td>0.217</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Moisture</td><td>R<sup>2</sup></td><td></td><td>0.700</td><td></td><td></td><td></td><td></td></tr><tr><td>P</td><td></td><td>0.400</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">pH</td><td>R<sup>2</sup></td><td>0.463</td><td></td><td></td><td>0.443</td><td>0.201</td><td>0.281</td></tr><tr><td>P</td><td>0.076</td><td></td><td></td><td>0.016</td><td>0.145</td><td>0.059</td></tr><tr><td rowspan="2">CPOM</td><td>R<sup>2</sup></td><td>0.334</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.131</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Cond.</td><td>R<sup>2</sup></td><td>0.327</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.169</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Oxygen</td><td>R<sup>2</sup></td><td>0.264</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.286</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
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 and Milner 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 colembolans, 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 and Graziano 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 cuculata* (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 temperature (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⁻²). Likewise, autochthonous organic matter was extremely low in the stream (< 0.11 µg cm⁻²). 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 colembolans (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 and Fü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 and Levesque 2006, Anthelme et al. 2021). Indeed, as precipitation gradient along elevations is often inverted in the tropics (Leuschner 2000, Anthelme and Dangles 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 glacier-influenced 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.

## Author contributions

**Pedro Rosero:** Data curation (equal); Formal analysis (equal); Methodology (equal); Writing – review and editing (equal). **Verónica Crespo-Pérez:** Data curation (equal); Writing – review and editing (equal). **Rodrigo Espinosa:** Data curation (equal); Writing – review and editing (equal).

**Patricio Andino:** Data curation (equal); Writing – review and editing (supporting). **Alvaro Barragan:** Data curation (equal); Writing – review and editing (supporting). **Pierre Moret:** Data curation (equal); Investigation (equal); Writing – review and editing (equal). **Mauro Gobbi:** Investigation (equal); Writing – review and editing (equal).

**Gentile Francesco Ficetola:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review and editing (equal).

**Ricardo Jaramillo:** Data curation (equal); Writing – review and editing (equal). **Priscilla Muriel:** Writing – review and editing (supporting). **Fabien Anthelme:** Investigation (equal); Writing – review and editing (equal). **Dean Jacobsen:** Data curation (equal); Investigation (equal); Writing – review and editing (equal). **Olivier Dangles:** Investigation (equal); Writing – review and editing (equal).

**Thomas Condom:** Data curation (equal); Writing – review and editing (equal). **Ludovic Gielly:** Data curation (equal); Methodology (equal); Writing – review and editing (supporting). **Jerôme Poulenard:** Writing – review and editing (equal). **Antoine Rabatel:** Data curation (equal); Writing – review and editing (equal). **Rubén Basantes:** Data curation (equal); Writing – review and editing (equal).

**Bolívar Cáceres Correa:** Data curation (equal); Writing – review and editing (equal). **Sophie Cauvy-Fraunié:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review and editing (equal).

## Transparent Peer Review

The peer review history for this article is available at https://publons.com/publon/10.1111/ecog.05478.

## Data availability statement

Data available from the Harvard Dataverse Repository: <https://doi.org/10.7910/DVN/K3VPZW> (Rosero et al. 2021).
## References

Alexander, J. M. et al. 2018. Lags in the response of mountain plant communities to climate change. – Global Change Biol. 24: 563–579.

Anthelme, F. and Dangles, O. 2012. Plant–plant interactions in tropical alpine environments. – Perspect. Plant Ecol. Evol. Syst. 14: 363–372.

Anthelme, F. et al. 2021. Living at the edge: increasing stress for plants 2–13 years after the retreat of a tropical glacier. – Front. Ecol. 9: in press.

Baselga, A. 2010. Partitioning the turnover and nestedness components of beta diversity. – Global Ecol. Biogeogr. 19: 134–143.

Bienert, F. et al. 2012. Tracking earthworm communities from soil DNA. – Mol. Ecol. 21: 2017–2030.

Borcard, D. et al. 2011. Numerical ecology with R. – Springer.

Boyer, F. et al. 2016. obitools: a unix-inspired software package for DNA metabarcoding. – Mol. Ecol. Res. 16: 176–182.

Breen, K. and Levesque, E. 2006. Proglacial succession of biological soil crusts and vascular plants: biotic interactions in the high Arctic. – Botany 84: 1714–1731.

Brown, L. E. and Milner, A. M. 2012. Rapid loss of glacial ice reveals stream community assembly processes. – Global Change Biol. 18: 2195–2204.

Brown, S. P. and Jumpponen, A. 2014. Contrasting primary successional trajectories of fungi and bacteria in retreating glacier soils. – Mol. Ecol. 23: 481–497.

Buytaert, W. et al. 2006. Spatial and temporal rainfall variability in mountainous areas: a case study from the south Ecuadorian Andes. – J. Hydrol. 329: 413–421.

Cáceres, B. 2015. Carihuayrazo SW a little Ecuadorian glacier in the way to extinction. – AGUFM, C13B-0812.

Calderón-Sanou, I. et al. 2020. From environmental DNA sequences to ecological conclusions: how strong is the influence of methodological choices? – J. Biogeogr. 47: 193–206.

Cauvy-Fraunié, S. and Dangles, O. 2019. A global synthesis of biodiversity responses to glacier retreat. – Nat. Ecol. Evol. 3: 1675-1680.

Cauvy-Fraunié, S. and Dangles, O. 2020. Reply to: Glacial ecosystems are essential to understanding biodiversity responses to glacier retreat. – Nat. Ecol. Evol. 4: 688-689.

Cauvy-Fraunié, S. et al. 2014. Relationships between stream macroinvertebrate communities and new flood-based indices of glacial influence. – Freshwater Biol. 59: 1916–1925.

Clements, F. E. 1916. Plant succession: an analysis of the development of vegetation. – Carnegie Inst. of Washington.

Cooper, W. S. 1923. The recent ecological history of Glacier Bay, Alaska: the present vegetation cycle. – Ecology 4: 223–246.

Crespo-Pérez, V. et al. 2020. Functional structure and diversity of invertebrate communities in a glacierised catchment of the tropical Andes. – Freshwater Biol. 65: 1348-1362.

Domínguez, E. and Fernández, H. R. 2009. Macroinvertebrados bentónicos sudamericanos. Sistemática y biología. – Fundación Miguel Lillo, Tucumán, Argentina.

Dullinger, S. et al. 2012. Extinction debt of high-mountain plants under twenty-first-century climate change. – Nat. Clim. Change 2: 619.

Elsen, P. R. and Tingley, M. W. 2015. Global mountain topography and the fate of montane species under climate change. – Nat. Clim. Change 5: 772–776.

Espinosa, R. et al. 2020. Diversity patterns of aquatic macroinvertebrates in a tropical high-Andean catchment. – Rev. Biol. Trop. 68: 29–53.

Fernández-Martínez, M. A. et al. 2017. Microbial succession dynamics along glacier forefield chronosequences in Tierra del Fuego (Chile). – Polar Biol. 40: 1939–1957.

Ficetola, G. F. et al. 2015. Replication levels, false presences and the estimation of the presence/absence from eDNA metabarcoding data. – Mol. Ecol. Res. 15: 543–556.

Ficetola, G. F. et al. 2021a. Comparison of markers for the monitoring of freshwater benthic biodiversity through DNA metabarcoding. – Mol. Ecol. in press, doi: 10.1111/mec.15674.

Ficetola, G. F. et al. 2021b. Dynamics of ecological communities following glacier retreats. – Annu. Rev. Ecol. Evol. Syst. in press.

Gobbi, M. et al. 2006. Epigean arthropod succession along a 154-year glacier foreland chronosequence in the Forni Valley (central Italian Alps). – Arct. Antarct. Alp. Res. 38: 357–362.

Gobbi, M. et al. 2007. Environmental features influencing carabid beetle (Coleoptera) assemblages along a recently deglaciated area in the Alpine region. – Ecol. Entomol. 32: 682–689.

Gobbi, M. et al. 2010. Plant adaptive responses during primary succession are associated with functional adaptations in ground beetles on deglaciated terrain. – Comm. Ecol. 11: 223–231.

Gobbi, M. et al. 2018. Hand searching versus pitfall trapping: how to assess biodiversity of ground beetles (Coleoptera: Carabidae) in high altitude equatorial Andes? – J. Insect Conserv. 22: 533–543.

Gottfried, M. et al. 2012. Continent-wide response of mountain vegetation to climate change. – Nat. Clim. Change 2: 111.

Granados, H. D. et al. 2007. Chronicle of a death foretold: extinction of the small-size tropical glaciers of Popocatépetl volcano (Mexico). – Global Planet. Change 56: 13–22.

Guardiola, M. et al. 2015. Deep-sea, deep-sequencing: metabarcoding extracellular DNA from sediments of marine canyons. – PLoS One 10: e0139633.

Guerrieri, A. et al. 2021. Effects of soil preservation for biodiversity monitoring using environmental DNA. – Mol. Ecol. in press, doi: 10.1111/mec.15674

Hágvar, S. 2012. Primary succession in glacier forelands: how small animals conquer new land around melting glaciers. – In: Young, S. (ed.), International perspectives on global environmental change. InTech, pp. 151–172.

Hágvar, S. et al. 2020. Ecosystem birth near melting glaciers: a review on the pioneer role of ground-dwelling arthropods. – Insects 11: 644.

Hersbach, H. and Dee, D. 2016. ERA5 reanalysis is in production, ECMWF Newsletter 147. – ECMWF, Reading, UK. Retrieved from <https://www.ecmwf.int/en/newsletter/147/news/era5-reanalysis-production>.

Jacobsen, D. 2020. The dilemma of altitudinal shifts: caught between high temperature and low oxygen. – Front. Ecol. Environment. 18: 211–218.

Jacobsen, D. et al. 2010. Longitudinal zonation of macroinvertebrates in an Ecuadorian glacier-fed stream: do tropical glacial systems fit the temperate model? – Freshwater Biol. 55: 1234–1248.

Jacobsen, D. et al. 2014. Temporal variability in discharge and benthic macroinvertebrate assemblages in a tropical glacier-fed stream. – Freshwater Sci. 33: 32–45.

Jomelli, V. et al. 2009. Fluctuations of glaciers in the tropical Andes over the last millennium and paleoclimatic implications: a review. – Palaeogeogr. Palaeoclimatol. Palaeoecol. 281: 269–282.

Jumpponen, A. et al. 2012. Twenty years of research on fungal–plant interactions on Lyman Glacier forefront–lessons learned and questions yet unanswered. – Fungal Ecol. 5: 430–442.
Kašťovská, K. et al. 2005. Microbial assemblages in soil microbial succession after glacial retreat in Svalbard (High Arctic). – *Microbial Ecol.* 50: 396–407.

Kaufmann, R. 2001. Invertebrate succession on an alpine glacier foreland. – *Ecology* 82: 2261–2278.

Kaufmann, R. et al. 2002. The soil fauna of an alpine glacier foreland: colonization and succession. – *Arct. Antarct. Alp. Res.* 34: 242–250.

Khamis, K. et al. 2016. Glacier-groundwater stress gradients control alpine river biodiversity. – *Ecohydrology* 9: 1263–1275.

Khedim, N. et al. 2020. Topsoil organic matter build-up in glacier forelands around the world. – *Global Change Biol.* 27: 1662-1677.

Laraque, A. et al. 2007. Heterogeneous distribution of rainfall and discharge regimes in the Ecuadorian Amazon basin. – *J. Hydro-meteorol.* 8: 1364–1381.

Lee, J. R. et al. 2017. Climate change drives expansion of Antarctic ice-free habitat. – *Nature* 547: 49.

Lencioni, V. and Gobbi, M. 2018. Do carabids (Coleoptera: Carabidae) and chironomids (Diptera: Chironomidae) exhibit similar diversity and distributional patterns along a spatio-temporal gradient on a glacier foreland? – *J. Limnol.* 77: 187–195.

Lencioni, V. et al. 2015. Cold adaptive potential of chironomids overwintering in a glacial stream. – *Physiol. Entomol.* 40: 43–53.

Leuschner, C. 2000. Are high elevations in tropical mountains arid environments for plants? – *Ecology* 81: 1425–1436.

Madsen, P. B. et al. 2015. Altitudinal distribution limits of aquatic macroinvertebrates: an experimental test in a tropical alpine stream. – *Ecol. Entomol.* 40: 629–638.

Matthews, J. A. 1992. The ecology of recently-deglaciated terrain: a geocological approach to glacier forelands. – *Cambridge Univ. Press*.

Meneses, R. I. et al. 2015. La Cordillera Real y sus plantas. – *IRD; 9Herbario Nacional de Bolivia*, 2015.

Milner, A. M. 1994. Colonization and succession of invertebrate communities in a new stream in Glacier Bay National Park, Alaska. – *Freshwater Biol.* 32: 387–400.

Milner, A. M. et al. 2001. Trends of macroinvertebrate community structure in glacier-fed rivers in relation to environmental conditions: a synthesis. – *Freshwater Biol.* 46: 1833–1847.

Milner, A. M. et al. 2008. Colonization and development of an Alaskan stream community over 28 years. – *Front. Ecol. Environment.* 6: 413–419.

Milner, A. M. et al. 2011. Evolution of a stream ecosystem in recently deglaciated terrain. – *Ecology* 92: 1924–1935.

Milner, A. M. et al. 2017. Glacier shrinkage driving global changes in downstream systems. – *Proc. Natl Acad. Sci. USA* 114: 9770–9778.

Moret, P. 2005. Los coleópteros Carabidae del páramo en los Andes del Ecuador. *Sistemática, ecología y bio-geografía*. – *Quito, Pontifia Universidad Católica del Ecuador, Centro de Biodiversidad y Ambiente, Monografía* 2: 306.

Moret, P. et al. 2020. When the ice has gone: colonisation of equatorial glacier forelands by ground beetles (*Coleoptera: Carabidae*). – *Neotrop. Entomol.* 49: 213-226.

Muhlfeld, C. C. et al. 2020. Specialized meltwater biodiversity persists despite widespread deglaciation. – *Proc. Natl Acad. Sci. USA* 117: 12208–12214.

Niedrist, G. H. and Füreder, L. 2018. When the going gets tough, the tough get going: the enigma of survival strategies in harsh glacial stream environments. – *Freshwater Biol.* 63: 1260–1272.

Oehl, F. et al. 2011. Succession of arbuscular mycorrhizal communities in the foreland of the retreating Morteratsch glacier in the central Alps. – *Pedobiologia* 54: 321–331.

Pansu, J. et al. 2015. Reconstructing long-term human impacts on plant communities: an ecological approach based on lake sediment DNA. – *Mol. Ecol.* 24: 1485–1498.

Parolo, G. and Graziano, R. 2008. Upward migration of vascular plants following a climate warming trend in the Alps. – *Basic Appl. Ecol.* 9: 100–107.

Permana, D. S. et al. 2019. Disappearance of the last tropical glaciers in the Western Pacific Warm Pool (Papua, Indonesia) appears imminent. – *Proc. Natl Acad. Sci. USA* 116: 26382–26388.

Rabatel, A. et al. 2013. Current state of glaciers in the tropical Andes: a multi-century perspective on glacier evolution and climate change. – *Cryosphere* 7: 81–102.

Rabatel, A. et al. 2018. Toward an imminent extinction of Colombian glaciers? – *Geografiska Annaler Ser A Phys Geogr.* 100: 75–95.

Raso, L. et al. 2014. Intraguild predation in pioneer predator communities of alpine glacier forelands. – *Mol. Ecol.* 23: 3744–3754.

Rivadeneira Barba, G. E. 2019. Hábitos alimenticios de macroinvertebrados en dos ríos de una cuenca glaciar. – *MSc thesis*, PUCE-Quito, Equador.

Robbins, W. W. 1918. Successions of vegetation in Boulder Park, Colorado. – *Bot. Gazette* 65: 493–525.

Rosero, P. et al. 2021. Data from: Multi-taxa colonisation along the foreland of a vanishing equatorial glacier. – *Harvard Dataverse Repository*: <https://doi.org/10.7910/DVN/K3VPZW>.

Schmidt, S. et al. 2008. The earliest stages of ecosystem succession in high-elevation (5000 metres above sea level): recently deglaciated soils. – *Proc. R. Soc. B* 275: 2793–2802.

Schmidt, S. et al. 2014. Do bacterial and fungal communities assemble differently during primary succession? – *Mol. Ecol.* 23: 254–258.

Simmons, B. L. et al. 2009. Terrestrial mesofauna in above- and below-ground habitats: Taylor Valley, Antarctica. – *Polar Biol.* 32: 1549–1558.

Stibal, M. et al. 2020. Glacial ecosystems are essential to understanding biodiversity responses to glacier retreat. – *Nat. Ecol. Evol.* 4: 686–687.

Taberlet, P. et al. 2012. Towards next-generation biodiversity assessment using DNA metabarcoding. – *Mol. Ecol.* 21: 2045–2050.

Taberlet, P. et al. 2018. Environmental DNA: for biodiversity research and monitoring. – *Oxford Univ. Press*.

Tampucci, D. et al. 2015. Plant and arthropod colonisation of a glacier foreland in a peripheral mountain range. – *Biodiversity 16*: 213–223.

Uehlinger, U. et al. 2010. The physico-chemical habitat template for periphyton in alpine glacial streams under a changing climate. – *Hydrobiologia* 657: 107–121.

Urban, M. C. 2018. Escalator to extinction. – *Proc. Natl Acad. Sci. USA* 115: 11871–11873.

Vater, A. E. and Matthews, J. A. 2015. Succession of pitfall-trapped insects and arachnids on eight Norwegian glacier forelands along an altitudinal gradient: patterns and models. – *Holocene* 25: 108–129.

Wilhelm, L. et al. 2013. Microbial biodiversity in glacier-fed streams. – *ISME J.* 7: 1651–1660.

Zemp, M. et al. 2019. Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. – *Nature* 568: 382.

Zimmer, A. et al. 2018. Time lag between glacial retreat and upward migration alters tropical alpine communities. – *Perspect. Plant Ecol. Evol. Syst.* 30: 89–102.