Esteban Suarez

Fabien Anthelme
email: fabien.anthelme@ird.fr

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

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

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 [START_REF] Matthews | The Ecology of Recently Deglaciated Terrain. A Geological Approach to Glacier Forelands and Primary Succession[END_REF][START_REF] Walker | Primary Succession and Ecosystem Rehabilitation[END_REF][START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF]. Although this type of primary succession depends on a high number of local, regional and latitudinal parameters [START_REF] D'amico | Primary vegetation succession and the serpentine syndrome: the proglacial area of the Verra Grande glacier, North-Western Italian Alps[END_REF][START_REF] Prach | Comparative Plant Succession Among Terrestrial Biomes of the World[END_REF] as well as stochastic effects [START_REF] Marteinsdóttir | An experimental test of the relationship between small scale topography and seedling establishment in primary succession[END_REF], 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 [START_REF] Losapio | Structure-dynamic relationship of plant-insect networks along a primary succession gradient on a glacier foreland[END_REF][START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF], increasing non-trophic interactions between plants (facilitation and competition [START_REF] Fastie | Causes and ecosystem consequences of multiple pathways of primary succession at Glacier Bay, Alaska[END_REF][START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF] and increasing plant cover and richness, as observed from field observations [START_REF] Matthews | The Ecology of Recently Deglaciated Terrain. A Geological Approach to Glacier Forelands and Primary Succession[END_REF][START_REF] Schumann | Factors affecting primary succession of glacier foreland vegetation in the European Alps[END_REF] and from satellite imagery analyses [START_REF] Fischer | Vegetation dynamics in Alpine glacier forelands tackled from space[END_REF]. 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 [START_REF] Tilman | Community diversity and succession: the roles of competition, dispersal, and habitat modification[END_REF].

Under the effects of accelerated warming for half a century, the pace of glacier retreat has dramatically increased across the globe [START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF]. 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 [START_REF] Heckmann | Retreating ice: research in proglacial areas matters[END_REF][START_REF] Lenoir | Climate-related range shifts-a global multidimensional synthesis and new research directions[END_REF]. At the same time, however, these new habitats are characterized by a very short temporal extent, which is expected to impact negatively their biodiversity [START_REF] Tilman | Community diversity and succession: the roles of competition, dispersal, and habitat modification[END_REF]. 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 [START_REF] Dullinger | Post-glacial migration lag restricts range filling of plants in the European Alps[END_REF]. 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 [START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF]. 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 [START_REF] Klanderud | The relative role of dispersal and local interactions for alpine plant community diversity under simulated climate warming[END_REF]. Shorter time may hinder site accessibility to seeds, especially those not dispersed by wind [START_REF] Makoto | When and where does dispersal limitation matter in primary succession?[END_REF] and is expected to increase abiotic constraints on plant development as the soils do not have enough time to develop [START_REF] Fastie | Causes and ecosystem consequences of multiple pathways of primary succession at Glacier Bay, Alaska[END_REF][START_REF] Sattin | Functional shifts in unvegetated, perhumid, recently-deglaciated soil do not correlate with shifts in soil bacterial community composition[END_REF]. 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; [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Anthelme | Facilitation among plants in alpine environments in the face of climate change[END_REF]. 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 [START_REF] Dullinger | Post-glacial migration lag restricts range filling of plants in the European Alps[END_REF][START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF].

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 [START_REF] Sklenár | Island biogeography of tropical alpine floras[END_REF][START_REF] Anthelme | Biogeography of South American Highlands[END_REF] and plants exhibit narrower thermal niches breadth compared to higher latitudes [START_REF] Cuesta | Thermal niche traits of high alpine plant species and communities across the tropical Andes and their vulnerability to global warming[END_REF]. Second, warming is usually thought to be more pronounced [START_REF] Pepin | Elevation-dependent warming in mountain regions of the world[END_REF][START_REF] Wang | Evidence of high-elevation amplification vs[END_REF], which translates into faster glacial retreat [START_REF] Vuille | Rapid decline of snow and ice in the tropical Andes-Impacts, uncertainties and challenges ahead[END_REF]. Inverted precipitations gradients, a common feature in the alpine tropics at the origin of tropical alpine deserts [START_REF] Leuschner | Are high elevations in tropical mountains arid environments for plants?[END_REF] 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 [START_REF] Vuille | Rapid decline of snow and ice in the tropical Andes-Impacts, uncertainties and challenges ahead[END_REF] which exposes plants to low temperature at night all year round (e.g., [START_REF] Rada | Low-temperature resistance in Polylepis tarapacana, a tree growing at the highest altitudes in the world[END_REF]. Finally, during the dry period the absence of snowmelt generates higher soil aridity than in temperate-arctic environments [START_REF] Vuille | Rapid decline of snow and ice in the tropical Andes-Impacts, uncertainties and challenges ahead[END_REF]. For all these reasons, and despite the high stochasticity previously observed in plant distribution patterns along proglacial chronosequences [START_REF] Marteinsdóttir | An experimental test of the relationship between small scale topography and seedling establishment in primary succession[END_REF][START_REF] Buma | 100 yr of primary succession highlights stochasticity and competition driving community establishment and stability[END_REF], 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 [START_REF] Palmer | Species-area curves and the geometry of nature[END_REF][START_REF] Storch | Scaling Biodiversity[END_REF]-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 cushionforming 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 [START_REF] Jomelli | Fluctuations of glaciers in the tropical Andes over the last millennium and palaeoclimatic implications: a review[END_REF]. 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 [START_REF] Maisincho | Glaciares del Ecuador: Antisana y Carihuayrazo: balance de masa, topografía, pluviometría, meteorología e hidrología: informe del[END_REF].

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 with the open access software ImageJ (http://imagej.nih.gov/ij/): strip 2-5 years (moraines 2008-2005; 984 m 2 ); strip 5-7 years (2005-2003; 1,382 m 2 ); strip 7-9 years ; 902 m 2 ); strip 9-13 years ; 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 [START_REF] Grime | Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory[END_REF]. 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; [START_REF] Sanchez | Properties and Management of Soils in the Tropics[END_REF].

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 [START_REF] Matthews | The Ecology of Recently Deglaciated Terrain. A Geological Approach to Glacier Forelands and Primary Succession[END_REF]. 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 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 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 [START_REF] Caccianiga | The functional basis of a primary succession resolved by CSR classification[END_REF] 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 [START_REF] Cavieres | Facilitative plant interactions and climate simultaneously drive alpine plant diversity[END_REF]. 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 [START_REF] Körner | Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems[END_REF]. 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"; [START_REF] Scherrer | Topographically controlled thermalhabitat differentiation buffers alpine plant diversity against climate warming[END_REF][START_REF] Moeslund | Topography as a driver of local terrestrial vascular plant diversity patterns[END_REF]. 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 [START_REF] Thioulouse | Multivariate Analysis of Ecological Data With ade4[END_REF] and factoextra 1.0.7 [START_REF] Kassambara | Package 'factoextra'. Extract and visualize the results of multivariate data analyses[END_REF] 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 [START_REF] Goreaud | Apports de l'analyse de la structure spatiale en forêt tempérée à l'étude et la modélisation des peuplements complexes[END_REF], 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; [START_REF] Harper | The evolution and ecology of closely related species living in the same area[END_REF]. Along the chronosequence, we calculated the number of patches without vegetation at the four deglaciated strips, taking into account different sizes of patches [START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF]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 [START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF]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 [START_REF] Scheiner | Six types of species-area curves[END_REF] 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 [START_REF] Storch | Scaling Biodiversity[END_REF]. These spatial analyses were performed using the SpatStat package of R [START_REF] Baddeley | Spatstat: an R package for analyzing spatial point pattens[END_REF].

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 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. 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 aggr. 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 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 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 [START_REF] Marteinsdóttir | An experimental test of the relationship between small scale topography and seedling establishment in primary succession[END_REF][START_REF] D'amico | Primary vegetation succession and the serpentine syndrome: the proglacial area of the Verra Grande glacier, North-Western Italian Alps[END_REF][START_REF] Sitzia | Landscape metrics as functional traits in plants: perspectives from a glacier foreland[END_REF][START_REF] Buma | 100 yr of primary succession highlights stochasticity and competition driving community establishment and stability[END_REF]. 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., [START_REF] Wiegand | Testing the independent species' arrangement assertion made by theories of stochastic geometry of biodiversity[END_REF]. 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 [START_REF] Frenot | Primary succession on glacier forelands in the subantarctic Kerguelen Islands[END_REF][START_REF] Caccianiga | The functional basis of a primary succession resolved by CSR classification[END_REF] and ( 2) a reduction in the allochthonous nutrient concentration in a few years only after glacial retreat, as being consumed by microbes [START_REF] Sattin | Functional shifts in unvegetated, perhumid, recently-deglaciated soil do not correlate with shifts in soil bacterial community composition[END_REF]. 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; [START_REF] Frenot | Primary succession on glacier forelands in the subantarctic Kerguelen Islands[END_REF]. Decreasing humidity along the chronosequence may be the result of increasing distance with melting ice [START_REF] Stöcklin | Seed rain, seedling establishment and clonal growth strategies on a glacier foreland[END_REF][START_REF] Frenot | Primary succession on glacier forelands in the subantarctic Kerguelen Islands[END_REF]. In the specific case of the alpine tropics, the absence of seasonal snowpack, which provides regularly water at soil surface [START_REF] Vuille | Rapid decline of snow and ice in the tropical Andes-Impacts, uncertainties and challenges ahead[END_REF], and decreasing precipitation at higher elevation [START_REF] Leuschner | Are high elevations in tropical mountains arid environments for plants?[END_REF] 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 [START_REF] Körner | Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems[END_REF]. 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., [START_REF] Kammer | Factors controlling species richness in alpine plant communities: an assessment of the importance of stress and disturbance[END_REF]. 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 [START_REF] Stöcklin | Seed rain, seedling establishment and clonal growth strategies on a glacier foreland[END_REF][START_REF] Caccianiga | The functional basis of a primary succession resolved by CSR classification[END_REF]. 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 [START_REF] Matthews | The Ecology of Recently Deglaciated Terrain. A Geological Approach to Glacier Forelands and Primary Succession[END_REF][START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF][START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF]. 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., [START_REF] Marteinsdóttir | An experimental test of the relationship between small scale topography and seedling establishment in primary succession[END_REF].

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 [START_REF] Grime | Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory[END_REF] to more stress-tolerant species occurred over a few years after glacial retreat [START_REF] Caccianiga | The functional basis of a primary succession resolved by CSR classification[END_REF][START_REF] Erschbamer | Glacier forelands: lessons of plant population and community development[END_REF][START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF]. 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 [START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF] 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 [START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF]. This is especially true for non-trophic, positive interactions among plants (facilitation: [START_REF] Crocker | Soil development in relation to vegetation and surface age at Glacier Bay, Alaska[END_REF][START_REF] Erschbamer | Colonization processes on a central alpine glacier foreland[END_REF]. Indeed, there is a strong relationship between abiotic stress and the intensity of facilitation (stress-gradient hypothesis; [START_REF] Bertness | Positive interactions in communities[END_REF], which is particularly well-documented in alpine regions [START_REF] Cavieres | Facilitative plant interactions and climate simultaneously drive alpine plant diversity[END_REF], 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; [START_REF] Erschbamer | Seed rain, soil seed bank, seedling recruitment, and survival of seedlings on a glacier foreland in the Central Alps[END_REF][START_REF] Scherrer | Topographically controlled thermalhabitat differentiation buffers alpine plant diversity against climate warming[END_REF]. 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 [START_REF] Suárez | Lava flows and moraines as corridors for early plant colonization of glacier forefronts on tropical volcanoes[END_REF]. Contrastingly, facilitation by nurse plants, which are able to provide better habitat ameliorations than rocks [START_REF] Filazzola | A systematic review and conceptual framework for the mechanistic pathways of nurse plants[END_REF], 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 [START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF]. 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 [START_REF] Reid | A systematic review of the recent ecological literature on cushion plants: champions of plant facilitation[END_REF][START_REF] Cavieres | Facilitative plant interactions and climate simultaneously drive alpine plant diversity[END_REF] 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 [START_REF] Makoto | When and where does dispersal limitation matter in primary succession?[END_REF]. 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 [START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF], could influence the composition and functioning of novel alpine ecosystems on the longer terms, just as local species pool does [START_REF] Schumann | Factors affecting primary succession of glacier foreland vegetation in the European Alps[END_REF]. The disproportionate abundance of wind-dispersed species in our sampled area was observable through the overwhelming majority of species being wind-dispersed (see also [START_REF] Erschbamer | Colonization processes on a central alpine glacier foreland[END_REF][START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF][START_REF] Franzén | Rapid plant colonization of the forelands of a vanishing glacier is strongly associated with species traits[END_REF] 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 [START_REF] Tovar | Plant dispersal strategies of high tropical alpine communities across the Andes[END_REF].

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, [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Anthelme | Facilitation among plants in alpine environments in the face of climate change[END_REF]. 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 [START_REF] Tilman | Community diversity and succession: the roles of competition, dispersal, and habitat modification[END_REF][START_REF] Gilman | A framework for community interactions under climate change[END_REF][START_REF] Bellard | Impacts of climate change on the future of biodiversity[END_REF]. 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 [START_REF] Suárez | Lava flows and moraines as corridors for early plant colonization of glacier forefronts on tropical volcanoes[END_REF].

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 [START_REF] Tilman | Community diversity and succession: the roles of competition, dispersal, and habitat modification[END_REF], resulting in a dispersal lag combined with an establishment lag [START_REF] Alexander | Lags in the response of mountain plant communities to climate change[END_REF]. 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 [START_REF] Vuille | Rapid decline of snow and ice in the tropical Andes-Impacts, uncertainties and challenges ahead[END_REF] and reduced precipitation is expected to increase stress for plants in comparison with alpine environments outside the tropics [START_REF] Leuschner | Are high elevations in tropical mountains arid environments for plants?[END_REF].

(IRD) and universities and institutions in Ecuador (PUCE) and Colombia (UniAndes).

FIGURE 1 |

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

FIGURE 3 |

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

FIGURE 4 |

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

FIGURE 5 |

 5 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 2 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.

FIGURE 6 |

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

TABLE 1 |

 1 Species composition, dispersal mode, growth form, and number observed along each of the four strips of the proglacial chronosequence.

	Species	Family	Dispersal mode Growth form	2-5 years 5-7 years 7-9 years 9-13 years Total
	A1	Asteraceae	A	Rosette	1	0	1	0	2
	A2	Asteraceae	A	Rosette	1	4	3	0	8
	Agrostis foliata Hook.f.	Poaceae	A	Tussock grass	373	922	1146	1224	3665
	Baccharis caespitosa (Ruiz and Pav.) Pers.	Asteraceae	A	Crawling shrub	3	5	6	0	14
	Calamagrostis fibrovaginata Laegaard	Poaceae	A	Tussock grass	1	4	6	2	13
	Calamagrostis intermedia (J.Presl) Steud.	Poaceae	A	Tussock grass	0	1	0	0	1
	Cerastium floccosum Benth.	Caryophyllaceae AZ	Rosette	70	474	412	513	1,469
	Chuquiraga jussieui J.F.Gmel.	Asteraceae	A	Erect shrub	0	1	0	1	2
	Ephedra americana Humb. and Bonpl. ex Willd. Ephedraceae	Z	Crawling shrub	1	2	2	1	6
	Erigeron ecuadoriensis Hieron.	Asteraceae	A	Rosette	0	2	0	0	2
	H1 (L. ovatus X S. nivalis)	A s t e r a c e a e	A	R o s e t t e	4	4	8	3	1 9
	H2 (W. nubigena X X. rigidum)	Asteraceae	A	Cushion/Rosette	0	1	0	0	1
	Hypochaeris sessiliflora aggr.	Asteraceae	A	Rosette	0	2	1	0	3
	Lasiocephalus ovatus Schltdl.	Asteraceae	A	Stem rosette	4	9	3	1	17
	Belloa kunthiana (DC.) Anderb. and Freire	Asteraceae	A	Rosette	0	1	2	0	3
	Plantago sericea Ruiz and Pav.	Plantaginaceae	AZ	Rosette	0	1	0	0	1
	Senecio nivalis (Kunth) Cuatrec.	Asteraceae	A	Rosette	9	36	143	90	278
	Taraxacum officinale L.	Asteraceae	A	Rosette	0	0	1	1	2
	Werneria nubigena Kunth	Asteraceae	A	Rosette	13	17	8	3	41
	Werneria pumila Kunth	Asteraceae	A	Cushion/Rosette	0	2	2	1	5
	Xenophyllum rigidum (Kunth) V.A.Funk	Asteraceae	A	Cushion	0	4	6	1	11
	Z1	-	?	?	1	2	0	0	3
	Z2	-	?	?	3	2	3	0	8
	Species richness (strip scale)				13	21	17	12	23
	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 2 plots ranged between 0.08 and 4.10%
	and was significantly variable among strips (Kruskall-Wallis:
	P < 0.001; Supplementary Material 3).

Frontiers in Ecology and Evolution | www.frontiersin.org

February 2021 | Volume 9 | Article 584872

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.

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

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.

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

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.