# Living at the Edge: Increasing Stress for Plants 2–13 Years After the Retreat of a Tropical Glacier

Fabien Anthelme, Sophie Cauvy-Fraunié, Bernard Francou, Bolívar Cáceres, Olivier Dangles

► To cite this version:

Fabien Anthelme, Sophie Cauvy-Fraunié, Bernard Francou, Bolívar Cáceres, Olivier Dangles. Living at the Edge: Increasing Stress for Plants 2–13 Years After the Retreat of a Tropical Glacier. Frontiers in Ecology and Evolution, 2021, 9, 10.3389/fevo.2021.584872. doi:10.3128706

HAL Id: hal-03128706

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

Submitted on 2 Feb 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Living at the Edge: Increasing Stress for Plants 2–13 Years After the Retreat of a Tropical Glacier

Fabien Anthelme<sup>1,2*</sup>, Sophie Cauvy-Fraunié<sup>2,3</sup>, Bernard Francou<sup>4</sup>, Bolívar Cáceres<sup>5</sup> and Olivier Dangles<sup>2,6</sup>

<sup>1</sup> Laboratory AMAP, IRD, University of Montpellier, CIRAD, CNRS, INRA, Montpellier, France, <sup>2</sup> Facultad de Ciencias Exactas y Naturales, Pontificia Universidad Católica del Ecuador, Quito, Ecuador, <sup>3</sup> INRAE, UR RiverLy, Centre de Lyon-Villeurbanne, Villeurbanne, France, <sup>4</sup> Laboratory IGE, IRD, Université Grenoble Alpes, CNRS, G-INP, Institut des Géosciences de l'Environnement, Quito, Ecuador, <sup>5</sup> Instituto Nacional de Meteorología e Hidrología (INAMHI), Quito, Ecuador, <sup>6</sup> Centre d'Ecologie Fonctionnelle et Evolutive, UMR 5175, CNRS, Université de Montpellier, Université Paul Valéry Montpellier, EPHE, IRD, Montpellier, France

Rapid warming is a major threat for the alpine biodiversity but, at the same time, accelerated glacial retreat constitutes an opportunity for taxa and communities to escape range contraction or extinction. We explored the first steps of plant primary succession after accelerated glacial retreat under the assumption that the first few years are critical for the success of plant establishment. To this end, we examined plant succession along a very short post-glacial chronosequence in the tropical Andes of Ecuador (2–13 years after glacial retreat). We recorded the location of all plant individuals within an area of 4,200 m² divided into plots of 1 m². This sampling made it possible to measure the responses of the microenvironment, plant diversity and plants traits to time since the glacial retreat. It also made it possible to produce species-area curves and to estimate positive interactions between species. Decreases in soil temperature, soil moisture, and soil macronutrients revealed increasing abiotic stress for plants between two and 13 years after glacial retreat. This increasing stress seemingly explained the lack of positive correlation between plant diversity and time since the glacial retreat. It might explain the decreasing performance of plants at both the population (lower plant height) and the community levels (lower species richness and lower accumulation of species per area). Meanwhile, infrequent spatial associations among plants indicated a facilitation deficit and animal-dispersed plants were almost absent. Although the presence of 21 species on such a small sampled area seven years after glacial retreat could look like a colonization success in the first place, the increasing abiotic stress may partly erase this success, reducing species richness to 13 species after 13 years and increasing the frequency of patches without vegetation. This fine-grain distribution study sheds new light on nature's responses to the effects of climate change in cold biomes, suggesting that faster glacial retreat would not necessarily result in accelerated plant colonization. Results are exploratory and require site replications for generalization.

**Keywords:** abiotic stress, dispersal limitation, facilitation, glacial retreat, chronosequence, primary succession, tropical andes, plant height

## OPEN ACCESS

**Edited by:**
Luis Daniel Llambi,
Universidad de Los Andes, Venezuela

**Reviewed by:**
Esteban Suarez,
Universidad San Francisco de
Quito, Ecuador
Gianalberto Losapio,
Stanford University, United States

**\*Correspondence:**
Fabien Anthelme
Fabien.anthelme@ird.fr

**Specialty section:**
This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

**Received:** 20 July 2020
**Accepted:** 12 January 2021
**Published:** 02 February 2021

**Citation:**
Anthelme F, Cauvy-Fraunié S, Francou B, Cáceres B and Dangles O (2021) *Living at the Edge: Increasing Stress for Plants 2–13 Years After the Retreat of a Tropical Glacier*. *Front. Ecol. Evol.* 9:584872.
doi: 10.3389/fevo.2021.584872
# INTRODUCTION

The spatial distribution of plant communities along dated chronosequences after glacial retreat since the period of the little ice age (LIA) have been extensively used to characterize plant primary succession in temperate-to-arctic environments (Matthews, 1992; Walker and Del Moral, 2003; Cauvy-Fraunié and Dangles, 2019). Although this type of primary succession depends on a high number of local, regional and latitudinal parameters (D'Amico et al., 2017; Prach and Walker, 2020) as well as stochastic effects (Marteinsdóttir et al., 2013), structured patterns of plant diversity and plant interactions have been consistently correlated with time since the glacial retreat, including increased soil development and trophic interactions (Losapio et al., 2015; Cauvy-Fraunié and Dangles, 2019), increasing non-trophic interactions between plants (facilitation and competition: Fastie, 1995; Zimmer et al., 2018) and increasing plant cover and richness, as observed from field observations (Matthews, 1992; Schumann et al., 2016) and from satellite imagery analyses (Fischer et al., 2019). This goes along with the assumption that the temporal extent (and also the spatial extent) of a habitat will be correlated with the diversity of the species composing it (Tilman, 1994).

Under the effects of accelerated warming for half a century, the pace of glacier retreat has dramatically increased across the globe (Zemp et al., 2019). Large areas become rapidly available for plant primary succession. These landscapes constitute remarkable opportunities for the alpine biodiversity, especially plants, to enlarge their leading-edge distributions upward whereas their trailing-edge distributions at lower elevation mostly retract because of the combined effects of thermophilization and increased interspecific competition with up-migrating lowland species (Heckmann et al., 2015; Lenoir and Svenning, 2015). At the same time, however, these new habitats are characterized by a very short temporal extent, which is expected to impact negatively their biodiversity (Tilman, 1994). It remains unclear, however, to what extent the recent accelerated warming may affect the general patterns of plant primary succession previously observed at the scale of the LIA (Dullinger et al., 2012). Indeed, most studies on plant primary succession after glacial retreat have been carried out along chronosequences longer than 50 years, and often including the little ice age (Cauvy-Fraunié and Dangles, 2019). We expect that the recent accelerated warming and consecutive faster glacial retreat modifies greatly the organization of novel plant communities at the very first years of primary succession.

The few quantitative studies focusing on the early arrival of plants after glacial retreat since the acceleration of warming, in the 70s, reveal a number of singularities that have not necessarily been observed in longer chronosequences. This includes the lack of significant trends in plant composition and in the value of intraspecific functional traits along with time since the glacial retreat (Sitzia et al., 2017: chronosequences 15–66 years), an overwhelming majority of wind-dispersed species and a deficit in facilitative interactions between plants (Zimmer et al., 2018: 0–38 years; Rosero et al., 2021: 0–61 years) and an accelerated pace of primary succession in comparison with measures made

before warming acceleration (Cannone et al., 2008: 1–25 years; Fickert et al., 2017: 1–35 years). Therefore, the accelerated pace of warming and glacial retreat is a game-changer for plants succession as it affords shorter time to perform upward migration (Klanderud and Totland, 2007). Shorter time may hinder site accessibility to seeds, especially those not dispersed by wind (Makoto and Wilson, 2019) and is expected to increase abiotic constraints on plant development as the soils do not have enough time to develop (Fastie, 1995; Sattin et al., 2009). Also, it may reduce the probability of positive interactions among species in stressful habitats where positive interactions have proven to be a central mechanism in the organization of plant communities (stress-gradient hypothesis; Bertness and Callaway, 1994; Anthelme et al., 2014). Consequently, to document the plant strategies designed to overcome the constraints limiting the first steps of primary succession is crucial to predict what types of alpine communities may develop in the next decades, and what species would face the highest extinction risks (Dullinger et al., 2012; Cauvy-Fraunié and Dangles, 2019).

Focusing on the study of plant primary succession after glacial retreat in the alpine, cold tropics is pertinent because they exhibit specific environmental features, with possibly unexpected impacts on successional processes. First, levels of biodiversity are globally higher due to specific climates, biogeography and isolation (Sklenár et al., 2014; Anthelme and Peyre, 2020) and plants exhibit narrower thermal niches breadth compared to higher latitudes (Cuesta et al., 2020). Second, warming is usually thought to be more pronounced (Pepin et al., 2015; Wang et al., 2016), which translates into faster glacial retreat (Vuille et al., 2018). Inverted precipitations gradients, a common feature in the alpine tropics at the origin of tropical alpine deserts (Leuschner, 2000) may increase stress for early colonists. A supplementary stress for the first colonizers may be the absence of a seasonal snowpack in the alpine tropics (Vuille et al., 2018) which exposes plants to low temperature at night all year round (e.g., Rada et al., 2001). Finally, during the dry period the absence of snowmelt generates higher soil aridity than in temperate-arctic environments (Vuille et al., 2018). For all these reasons, and despite the high stochasticity previously observed in plant distribution patterns along proglacial chronosequences (Marteinsdóttir et al., 2013; Buma et al., 2019), we expect to observe non-stochastic trends in plant distribution over the first years following glacial retreat in the alpine tropics.

Here we examined the first steps of plant primary succession, at a fine grain, in the foreland of a rapidly retreating glacier in the tropical high Andes, with the aim to reveal trends in plant distribution patterns. We used a 2–13 years chronosequence sequenced into four dated glacial retreat strips (Mount Antisana, Ecuador). In a 4,200 m² area, we measured plant traits at the community level and at the individual level and we compared them with data on the abiotic environment measured at each strip. We tested the following hypotheses: (1) the general trend of increasing plant diversity and productivity during early primary succession may be altered because of increasing stress a few years only after glacial retreat -for this purpose, we tested to what extent changes in soil macronutrients, temperature and soil humidity along the chronosequence modified plant
metrics, plant aggregation patterns and species-area relationships (Palmer, 2007; Storch et al., 2007)- and (2) facilitation among species are infrequent early after glacial retreat and successful colonizers are mainly dispersed by wind, two effects that may slow primary succession.

# METHODS

## Study Site

Our study site was located at 4870 m a.s.l. at the forefront of the glacier 15-alpha, on the western side of Antisana volcano (5,752 m), in the Eastern cordillera of the Ecuadorian Andes (0°29'06"S, 78°08'31"W; **Figure 1A**). The site lies in the "superpáramo" vegetation belt, where plant cover is discontinuous and plant communities are dominated by cushion-forming species [Azorella aretioides (Kunth) Willd. ex DC., Xenophyllum rigidum (Kunth) V.A.Funk], shrubs (Chuquiraga jussieui J.F.Gmel.) and tussock grasses [Calamagrostis intermedia (J.Presl) Steud.; Anthelme et al., 2012].

The glacier 15-alpha displayed a maximum extension during the Little-Ice Age, and reduced since with a severe acceleration of this reduction in the last 50 years (Jomelli et al., 2009). Inter-annual observations between 1997 and 2008 on glacial retreat using topographical measurements and differential GPS provided data for drawing a precise chronological map of recent glacial retreat (**Figure 1B**). An ORE meteorological station stands at only a few decametres of the sampling site. Data showed an annual rainfall amount of 907 mm in 2007 with a pronounced humid period between March and August (714 mm) and a monthly maximum of 276 mm in June. Between 2005 and 2007 annual mean temperature was 1.15°C and relative air humidity was 79.46%, without significant variations among months (Maisincho et al., 2009).

## Study Design

Our study site consisted of a 2–13 years proglacial chronosequence divided into four strips delimited by precisely-dated and marked lateral moraines of the glacier (**Figure 1C**). We selected a site on the lateral border of the glacier in order to reduce the effects of meltwater flows into our plots. The area considered in each strip was calculated from **Figure 1C** with the open access software ImageJ (http://imagej.nih.gov/ij/): strip 2–5 years (moraines 2008–2005; 984 m²); strip 5–7 years (2005–2003; 1,382 m²); strip 7–9 years (2003–2001; 902 m²); strip 9–13 years (2001–1997; 932 m²). Between moraines, the mineral granulometry at soil surface was equally dominated by small rocks up to 10 cm in diameter and sandy material. Rocks over 1 m in diameter were almost absent. Slope gradient followed the lines drawn by moraines and was about 10–15% (clinometer Suunto Tandem, Vantaa, Finland). As a result, both topography and granulometry between moraines were relatively similar within each strip, reducing the potential effect of local heterogeneity between strips on plant distribution patterns. Moraines were relatively small (between 0.1 and 2 m high) and distributed in all the glacial retreat strips (**Figure 1C; Supplementary Material 1**). Slopes on moraines were steeper, but with a similar granulometry as in the rest of the sampling (**Supplementary Material 1**). We

took them into account in our sampling as natural elements of the environment. Nevertheless, they were excluded from analyses related with species-area accumulation and from abiotic data analyses in order to focus on the effect of time since the glacial retreat on these variables.

## Abiotic Data

We define abiotic stress as the effects of external constraints on the limitation of plant dry biomass (Grime, 1977). Soil samples were collected 30 plots, each of 1 m²) for which we recorded the coordinates (**Figure 1**). Selection was random except they had similar slope and granulometry (coarse sand) at each of the four deglaciated strips. The aim was to study variation in soil moisture and soil macronutrient content along the chronosequence. Soil moisture was calculated as the difference between fresh and oven dried mass (48h, 70°C; Balance TU-OI FA-2104, Mhand, Fuzhou, China). In each strip, 10 of these soil samples were sent to the Instituto Nacional Autonomo de Investigaciones Agropecuarias (INIAP, Quito, Ecuador) for analyses of soil macronutrients (NH4, P, K, Ca, and Mg; modified Olsen methodology; Sanchez, 1976).

At each soil sample location (30 in each strip) we measured soil temperature at 5 cm depth, considering that temperature stress at this depth will be especially harmful for plants (Matthews, 1992). We made a series of 30 measurements between 9:00 and 10:00 AM during a cloudy day, using probe thermometer (Fluke 574CF, Everett, USA). In parallel, we recorded temperature each 15 min over 3 months (March-June 2010), using 10 TidbiT v2 Temp data loggers (two or three loggers at each strip) from which we extracted daily minimum temperature as a proxy of stress for plants. The location of the loggers was different from that of the previous temperature measurements but was on similar slope and granulometry. These two types of temperature measurements allowed us to take into account both the temporal variability and the spatial variability of temperatures among the four strips.

## Vegetation Data

### Vegetation Sampling

A 60 x 70 m² rectangular grid was superimposed on the study site across the four glacial retreat strips (**Figure 1C**). Unlike scattered plots, the grid method allowed us to spatialize collected information on plants and to document plant diversity and plant spatial associations patterns at various spatial scales, from 1 to 100 m². This grid was composed by 42 squares, each one divided into 100 plots of 1 m² (**Figure 1D**).

### Plant Metrics

In each plot of 1 m² (N = 4,200) all individuals of vascular plant species (N = 5,574) were recorded with XY coordinates (1–60 and 1–70 m, respectively) and identified in 2010. They were subject to various types of measurements (see below). Mosses and lichens were rare in the site and were not taken into account in the analysis. Two vascular species might be hybrids of two other Asteraceae species: *Senecio nivalis* (Kunth) Cuatrec. x *Lasiocephalus ovatus* Schltdl. (three individuals), and *Werneria nubigena* Kunth x X. *rigidum* (one individual); they
<figure><img src="image_3.png" /></figure>

**FIGURE 1** | Study area, study site and sampling design. (A) Volcano Antisana in the Ecuadorian Andes (extracted from Google map) and location of glacier 15-alpha within the Antisana ice-cap. (B) Map of the 15-alpha glacier between 1997 and 2008. Datum: WGMS 84; unit: UTM. (C) Sampling design: black areas represent the three lateral moraines of years 2001, 2003, and 2005, respectively; these areas were removed the species-area analyses because providing undesirable heterogeneity in the topography and granulometry. (D) Distribution of plots (1 m²).
were coded H1 and H2, respectively. Two species did not display flowering material and were identified at the family level only (Asteraceae: A1 and A2; two and eight individuals, respectively). Eleven juvenile individuals belonging to two unidentified morphospecies were classified as Z1 and Z2. Plants were segregated into mature individuals and seedlings; the latter category not being taken into account when measuring plant traits in order to avoid ontogenic effects. Seedlings were defined as individuals without sexual reproductive attributes (flowers, fruits, or seeds) and with a cover at soil surface not exceeding 1 cm².

Plant community changes were characterized based on plant diversity, species richness, plant cover and plant functional traits along the chronosequence. Plant cover provided information on general primary succession process. It was calculated inside each 1 m² plot by summing the surface of all individuals (N = 5,574), as the product between plant area (extrapolated by multiplying length and width of each individual; cm) and its relative density (visually estimated as the % of vegetal material within the area). In addition, we measured plant height (height of the highest leaf, cm), the proportion of flowering individuals, the number of ramets per individuals (proxy of vegetative reproduction) and plant necromass (visual estimation of dead tissues, %, following Caccianiga et al., 2006) to characterize the responses of the plant community to the increasing time after glacial retreat. These traits were measured for the 1,539 mature individuals. Each of these traits were studied at the community level by averaging the values of all individuals, and at the population level for the three most abundant species: *Agrostis foliata* Hook.f., *Cerastium floccosum* Benth. and *Senecio nivalis* (Kunth) Cuatrec.

## Spatial Associations Plant-Plant-Rock

We estimated the frequency of interspecific interactions among plants by observing the proportion of pairs of individuals spatially associated at each strip. We assumed associations to be a good proxy of positive interactions among plants because in alpine environments competition between plants is generally low (Cavieres et al., 2014). For this purpose, at each plot, we recorded each pair of spatially associated species. We assumed that spatial root interactions did not exceed aerial interactions among plants because the age of the target individuals did not exceed 13 years in our sampling, not giving sufficient time for plants to develop large root systems, as expected in mature alpine plants (Körner, 2003). Therefore, plants were considered spatially associated as soon as their aerial canopy overlapped. This method avoided the destruction of individuals to extract the roots. Nevertheless, we acknowledge that it may result in a slight underestimation of positive interactions. We also measured the proportion of individuals protected by an abiotic refuge, potentially providing protection from wind, moister conditions or temperature buffering to plants ("spatially associated with a nurse rock"; Scherrer and Körner, 2011; Moeslund et al., 2013). We considered plants as protected by a rock as soon as the rock exceeded the plant's height and the distance between plant and rock was not superior to the plant height.

## Data Analysis

### Chronosequence

Changes in plant parameters and abiotic parameters along the chronosequence (four strips) were first analyzed with simple linear regressions, plotting each variable against the distance to the glacier. The distance was the coordinate Y in our sampling as the sampling grid was parallel to the glacier border (Figure 1C). Thus, the distance to the glacier was used as a proxy of age of deglaciated terrains, which increased in the same direction as the distance. As an exception, data on Tmin were assigned to each strip without having XY coordinates. Therefore, we used the variable “strip” to provide a regression, by averaging the age of each strip. In the specific case of nurse plants and nurse rocks (modalities: protected or not protected) we used logistic regressions. Second, comparisons between means of each variable in each strip were provided with one-way ANOVA and post-hoc Tukey tests, or Kruskall-Wallis tests as soon as data did not follow normal distribution. This helped providing simple values along the chronosequence.

Correlations between two microclimatic variables (Tsoil 9-10 AM and soil moisture), three diversity variables (species richness, abundance of individuals and relative plant cover) and the time since glacier retreat were possible on the subsample of 120 plots used for the microclimatic measurements (see above). For this purpose, we conducted a principal component analysis on these data and we used “strip” as a supplementary variable to visualize the belonging of plots to each strip on the biplot graph. PCA was designed with ade4 1.7-16 (Thioulouse et al., 2018) and factoextra 1.0.7 (Kassambara and Mundt, 2017) in R.

### Spatial Analyses

We plotted all plant individuals as a marked point pattern at a one-meter resolution, with species identity as a mark [see Goreaud (2000), for a detailed presentation of the use of point pattern analysis in plant community ecology]. First, this spatial design allowed inferring species-area relationships (1) by calculating the frequency of patches without vegetation (hereafter termed “empty patches”) at various spatial scales, and (2) by providing species-area curves. Empty patches permitted exploring the frequency of suitable sites for plant establishment (hereafter referred to as safe sites; Harper et al., 1961). Along the chronosequence, we calculated the number of patches without vegetation at the four deglaciated strips, taking into account different sizes of patches (1, 4, 49, and 100 m²). For each strip, we calculated the proportion of empty patches (as a measure of environmental heterogeneity) for different patch sizes: 1, 4, 49, and 100 m², respectively. Moraine habitats, which were poorly colonized and under the influence of distinctive topographical parameters, were excluded from the analysis by making masks (see the R spatstat reference manual for details, http://cran.rproject.org/web/packages/spatstat/spatstat.pdf). For each patch size, we ran 100 simulations randomly located in each strip using a complete spatial randomness model, which is a completely random point set obtained from a homogeneous Poisson process. Using the same scale variation (1, 4, 49, and 100 m²), we drew species-area curves for each glacial retreat strip. We used the slope of each curve to generate a Preston model, which was the
most appropriate for such type IIa curves [type IIa: plot arrayed in a contiguous grid, see Scheiner (2003) for a description of the different species-area curve types]. The equation can be written as follows:

$$S = cA^z$$

with S the number of species, A the patch area, and c and z constants. The z parameter (slope of the curve on log-log axis) has been widely used by ecologists to make comparisons among different habitats and communities (Storch et al., 2007). These spatial analyses were performed using the SpatStat package of R (Baddeley and Turner, 2004).

# RESULTS

## Increasing Abiotic Stress Along the Chronosequence

Soil temperature between 9:00 and 10:00 AM, daily minimum soil temperature along a three-month period and soil moisture all decreased significantly with increasing distance to the glacier (linear regressions: R² = 0.41, R² =0.22, and R² = 0.19,

respectively; **Figure 2**). Soil temperature 9–10 AM changed significantly along the chronosequence, decreasing from 5.88°C ± 0.13 in the strip 2–5 years after glacial retreat to 4.54°C ± 0.11 in the strip 9–13 (post-hoc Tukey tests *P* < 0.05; **Supplementary Material 2**). Minimum soil temperatures followed the same trend, decreasing from 1.48°C ± 0.08 in the strip 2–5 years to 0.39°C ± 0.07 in the strip 9–13 (*P* < 0.05). Soil moisture was also significantly variable between strips, shifting from 12% ± 0.34 in strip 2–5 to 8.8% ± 0.37 in strip 9–13 (*P* < 0.05; **Supplementary Material 2**). For soil macronutrients, the concentration of Ca, K, and Mg decreased significantly along the chronosequence (linear regressions: R² = 0.43, R² =0.50, and R² = 0.08, respectively; **Figure 2**). As for microclimatic data, these nutrient values were significantly variables between strips, with higher values at the start of the chronosequence than at its end (*P* < 0.05).

## Plant Diversity and Composition

A total of 23 vascular plant species and 5,574 individuals were recorded in our sampling area (Table 1). Overall, three species dominated the community: Agrostis foliata (66% of observed

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Effects of the time since the glacial retreat on the microclimate and the relative content of three soil macronutrients (linear regressions; ***R² significant at *P* < 0.001). Time since the glacial retreat estimated through increasing distance from the glacier, except for Tsoil min (mean age of each strip). Different colors represent the belonging of plots to the four dated strips along the chronosequence. See also **Supplementary Material 2** for comparisons among strips and data on other soil nutrients.</figcaption></figure>
<table><caption>TABLE 1 | Species composition, dispersal mode, growth form, and number observed along each of the four strips of the proglacial chronosequence.</caption><thead><tr><th>Species</th><th>Family</th><th>Dispersal mode</th><th>Growth form</th><th>2-5 years</th><th>5-7 years</th><th>7-9 years</th><th>9-13 years</th><th>Total</th></tr></thead><tbody><tr><td>A1</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>1</td><td>0</td><td>1</td><td>0</td><td>2</td></tr><tr><td>A2</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>1</td><td>4</td><td>3</td><td>0</td><td>8</td></tr><tr><td>Agrostis foliata Hook.f.</td><td>Poaceae</td><td>A</td><td>Tussock grass</td><td>373</td><td>922</td><td>1146</td><td>1224</td><td>3665</td></tr><tr><td>Baccharis caespitosa (Ruiz and Pav.) Pers.</td><td>Asteraceae</td><td>A</td><td>Crawling shrub</td><td>3</td><td>5</td><td>6</td><td>0</td><td>14</td></tr><tr><td>Calamagrostis fibrovaginata Laegaard</td><td>Poaceae</td><td>A</td><td>Tussock grass</td><td>1</td><td>4</td><td>6</td><td>2</td><td>13</td></tr><tr><td>Calamagrostis intermedia (J.Presl) Steud.</td><td>Poaceae</td><td>A</td><td>Tussock grass</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><td>Cerastium floccosum Benth.</td><td>Caryophyllaceae</td><td>AZ</td><td>Rosette</td><td>70</td><td>474</td><td>412</td><td>513</td><td>1,469</td></tr><tr><td>Chuquiraga jussieui J.F.Gmel.</td><td>Asteraceae</td><td>A</td><td>Erect shrub</td><td>0</td><td>1</td><td>0</td><td>1</td><td>2</td></tr><tr><td>Ephedra americana Humb. and Bonpl. ex Willd.</td><td>Ephedraceae</td><td>Z</td><td>Crawling shrub</td><td>1</td><td>2</td><td>2</td><td>1</td><td>6</td></tr><tr><td>Erigeron ecuadoriensis Hieron.</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>0</td><td>2</td><td>0</td><td>0</td><td>2</td></tr><tr><td>H1 (L. ovatus X S. nivalis)</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>4</td><td>4</td><td>8</td><td>3</td><td>19</td></tr><tr><td>H2 (W. nubigena X X. rigidum)</td><td>Asteraceae</td><td>A</td><td>Cushion/Rosette</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><td>Hypochaeris sessiliflora agrgr.</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>0</td><td>2</td><td>1</td><td>0</td><td>3</td></tr><tr><td>Lasiocephalus ovatus Schltdl.</td><td>Asteraceae</td><td>A</td><td>Stem rosette</td><td>4</td><td>9</td><td>3</td><td>1</td><td>17</td></tr><tr><td>Belloa kunthiana (DC.) Anderb. and Freire</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>0</td><td>1</td><td>2</td><td>0</td><td>3</td></tr><tr><td>Plantago sericea Ruiz and Pav.</td><td>Plantaginaceae</td><td>AZ</td><td>Rosette</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><td>Senecio nivalis (Kunth) Cuatrec.</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>9</td><td>36</td><td>143</td><td>90</td><td>278</td></tr><tr><td>Taraxacum officinale L.</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>0</td><td>0</td><td>1</td><td>1</td><td>2</td></tr><tr><td>Werneria nubigena Kunth</td><td>Asteraceae</td><td>A</td><td>Rosette</td><td>13</td><td>17</td><td>8</td><td>3</td><td>41</td></tr><tr><td>Werneria pumila Kunth</td><td>Asteraceae</td><td>A</td><td>Cushion/Rosette</td><td>0</td><td>2</td><td>2</td><td>1</td><td>5</td></tr><tr><td>Xenophyllum rigidum (Kunth) V.A.Funk</td><td>Asteraceae</td><td>A</td><td>Cushion</td><td>0</td><td>4</td><td>6</td><td>1</td><td>11</td></tr><tr><td>Z1</td><td>-</td><td>?</td><td>?</td><td>1</td><td>2</td><td>0</td><td>0</td><td>3</td></tr><tr><td>Z2</td><td>-</td><td>?</td><td>?</td><td>3</td><td>2</td><td>3</td><td>0</td><td>8</td></tr><tr><td>Species richness (strip scale)</td><td></td><td></td><td></td><td>13</td><td>21</td><td>17</td><td>12</td><td>23</td></tr></tbody></table>

A1, A2, Z1, Z2: (partially) unidentified morphospecies. Dispersal modes, A: anemochory; Z: zoochory.

individuals), *Cerastium floccosum* (26%), and *Senecio nivalis* (5%). Rosette was the overwhelming dominant growth form before tussock grasses whereas cushion, shrub and stem rosette growth forms were represented by a limited number of species and individuals, especially in the older glacial retreat strip. Almost all species observed in the sampling were dispersed by wind, with the exception of *Ephedra americana* Humb. & Bonpl. ex Willd., dispersed only by animals. Relative plant cover at the community level in 1 m² plots ranged between 0.08 and 4.10% and was significantly variable among strips (Kruskall-Wallis: *P* < 0.001; **Supplementary Material 3**).

The two first axes of the PCA extracted from a subsample of 120 plots combining microclimatic and diversity data gathered 69.5% of the plot ordination (Figure 3). Axis 1 of the biplot confirmed that soil moisture and soil temperature were negatively correlated with time since glacier retreat (here estimated as the distance of plots from the glacier). The variability of plant diversity (species richness, abundance, cover) is expressed on axis 2, but not on axis 1 of the PCA, which means that these indices were not linearly correlated with time since the glacial retreat.

Species richness at the strip scale reached 13 species in strip 2–5, peaked in strip 5–7 with 21 species observed, and then decreased to 12 species in strip 9–13. Nine species that were present in the younger strips were absent in strip 9–13, among

which *Baccharis caespitosa* (Ruiz & Pav.) Pers., *C. intermedia*, *Hypochaeris sessiliflora* agrgr. and *Belloa kunthiana* (DC.) Anderb. & Freire. On the contrary, no species were exclusively observed in strip 9-13 (Table 1).

## Species-Area Relationship

Species richness increased with increasing area in each of the four strips of glacial retreat. However, the accumulation of species across area was reduced by half along the chronosequence [Preston coefficients of species-area curves (Z): 0.66 in strip 2–5 and 0.33 in strip 9–13; Spearman rank test *P* < 0.001; **Figure 4**]. This was the result of a higher frequency of large patches without vegetation in the older deglaciated strips (patches of 49 m²: 11% ± 0 in strip 2–5 and 28% ± 0 in strip 9–13; patches of 100 m²: 4% ± 0 in strip 2–5 and 23% ± 0 in strip 9–13; **Supplementary Material 4**).

## Individual Traits

A total of 1,539 mature individuals were observed in the sampling and were subject to individual trait measurements. Plant height at the community level decreased significantly along the chronosequence (linear regression: $R^2$ = 0.03; *P* < 0.001; **Figure 5**). It was significantly different between strip 2–5 and strip 9–13 (6.0 cm ± 0.3 and 4.5 cm ± 0.1, respectively; **Supplementary Material 5**). A similar significant trend was
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | PCA biplot combining diversity variables with microclimatic variables and the distance from the glacier on 120 plots. Each plot represented with a symbol assigned to one glacial retreat strip. Larger symbols represent the mean location of plots in each strip.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Species-area curves at each glacial retreat strip along the chronosequence. Z values (see Preston equation in the section Methods) displayed for each strip. Symbols represent mean values ± SD based on 100 random simulations.</figcaption></figure>

observed for *Agrostis foliata* and *Cerastium floccosum* ($R^2$ = 0.07 and $R^2$ = 0.05, respectively; *P* < 0.001) and, to a lesser extent, for *S. nivalis* ($R^2$ = 0.03; *P* < 0.05). The relative necromass in plants increased significantly along the chronosequence at

the community level (linear regression: $R^2$ = 0.09; *P* < 0.001). This trend was also observed in *A. foliata* and *C. floccosum* ($R^2$ = 0.09 and $R^2$ = 0.11, respectively; *P* < 0.001) but not in *S. nivalis* for which relative necromass remained stable at ~20% at each strip of the chronosequence ($R^2$ = 0.01, not significant). The proportion of flowering individuals increased at the community level and for *A. foliata* and *C. floccosum* (Logistic regressions: *P* < 0.001) whereas it remained stable -and lower than for other species- for *S. nivalis*. The number of ramets per individual did not change at the community level along the chronosequence (linear regression: $R^2$ = 0.00) and it remained low and stable across all strips (1.04 ± 0.01). This pattern was also observed in *A. foliata* and *C. floccosum* whereas *S. nivalis* developed more ramets than the rest of the community (1.55 ± 0.10; **Supplementary Material 5**).

## Spatial Associations

In the whole study site, almost half of the individuals were spatially associated with nurse rocks (722 out of 1,539), and this number significantly decreased along the chronosequence (Logistic regression: *z* = -3.32; *P* < 0.001; **Figure 6**). In contrast, only 68 individuals out of 1,539 were spatially associated with other plants. Still, the frequency of spatial associations among plants increased significantly along the chronosequence (Logistic regression: *z* = 6.17; *P* < 0.001), shifting from 0% of the individuals in strip 2-5 to 5.4% ± 1.0 in strip 9-13.
<figure><img src="image_7.png" /></figure>

**FIGURE 5** | Effect of time after glacial retreat on the plant height and the relative necromass of individuals along the chronosequence, at the community level (linear regressions; ***R² significant at *P* < 0.001). Time effect estimated through the distance from the glacier. Different colors represent the belonging of plots to the four dated strips along the chronosequence. See also **Supplementary Material 5** for comparisons among strips and data on complementary traits.

## DISCUSSION

Revealing the drivers of plant distribution and plant individual traits during the first steps primary succession after glacial retreat has been proven challenging because of high stochasticity, complex environmental gradients, small-scale topography, specific substrates slowing succession and, possibly, unsuitable sampling with high environmental heterogeneity among plots or insufficiently precise time delimitations between glacial retreat strips (Marteinsdóttir et al., 2013; D'Amico et al., 2017; Sitzia et al., 2017; Buma et al., 2019). The finesse of our grid-sampling (location of each of the 5,574 plant individuals at a one-meter resolution) is uncommon but not new as a number of other studies had a similar approach (e.g., Wiegand et al., 2012). However, combining this approach with a proglacial chronosequence distributed over such a short period of time is original. This allowed revealing increasing abiotic stress and a shift toward greater stress-tolerance by plants on a time scale of only 11 years. Although requiring further site replications, this result is a significant step toward a better understanding of the first steps in primary plant succession after glacial retreat.

## Plant Stress Increases During the First Steps of Succession

During primary succession, at the scale of the little ice age (250–350 years), it is commonly accepted that plants are gradually less constrained by the environment, as materialized by soil enrichment and stabilization, increase in both trophic and non-trophic interactions, increase in diversity and increase in local temperature (Cauvy-Fraunié and Dangles, 2019). However, during the first steps of glacial retreat, our results revealed that environmental constraints significantly increased between two and 13 years of glacial retreat; with lower soil temperature, soil moisture and soil macronutrient contents. This trend is in line with fragments of existing literature reporting (1) relative abundance of allochthonous nutrients just after glacial retreat as released by the glacier itself (Frenot et al., 1998; Caccianiga et al., 2006) and (2) a reduction in the allochthonous nutrient concentration in a few years only after glacial retreat, as being consumed by microbes (Sattin et al., 2009). At the same time enrichment in organic matter by plants themselves is likely negligible given the extremely low vegetation cover on the glacier foreland (Supplementary Material 3; Frenot et al., 1998). Decreasing humidity along the chronosequence may
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | Frequency of plant individuals developing under the shelter of abiotic refuges (nurse rocks) or heterospecific plants (nurse plants) at the community level, along the chronosequence. Mean values and Error bars: 95% CI are showed for each glacial retreat strip; different letters between strips indicate significant variations (post-hoc Tukey tests). The significance and the sign of the two trends were tested with logistic regressions: ***significant at *P* < 0.001.</figcaption></figure>

be the result of increasing distance with melting ice (Stöcklin and Bäumler, 1996; Frenot et al., 1998). In the specific case of the alpine tropics, the absence of seasonal snowpack, which provides regularly water at soil surface (Vuille et al., 2018), and decreasing precipitation at higher elevation (Leuschner, 2000) may exacerbate this trend as the melting glacier becomes the only source of water availability. Reduced humidity has itself a negative impact on temperature amplitude and soil nutrients thus explaining possibly the pronounced abiotic gradients that we observed along such a short chronosequence. The reduced local environmental variations (slope, granulometry) between the plots and the strips of our sampling contribute to increase the reliability of the observed changes. It is necessary, however, to confirm them in future studies by considering other sites at various latitudes to gain representativeness.

In parallel, increasing stress for plants was observable between two and 13 years after glacial retreat, in several ways. First, this has been observed by the fact that plant communities reduce their vegetative height, a classic response to the effects of alpine stress (Körner, 2003). Second, the plant richness at the scale of the strip was also lower at the end of the chronosequence, with especially the absence of nine species previously recorded in the succession, this trend being possibly attributable to increasing stress for plants since the negative correlation between species richness and abiotic stress is commonly pattern in the alpine zone (e.g., Kammer and Möhl, 2002). Third, the fact that the

most efficient clonal species in the sampling, *Senecio nivalis* (high abundance, high number of ramets and low investment in sexual reproduction) has not suffered an increase in necromass whereas the plant community does along the chronosequence goes with higher abiotic stress because clonal plants are generally more resistant to the alpine abiotic constraints (Stöcklin and Bäumler, 1996; Caccianiga et al., 2006). Fourth, the spatial distribution of plants at various spatial scales also evidenced diversity erosion along the chronosequence as both the reduction in the accumulation of species with increasing area (Preston curves) and in the proportion of safe sites (as evidenced by higher frequency of large empty patches) decreased with increasing time since the glacial retreat.

Focusing on the correlations between diversity, time since the glacial retreat and abiotic variables showed that diversity indices (abundance, species richness, plant cover) were not correlated linearly with time since they were variable along axis 2 of the PCA whereas time since the glacial retreat was variable along axis 1 (see **Figure 3**). This result contrasts with the positive correlation between plant diversity and time that is commonly observed during primary succession (Matthews, 1992; Zimmer et al., 2018; Cauvy-Fraunié and Dangles, 2019). An interpretation to be further tested is that the absence of linear relationship is the result of the antagonistic effects of time and the soil microclimate on diversity, which are opposite along axis 1 of the PCA. Such an antagonistic gradient may generate a unimodal relationship between time and diversity, which is not observable on the PCA biplot (based on linear correlations) but becomes visible when looking at the species richness at strip scale (see above). It could contribute to explain the lack of linear increase in diversity during early primary succession, as observed in some studies (e.g., Marteinsdóttir et al., 2013).

All in all, these results shed new light on recent studies observing a shift from relatively opportunistic to more stress-tolerant plant communities during the first steps of succession after glacial retreat. Indeed, three studies suggested that a shift from species with relative rapid growth rate (ruderals *sensu* Grime, 1977) to more stress-tolerant species occurred over a few years after glacial retreat (Caccianiga et al., 2006; Erschbamer and Caccianiga, 2016; Zimmer et al., 2018). By focusing on the first steps of colonization after glacial retreat (2–13 years), and by providing precise measurements on the abiotic environment (temperature, nutrients, soil moisture) our data provide that this trend can be explained with increasing abiotic stress. A decrease from 12 to 9% of the soil moisture between 2–5 years and 9–13 years, even statistically significant, might seem trivial to justify higher negative effects on plants. However, it represents a reduction by 25%, which, associated with decreases in temperature and nutrients, may be sufficient to exclude individuals of species that do not develop pronounced stress-tolerant traits. Such abiotic changes in longer chronosequences were not detected despite shifts from ruderal to stress-tolerant species (Zimmer et al., 2018, four chronosequences). This suggests (1) that increasing stress is observable only with a fine grain chronosequence, (2) that the consequences of this increasing stress after a few years can have an impact on plant traits on a longer term (see Zimmer et al., 2018) and (3) that
more evidence on this increasing stress is needed in other sites and regions.

## Deficit in Facilitation and Dispersal Filtering

The first steps of plant primary succession after glacial retreat are generally assumed to rely strongly on biotic interactions (Cauvy-Fraunié and Dangles, 2019). This is especially true for non-trophic, positive interactions among plants (facilitation: Crocker and Major, 1955; Erschbamer et al., 2008). Indeed, there is a strong relationship between abiotic stress and the intensity of facilitation (stress-gradient hypothesis; Bertness and Callaway, 1994), which is particularly well-documented in alpine regions (Cavieres et al., 2014), including in the alpine tropics (Anthelme and Dangles, 2012), and is sustained by habitat amelioration provided by the nurse/facilitator. The high frequency of plants associated with volcanic rocks in our sampling indicated that plants in recently deglaciated terrains do require habitat amelioration (temperature, soil moisture, macronutrients, protection from wind; Erschbamer et al., 2001; Scherrer and Körner, 2011). This result is in line with another work in the Ecuadorian superpáramo evidencing lava rocks as much better refuges for plants than sandy areas soon after glacial retreat (Suárez et al., 2015). Contrastingly, facilitation by nurse plants, which are able to provide better habitat ameliorations than rocks (Filazzola and Lortie, 2014), was particularly infrequent in our sampling (1.7% of heterospecific plants spatially associated). This contrast between strong facilitation by rocks and poor facilitation by plants reveals a severe deficit in nurse plants, a pattern that was also observed in the dry tropical high Andes (Zimmer et al., 2018). Our results support this hypothesis with the fact that giant cushion-forming species, probably the most efficient nurse life form in the (tropical) alpine regions (Reid et al., 2010; Cavieres et al., 2014) are absent in our sampling site although they are present in older, adjacent deglaciated terrains (Anthelme et al., 2012).

In combination with facilitation deficit, the difficulties encountered by plants not dispersed by the wind to colonize the recently deglaciated site may inhibit the first steps of primary succession (Makoto and Wilson, 2019). In a context of unprecedented velocity of temperature warming, the absence of keystone, cushion-forming species dispersed by animals during the first steps of primary succession, such as *A. aretioides* in our sampling or *D. muscoides* and *O. andina* in the southern tropical high Andes (Zimmer et al., 2018), could influence the composition and functioning of novel alpine ecosystems on the longer terms, just as local species pool does (Schumann et al., 2016). The disproportionate abundance of wind-dispersed species in our sampled area was observable through the overwhelming majority of species being wind-dispersed (see also Erschbamer et al., 2008; Zimmer et al., 2018; Franzén et al., 2019 in other recently deglaciated terrains) whereas “only” 39% of the species in the tropical alpine regions of the Andes are expected to be dispersed by wind (Tovar et al., 2020).

The deficit in facilitative interactions among plants combined with strong dispersal filtering in novel alpine plant communities will likely slow primary succession in the context of rapid climate change/rapid glacial retreat. As confronted to an increasing

abiotic stress during the few years after glacial retreat, plants may suffer even more from the absence of facilitation (stress-gradient hypothesis, Bertness and Callaway, 1994; Anthelme et al., 2014). This supports the already well-documented hypothesis that the collapse of positive interactions among organisms indebted to rapid climate change is an important driver of species vulnerability (Tilman, 1994; Gilman et al., 2010; Bellard et al., 2012). Nevertheless, the type of substrate, especially the presence of lava rocks, might also be an important driver of the plant colonization success, which deserves to be further tested (Suárez et al., 2015).

## CONCLUSION

Under the effects of an unprecedented velocity of warming, exploring the first steps of primary succession following glacial retreat is key to build scenarios on the organization and dynamic of novel alpine plant communities. The increasing abiotic stress for plants we found between two and 13 years of glacial retreat seems to alter the commonly observed trend that plant diversity increases along with time early after glacial retreat. It reveals that some of the first colonizing species, brought by wind, could be present in vain as being seemingly unable to withstand the harsher environment farther from the glacier and disappear. This interpretation is consistent with the view that, at a time when global changes are extremely rapid, it is greatly difficult for plants to find successful trade-offs between high dispersal capacities and specific adaptations to the environmental constraints of the new habitats, leading to impoverished biodiversity (Tilman, 1994), resulting in a dispersal lag combined with an establishment lag (Alexander et al., 2018). The deficit of nurse plants that we observed in our sampling makes the colonization process even more difficult for the plant communities. These results are obviously exploratory and require further site replications to test if these results are generalizable to the alpine tropics, where the absence of snowpack (Vuille et al., 2018) and reduced precipitation is expected to increase stress for plants in comparison with alpine environments outside the tropics (Leuschner, 2000).

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## AUTHOR CONTRIBUTIONS

FA, OD, SC-F, and BF designed the study. SC-F and FA collected the field data. SC-F, FA, and OD analyzed the results. BC produced the glacial retreat map of glacier 15 alpha. All authors contributed to the writing of the manuscript.

## FUNDING

The work was conducted within the frameworks of the International Joint Laboratories GREAT-ICE and BIOINCA, initiatives of the French Institute of Research for Development
(IRD) and universities and institutions in Ecuador (PUCE) and
Colombia (UniAndes).

ACKNOWLEDGMENTS

We thank Patricio Andino and Rodrigo Espinosa for their help in
collecting data in the field, Antoine Rabatel for discussions on the
dynamics of Andean glaciers, Gilles Le Moguedec for statistical

advices and two reviewers for their constructive comments on an
earlier version.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found
online at: https://www.frontiersin.org/articles/10.3389/fevo.
2021.584872/full#supplementary-material

REFERENCES

Alexander, J. M., Chalmandrier, L., Lenoir, J., Burgess, T. I., Essl, F., Haider, S., et al. (2018). Lags in the response of mountain plant communities to climate change. *Glob Change Biol.* 24, 563–579. doi: 10.1111/gcb.13976

Anthelme, F., Buendia, B., Mazoyer, C., and Dangles, O. (2012). Unexpected mechanisms sustain the stress gradient hypothesis in a tropical alpine environment. *J. Veg. Sci.* 23, 62–72. doi: 10.1111/j.1654-1103.2011.01333.x

Anthelme, F., Cavieres, L. A., and Dangles, O. (2014). Facilitation among plants in alpine environments in the face of climate change. *Front. Plant Sci.* 5:387. doi: 10.3389/fpls.2014.00387

Anthelme, F., and Dangles, O. (2012). Plant-plant interactions in tropical alpine environments. *Perspect. Plant Ecol. Evol. Syst.* 14, 363–372. doi: 10.1016/j.ppees.2012.05.002

Anthelme, F., and Peyre, G. (2020). "Biogeography of South American Highlands," in Reference Module in Earth Systems and Environmental Sciences, Encyclopedia of the World's Biomes, eds M. I. Goldstein and D. A. DellaSala (Amsterdam: Elsevier), 518–529. doi: 10.1016/B978-0-12-409548-9.11811-1

Baddeley, A. J., and Turner, R. (2004). "Spatstat: an R package for analyzing spatial point patterns," in R. Development Core Team R - A Language and Environment for Statistical Computing (Vienna: R Foundation for Statistical Computing). Available online at: http://www.r-project.org/. doi: 10.18637/jss.v012.i06 (accessed May, 2016).

Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W., and Courchamp, F. (2012). Impacts of climate change on the future of biodiversity. *Ecol. Lett.* 15, 365–377. doi: 10.1111/j.1461-0248.2011.01736.x

Bertness, M. D., and Callaway, R. (1994). Positive interactions in communities. *Trends Ecol. Evol.* 9, 191–193. doi: 10.1016/0169-5347(94)90088-4

Buma, B., Bisbing, S. M., Wiles, G., and Bidlack, A. L. (2019). 100 yr of primary succession highlights stochasticity and competition driving community establishment and stability. *Ecology* 100:e02885. doi: 10.1002/ecy.2885

Caccianiga, M., Luzzaro, A., Pierce, S., Ceriani, R. M., and Cerabolini, B. (2006). The functional basis of a primary succession resolved by CSR classification. *Oikos* 112, 10–20. doi: 10.1111/j.0030-1299.2006.14107.x

Cannone, N., Diolaiuti, G., Guglielmin, M., and Smiraglia, C. (2008). Accelerating climate change impacts on alpine glacier forefield ecosystems in the European Alps. *Ecol. Appl.* 18, 637–648. doi: 10.1890/07-1188.1

Cauvy-Fraunié, S., and Dangles, O. (2019). A global synthesis of biodiversity responses to glacier retreat. *Nat. Ecol. Evol.* 3, 1675–1685. doi: 10.1038/s41559-019-1042-8

Cavieres, L. A., Brooker, R. W., Butterfield, B. J., Cook, B. J., Kikvidze, Z., Lortie, C. J., et al. (2014). Facilitative plant interactions and climate simultaneously drive alpine plant diversity. *Ecol. Lett.* 17, 193–202. doi: 10.1111/ele.12217

Crocker, R. L., and Major, J. (1955). Soil development in relation to vegetation and surface age at Glacier Bay, Alaska. *J. Ecol.* 43, 427–448. doi: 10.2307/2257005

Cuesta, F., Tovar, C., Llambi, L. D., Gosling, W. D., Halloy, S., Carilla, J., et al. (2020). Thermal niche traits of high alpine plant species and communities across the tropical Andes and their vulnerability to global warming. *J. Biogeogr.* 47, 408–420. doi: 10.1111/jbi.13759

D'Amico, M. E., Freppaz, M., Zanini, E., and Bonifacio, E. (2017). Primary vegetation succession and the serpentine syndrome: the proglacial area of the Verra Grande glacier, North-Western Italian Alps. *Plant Soil* 415, 283–298. doi: 10.1007/s11104-016-3165-x

Dullinger, S., Willner, W., Plutzar, C., Englisch, T., Schratt-Ehrendorfer, L., Moser, D., et al. (2012). Post-glacial migration lag restricts range filling

of plants in the European Alps. *Global Ecol. Biogeogr.* 21, 829–840. doi: 10.1111/j.1466-8238.2011.00732.x

Erschbamer, B., and Caccianiga, M. S. (2016). “Glacier forelands: lessons of plant population and community development,” in *Progress in Botany*, eds F. Canovas, U. Lüttge, and R. Matyssek (Cham: Springer), 259–284. doi: 10.1007/124_2016_4

Erschbamer, B., Kneringer, E., and Niederfriniger Schlag, R. (2001). Seed rain, soil seed bank, seedling recruitment, and survival of seedlings on a glacier foreland in the Central Alps. *Flora* 196, 304–312. doi: 10.1016/S0367-2530(17)30059-2

Erschbamer, B., Niederfriniger Schlag, R., and Winkler, E. (2008). Colonization processes on a central alpine glacier foreland. *J. Veg. Sci.* 19, 855–862. doi: 10.3170/2008-8-18464

Fastie, C. L. (1995). Causes and ecosystem consequences of multiple pathways of primary succession at Glacier Bay, Alaska. *Ecology* 76, 1899–1916. doi: 10.2307/1940722

Fickert, T., Grüninger, F., and Damm, B. (2017). Klebensberg revisited: did primary succession of plants in glacier forelands a century ago differ from today? *Alp. Bot.* 127, 17–29. doi: 10.1007/s00035-016-0179-1

Filazzola, A., and Lortie, C. J. (2014). A systematic review and conceptual framework for the mechanistic pathways of nurse plants. *Global Ecol. Biogeog.* 23, 1335–1345. doi: 10.1111/geb.12202

Fischer, A., Fickert, T., Schwaizer, G., Patzelt, G., and Groß, G. (2019). Vegetation dynamics in Alpine glacier forelands tackled from space. *Sci. Rep.* 9, 1–13. doi: 10.1038/s41598-019-50273-2

Franzén, M., Dieker, P., Schrader, J., and Helm, A. (2019). Rapid plant colonization of the forelands of a vanishing glacier is strongly associated with species traits. *Arctic Antarc. Alp. Res.* 51, 366–378. doi: 10.1080/15230430.2019.1646574

Frenot, Y., Gloaguen, J., Cannavacchiolo, M., and Bellido, A. (1998). Primary succession on glacier forelands in the subantarctic Kerguelen Islands. *J. Veg. Sci.* 9, 75–84. doi: 10.2307/3237225

Gilman, S., Urban, M. C., Tewksbury, J., Gilchrist, G. W., and Holt, R. D. (2010). A framework for community interactions under climate change. *Trends Ecol. Evol.* 25, 325–331. doi: 10.1016/j.tree.2010.03.002

Goreaud, F. (2000). *Apports de l'analyse de la structure spatiale en forêt tempérée à l'étude et la modélisation des peuplements complexes* (PhD thesis), ENGREF, Montpellier, France, 524.

Grime, J. P. (1977). Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory. *Am. Nat.* 111, 1169–1194. doi: 10.1086/283244

Harper, J. L., Clatworthy, J. N., McNaughton, I. H., and Sagar, G. R. (1961). The evolution and ecology of closely related species living in the same area. *Evolution* 15, 209–227. doi: 10.1111/j.1558-5646.1961.tb03144.x

Heckmann, T., McColl, S. T., and Morche, D. (2015). Retreating ice: research in proglacial areas matters. *Earth Surf. Process. Landforms* 41, 271–276. doi: 10.1002/esp.3858

Jomelli, V., Favier, V., Rabatel, A., Brunstein, D., Hoffmann, G., and Francou, B. (2009). Fluctuations of glaciers in the tropical Andes over the last millennium and palaeoclimatic implications: a review. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* 281, 269–282. doi: 10.1016/j.palaeo.2008.10.033

Kammer, P. M., and Möhl, A. (2002). Factors controlling species richness in alpine plant communities: an assessment of the importance of stress and disturbance. *Arctic Antarc. Alp. Res.* 34, 398–407. doi: 10.1080/15230430.2002.1203510

Kassambara, A., and Mundt, F. (2017). *Package 'factoextra'. Extract and visualize the results of multivariate data analyses, CRAN Repository.
Klanderud, K., and Totland, Ø. (2007). The relative role of dispersal and local interactions for alpine plant community diversity under simulated climate warming. *Oikos* 116, 1279–1288. doi: 10.1111/j.0030-1299.2007.15906.x

Körner, C. (2003). *Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems*. Berlin: Heidelberg: Springer. doi: 10.1007/978-3-642-18970-8

Lenoir, J., and Svenning, J. C. (2015). Climate-related range shifts—a global multidimensional synthesis and new research directions. *Ecography* 38, 15–28. doi: 10.1111/ecog.00967

Leuschner, C. (2000). Are high elevations in tropical mountains arid environments for plants? *Ecology* 81, 1425–1436. doi: 10.1890/0012-9658(2000)081[1425:ANEITM]2.0.CO;2

Losapio, G., Jordán, F., Caccianiga, M., and Gobbi, M. (2015). Structure-dynamic relationship of plant-insect networks along a primary succession gradient on a glacier foreland. *Ecol. Model.* 314, 73–79. doi: 10.1016/j.ecolmodel.2015.07.014

Maisincho, L., Cáceres, B., Manciati, C., Loyo, C., Cuenca, E., Arias, M., et al. (2009). *Glaciares del Ecuador: Antisana y Carihuayrazo: balance de masa, topografía, pluviometría, meteorología e hidrología: informe del año 2007*.

Makoto, K., and Wilson, S. D. (2019). When and where does dispersal limitation matter in primary succession? *J. Ecol.* 107, 559–565. doi: 10.1111/1365-2745.12988

Marteinsdóttir, B., Thórhallsdóttir, T. E., and Svavarsdóttir, K. (2013). An experimental test of the relationship between small scale topography and seedling establishment in primary succession. *Plant Ecol.* 214, 1007–1015 doi: 10.1007/s11258-013-0226-6

Matthews, J. A. (1992). *The Ecology of Recently Deglaciated Terrain. A Geological Approach to Glacier Forelands and Primary Succession*. Cambridge: Cambridge University Press.

Moeslund, J. E., Arge, L., Bøcher, P. K., Dalgaard, T., and Svenning, J. C. (2013). Topography as a driver of local terrestrial vascular plant diversity patterns. *Nordic J. Bot.* 31, 129–144. doi: 10.1111/j.1756-1051.2013.00082.x

Palmer, M. W. (2007). "Species-area curves and the geometry of nature," in *Scaling Biodiversity*, eds D. Storch, P. A. Marquet, and J. H. Brown (Cambridge: Cambridge University Press), 15–31. doi: 10.1017/CBO9780511814938.004

Pepin, N., Bradley, R., Diaz, H., Baraer, M., Caceres, E., Forsythe, N., et al. (2015). Elevation-dependent warming in mountain regions of the world. *Nat. Clim. Change* 5, 424–430. doi: 10.1038/nclimate2563

Prach, K., and Walker, L. R. (2020). *Comparative Plant Succession Among Terrestrial Biomes of the World*. Cambridge: Cambridge University Press. doi: 10.1017/9781108561167

Rada, F., García-Núñez, C., Boero, C., Gallardo, M., Hilal, M., Gonzalez, J., et al. (2001). Low-temperature resistance in Polylepis tarapacana, a tree growing at the highest altitudes in the world. *Plant Cell Environ.* 24, 377–381. doi: 10.1046/j.1365-3040.2001.00685.x

Reid, A. M., Lamarque, L. J., and Lortie, C. J. (2010). A systematic review of the recent ecological literature on cushion plants: champions of plant facilitation. *Web Ecol.* 10, 44–49. doi: 10.5194/we-10-44-2010

Rosero, P., Crespo-Pérez, V., Espinosa, R., Andino, P., Barragán, A., Moret, P., et al. (2021). Multi-taxa colonisation along the foreland of a vanishing equatorial glacier. *Ecography*.

Sanchez, P. A. (1976). *Properties and Management of Soils in the Tropics*. Raleigh, NC: Soil Science Department, NCSU, Wiley Interscience, 618.

Sattin, S. R., Cleveland, C., Hood, E., Reed, S. C., King, A. J., Schmidt, S. K., et al. (2009). Functional shifts in unvegetated, perhumid, recently-deglaciated soil do not correlate with shifts in soil bacterial community composition. *J. Microb.* 47, 673–681. doi: 10.1007/s12275-009-0194-7

Scheiner, S. M. (2003). Six types of species-area curves. *Global Ecol. Biogeogr.* 12, 441–447. doi: 10.1046/j.1466-822X.2003.00061.x

Scherrer, D., and Körner, C. (2011). Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming. *J. Biogeogr.* 38, 406–416. doi: 10.1111/j.1365-2699.2010.02407.x

Schumann, K., Gewolf, S., and Tackenberg, O. (2016). Factors affecting primary succession of glacier foreland vegetation in the European Alps. *Alp. Bot.* 126, 105–117. doi: 10.1007/s00035-016-0166-6

Sitzia, T., Dainese, M., Krüsi, B. O., and McCollin, D. (2017). Landscape metrics as functional traits in plants: perspectives from a glacier foreland. *PeerJ* 5:e3552. doi: 10.7717/peerj.3552

Sklenár, P., Hedberg, I., and Cleef, A. M. (2014). Island biogeography of tropical alpine floras. *J. Biogeogr.* 41, 287–297. doi: 10.1111/jbi.12212

Stöcklin, J., and Bäumler, E. (1996). Seed rain, seedling establishment and clonal growth strategies on a glacier foreland. *J. Veg. Sci.* 9, 45–56. doi: 10.2307/3236415

Storch, D., Marquet, P. A., and Brown, J. H. (2007). *Scaling Biodiversity*. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511814938

Suárez, E., Orndahl, K., and Goodwin, K. (2015). Lava flows and moraines as corridors for early plant colonization of glacier forefronts on tropical volcanoes. *Biotropica* 47, 645–649. doi: 10.1111/btp.12260

Thioulouse, J., Dray, S., Dufour, A. B., Siberchicot, A., Jombart, T., and Pavoine, S. (2018). *Multivariate Analysis of Ecological Data With ade4*. New York, NY: Springer. doi: 10.1007/978-1-4939-8850-1

Tilman, D. (1994). "Community diversity and succession: the roles of competition, dispersal, and habitat modification," in *Biodiversity and Ecosystem Function*, eds E. D. Schulze, and H. A. Mooney (Berlin, Heidelberg: Springer), 327–344. doi: 10.1007/978-3-642-58001-7_15

Tovar, C., Melcher, I., Kusumoto, B., Cuesta, F., Cleef, A., Meneses, R. I., et al. (2020). Plant dispersal strategies of high tropical alpine communities across the Andes. *J. Ecol.* 108, 1910–1922. doi: 10.1111/1365-2745.13416

Vuille, M., Carey, M., Huggel, C., Buytaert, W., Rabatel, A., Jacobsen, D., et al. (2018). Rapid decline of snow and ice in the tropical Andes- Impacts, uncertainties and challenges ahead. *Earth Sci. Rev.* 176, 195–213. doi: 10.1016/j.earscirev.2017.09.019

Walker, L. R., and Del Moral, R. (2003). *Primary Succession and Ecosystem Rehabilitation*. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511615078

Wang, Q., Fan, X., and Wang, M. (2016). Evidence of high-elevation amplification vs. Arctic amplification. *Sci. Rep.* 6, 1–8. doi: 10.1038/srep19219

Wiegand, T., Huth, A., Getzin, S., Wang, X., Hao, Z., Gunatilleke, C. S., et al. (2012). Testing the independent species' arrangement assertion made by theories of stochastic geometry of biodiversity. *Proc. Royal Soc. B: Biol. Sci.* 279, 3312–3320. doi: 10.1098/rspb.2012.0376

Zemp, M., Huss, M., Thibert, E., Eckert, N., McNabb, R., Huber, J., et al. (2019). Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. *Nature* 568, 382–386. doi: 10.1038/s41586-019-1071-0

Zimmer, A., Meneses, R. I., Rabatel, A., Soruco, A., Dangles, O., and Anthelme, F. (2018). Time lag between glacial retreat and upward migration alters tropical alpine communities. *Perspect. Plant Ecol. Evol. Syst.* 30, 89–102. doi: 10.1016/j.ppees.2017.05.003

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2021 Anthelme, Cauvy-Fraunié, Francou, Cáceres and Dangles. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.