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

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 m2 divided into plots of 1 m2. 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.

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 2 divided into plots of 1 m 2 . 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.

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 . 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 2 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.
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 preciselydated 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   ; 932 m 2 ). 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 2 ) 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 Sampling
A 60 x 70 m 2 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 2 . This grid was composed by 42 squares, each one divided into 100 plots of 1 m 2 ( Figure 1D).

Plant Metrics
In each plot of 1 m 2 (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 Frontiers in Ecology and Evolution | www.frontiersin.org 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 2 .
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 2 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 . 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 T min 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 2 ). 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 2 , 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 2 ), 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).

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 2 = 0.41, R 2 =0.22, and R 2 = 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 2 = 0.43, R 2 =0.50, and R 2 = 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 2 | Effects of the time since the glacial retreat on the microclimate and the relative content of three soil macronutrients (linear regressions; ***R 2 significant at P < 0.001). Time since the glacial retreat estimated through increasing distance from the glacier, except for T soil 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. 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 2 plots ranged between 0.08 and 4.10% and was significantly variable among strips (Kruskall-Wallis: 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-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 2 : 11% ± 0 in strip 2-5 and 28% ± 0 in strip 9-13; patches of 100 m 2 : 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 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. 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.

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 gridsampling (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 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 stresstolerant 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 stresstolerant 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 , including in the alpine tropics , 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 .
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 .
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.