# The importance of species addition ‘versus’ replacement varies over succession in plant communities after glacier retreat

Isabel Cantera, Alexis Carteron, Alessia Guerrieri, Silvio Marta, Aurélie Bonin, Roberto Ambrosini, Fabien Anthelme, Roberto Sergio Azzoni, Peter Almond, Pablo Alviz Gazitúa, et al.

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

## To cite this version:

Isabel Cantera, Alexis Carteron, Alessia Guerrieri, Silvio Marta, Aurélie Bonin, et al. The importance of species addition ‘versus’ replacement varies over succession in plant communities after glacier retreat. Nature Plants, 2024, 10, pp.256-267. 10.1038/s41477-023-01609-4. hal-04411640

HAL Id: hal-04411640

https://hal.inrae.fr/hal-04411640v1

Submitted on 8 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The importance of species addition ‘versus’ replacement varies over succession in plant communities after glacier retreat

Received: 16 January 2023

Accepted: 6 December 2023

Published online: 17 January 2024

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

Check for updates

A list of authors and their affiliations appears at the end of the paper

The mechanisms underlying plant succession remain highly debated. Due to the local scope of most studies, we lack a global quantification of the relative importance of species addition ‘versus’ replacement. We assessed the role of these processes in the variation (β-diversity) of plant communities colonizing the forelands of 46 retreating glaciers worldwide, using both environmental DNA and traditional surveys. Our findings indicate that addition and replacement concur in determining community changes in deglaciated sites, but their relative importance varied over time. Taxa addition dominated immediately after glacier retreat, as expected in harsh environments, while replacement became more important for late-successional communities. These changes were aligned with total β-diversity changes, which were more pronounced between early-successional communities than between late-successional communities (>50 yr since glacier retreat). Despite the complexity of community assembly during plant succession, the observed global pattern suggests a generalized shift from the dominance of facilitation and/or stochastic processes in early-successional communities to a predominance of competition later on.

Ecological successions—how communities change or replace one another over time—have been a cornerstone of ecology since its inception¹,². Primary successions, that is, the development of ecosystems where a severe disturbance opens up large areas lacking most lifeforms³, start when a given species or a set of species colonize a newly exposed surface, which would then be further colonized by other species to reach complex communities³. Plant communities have been a major focus of primary succession studies for over a century since Clements’s⁴ work. Still, despite decades of work, the mechanisms that drive plant primary succession remain not fully understood²,⁵. In the pioneering Clementsian deterministic view, succession occurs in a progressive, directional and homogeneous manner, with a relatively stable and predictable community structure (the climax stage) reached after some time without punctuated changes⁴,⁵. This deterministic predictability was rapidly questioned and debated⁶. Over the past century, some studies showed homogeneous responses among plant successions⁷,⁸, while others showed that successions are not always convergent⁹,¹⁰, barely reach an equilibrium³ and are largely determined by stochastic processes¹¹.

While there is no clear concession about their determinism or stochasticity, successions are generally characterized by an increase in the number of plant species over time, accompanied by changes in community composition²,³ (but see ref. 12). These changes can be the result of two non-exclusive mechanisms: changes in species richness (due to species addition) and species replacement. Under a mechanism of species addition during succession, a community at time *t* is a subset of the species assemblage at *t + 1* because of the persistence of early colonizers (that is, pioneer species) and the addition of new species. In contrast, replacement involves the substitution over time of early colonizers by other species¹,¹³. The relative importance of these two key mechanisms in plant succession is difficult to quantify as it may vary over the succession period, but also as a function of the local topographic and environmental characteristics²,⁸,¹²,¹⁴. Thus, addressing the importance of species addition versus replacement requires the assessment of plant successions in multiple environmental settings.

Ongoing climate change is dramatically accelerating the retreat of glaciers worldwide¹⁵,¹⁶, exposing new terrains to the development of

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

e-mail: isa_cantera@hotmail.com
plant successions¹³. As such, glacier forelands provide opportunities to study plant primary successions. Recently deglaciated terrains are typically isolated and characterized by harsh conditions, including a cold climate, unstable substrate and limited nutrient availability¹⁷. Thus, these emerging ecosystems typically show sparse vegetation¹⁸ dominated by few species that exhibit a high degree of specialization to live in these conditions¹⁹. In this context, successions are generally dominated by a gradual addition of species during ecosystem development²⁰,²¹, while replacement is often hypothesized to be weak or absent²,¹²,¹⁴. It is, nonetheless, possible that the relative importance of species addition versus replacement changes over time, and the short time window covered by most studies²⁰,²¹ could have overemphasized the importance of addition.

In principle, species accumulation could result from the joint effect of dispersal-related processes (that is, species require time to disperse in a new area) and facilitative interactions⁷,¹²,²²⁻²⁴, which have a key role under severe environmental conditions²⁵⁻²⁷. Plant species that initially colonize recently deglaciated terrains can modify the environment through the accumulation of nutrients and organic matter, and can create new micro-environmental conditions²⁵,²⁶,²⁸⁻³⁰ that are more suitable to the establishment and development of subsequent colonizers¹. Nevertheless, competitive exclusions and facilitation were found to jointly affect biodiversity in alpine plant communities³¹. Their relative contributions vary along environmental gradients¹² and are expected to change over time along the glacier foreland. For instance, the increase in nutrient availability and the reduction in physical constraints over succession may progressively reduce limiting conditions, stimulating growth and reproduction³,¹⁹,³². In the long term, this can intensify competition for space and light²⁴,³³,³⁴. Moreover, nutrient availability generally increases during the first stages of succession, but can reach saturation in late stages²⁵,³⁴,³⁵. The resulting nutrient limitation can favour late-successional species that are efficient nutrient users²⁵,³⁴,³⁶. Therefore, species addition is expected to play a dominant role in plant succession soon after glacier retreat, with its importance decreasing as replacement becomes more important over time. However, a proper quantification of these expectations has never been carried out so far, as it requires a global dataset that spans large temporal and spatial scales.

In this study, we addressed this challenge by measuring and decomposing the β-diversity (that is, compositional variation between assemblages³⁷) of plant communities colonizing 46 glacier forelands distributed worldwide (Fig. 1a). Glacier forelands include chronosequences of progressively older terrains at growing distances from the glacier forefront³⁸. We used this space-for-time substitution to encompass communities covering a wide range of time since glacier retreat (sites deglaciated from 1 to >400 yr); each site represents the past position of the glacier at known dates (Fig. 1b). Along these chronosequences, we collected environmental DNA (eDNA) from soil to reconstruct 266 communities of vascular plants, that is, assemblages of taxa inhabiting a deglaciated site with a specific age since glacier retreat. eDNA metabarcoding allows the production of inventories consistent with traditional inventories of aboveground plants³⁹⁻⁴¹ and enables the rapid assessment of communities over broad geographic scales and from remote areas, yielding data that would have been challenging to assemble with traditional methods¹²,⁴³. Although soil eDNA is increasingly used to analyse plant communities, some uncertainty persists on the correspondence between patterns obtained with this approach and those obtained with traditional sampling³⁹,⁴⁴,⁴⁵. Thus, we compared eDNA-based patterns with the ones obtained through morphological identification of species ('traditional sampling') for a subset of forelands (Fig. 1a).

To understand how compositional variation changes over plant succession after glacier retreat, we first quantified the total β-diversity (β-total) between communities within each foreland and assessed whether it decreases over the succession²². Second, we tested the

hypothesis that taxa addition prevails over taxa replacement soon after glacier retreat, but its importance decreases along the succession. To this end, β-total was decomposed into its β-richness (that is, richness differences between communities due to species gain) and β-replacement components (that is, the substitution of one taxon by another, without affecting species richness), and the temporal patterns of these components were assessed. Finally, we assessed whether the detected patterns are progressive or show breakpoints over successions²².

## Diversity of plant communities in glacier forelands

After data filtering and removal of likely contaminants, 519 molecular operational taxonomic units (MOTUs) of vascular plants were detected with eDNA metabarcoding across the dated sites of the 46 glacier forelands (Supplementary Tables 1 and 2). eDNA metabarcoding detected 0–60 MOTUs per site (mean ± s.e. = 11 ± 0.6; Supplementary Fig. 1a) and 5–150 MOTUs per foreland (mean = 36 ± 4.4). No plant MOTUs were detected at 17 sites, all of which were <32 yr old (Supplementary Table 1). Using traditional floristic surveys, 365 vascular plant taxa were detected across 57 dated sites from 13 glacier forelands (Supplementary Tables 1 and 3), with 0–133 (mean = 24 ± 3) taxa per site (Supplementary Fig. 1a) and 15–162 taxa (mean = 48 ± 10) per foreland. Traditional sampling did not show any vascular plant species only in one site, which was deglaciated since just 2 yr (Carihuairazo glacier; Supplementary Table 1).

To assess the variation in community composition over time, as well as the relative contribution of taxa addition versus replacement to successions, we quantified and partitioned the β-total between pairs of communities. We compared communities within each foreland (excluding comparisons between communities belonging to different forelands), also including comparisons that do not correspond to directly subsequent age classes (Fig. 1b). This resulted in 771 and 102 comparisons of communities for the eDNA and the traditional data, respectively. β-total was partitioned into taxa replacement (β-replacement) and richness differences (β-richness) following ref. 46 (see Methods). Overall, the contributions of β-richness and β-replacement to β-total were comparable (Fig. 2; eDNA data: mean contribution of β-richness and β-replacement was 54% and 46%, respectively; traditional data: mean contribution of β-richness and β-replacement was 49.8% and 50.2%, respectively). β-richness values tended to be significantly higher than β-replacement values for eDNA data (randomization test for paired samples: *P* < 0.001) but not for traditional data (*P* = 0.773).

## Variation in β-diversity components over succession

Bayesian generalized linear mixed models (GLMMs) were used to assess how the β-diversity components varied over succession, considering two independent variables: 'age differences' between communities and their 'mean age' (Fig. 1b). The age differences indicate how different the communities are from each other in terms of age since deglaciation. Low values represent comparisons between communities in similar successional stages (for example, early versus early or late versus late), while high values represent comparisons between communities at very different successional stages (for example, early versus late). The mean age is the averaged time since glacier retreat between the compared communities; low values represent comparisons between young communities, while high values represent comparisons between late-successional communities.

For both eDNA and traditional data, β-total increased with age differences (the 95% credible intervals (CIs) of this effect were consistently positive; Fig. 3a–c and Table 1). Furthermore, the dissimilarity between old communities was generally smaller than the dissimilarity between young communities with similar age differences (the mean age of compared communities showed effects with 95% CIs consistently
<figure><img src="image_5.png" /></figure>

**Fig. 1 | Sampling design.** a, Global distribution of the 46 glacier forelands where plant communities were sampled in dated sites along chronosequences with eDNA from the soil (all circles). In 13 of these glacier forelands, we also gathered traditional plant inventories (purple and yellow circles). The background blue grid represents the number of glaciers for each 1 × 1° cell (www.glims.org) and ranges from 1 (pale blue) to 3,500 glaciers (darkest blue) (see colour bar). b, Sampling scheme used for the eDNA approach. For each of the 46 forelands, we identified 3–17 sites along the chronosequences. Each site represents the past position of the glacier at known dates (each colour indicates a site and its corresponding age class). For each site, we established -5 plots (diamonds). Within each plot, we collected 5 soil subsamples within 1 m² (subsamples

distribution is shown by the pink inset); subsamples were pooled into one composite sample per plot. The taxa detected in the different plots of the same site were combined to inventor plant communities in each of the 266 sites. For each pair of sites from the same glacier foreland, we calculated the total β-diversity and its β-replacement and β-richness components, as well as the mean age and the age differences between the compared communities. Comparisons between sites from different forelands were not assessed. We show examples of the calculation of age differences and mean age variables for a subset of sites. We considered all the pairwise comparisons, including comparisons that do not correspond to directly subsequent age classes in the chronosequence.

negative for eDNA data and slightly overlapping zero for traditional data; Fig. 3b-d and Table 1).

β-richness and β-replacement showed distinct responses to mean age and age differences. β-richness increased with age differences between compared communities (consistently positive CIs; Fig. 3a-c and Table 1), while it decreased with the mean age of compared communities (Fig. 3b-d and Table 1), with similar patterns between eDNA and traditional data. The relationship between β-replacement and age differences was negative but generally weak, with CIs overlapping zero for eDNA and marginally so for traditional data (Fig. 3a-c and Table 1). Finally, β-replacement tended to increase with the mean age of communities (CIs consistently positive for traditional data and slightly overlapping zero for eDNA; Fig. 3b-d and Table 1). Despite some differences across forelands, the relationships remained consistent when we excluded the forelands with a small number of sites (<8 sites per foreland; Supplementary Figs. 2 and 3). The pattern also remained consistent when we repeated the analyses using Sørensen's instead of

Jaccard's index to calculate dissimilarities (Supplementary Fig. 4) and when we excluded rare MOTUs (Supplementary Fig. 5).

β-total, β-richness and β-replacement changed steadily over time, given that segmented regressions did not reveal significant breakpoints for the relationships between these variables and mean age (all *P* > 0.05). Models with breakpoints also showed higher Bayesian information criterion values compared with the linear ones (Supplementary Table 4).

Discussion

Predicting ecosystem responses to disturbance events and environmental changes requires understanding of the mechanisms that govern community assembly during primary successions and thus modulate biodiversity. According to our results, compositional changes during successions after the retreat of glaciers are shaped by both the addition and the replacement of taxa. Within 483 yr after deglaciation, both mechanisms provided an overall similar contribution to compositional
<figure><img src="image_6.png" /><figcaption>**Fig. 2 | β-diversity components obtained with eDNA and traditional methods.** Left: components obtained with eDNA (purple, 771 comparisons between communities within the same glacier foreland). Right: components obtained with traditional methods (yellow, 102 comparisons). Boxplots indicate median (middle line), 25th and 75th percentiles (box), as well as ranges of 1.5×</figcaption></figure>

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

interquartile range (whiskers) and outliers (dots). Diamonds indicate the mean values. *P* values were obtained using two-sided randomization tests for paired samples assessing whether the differences between β-replacement and β-richness are significant.

<figure><img src="image_8.png" /><figcaption>**Fig. 3 | Patterns of β-diversity components over plant succession following glacier retreat.** Changes in total dissimilarity (β-total, grey), taxa replacement (β-replacement, blue) and taxa addition (β-richness, green) with age differences (a,c) and mean age (b,d) between compared communities. a,b, Patterns obtained with eDNA data (771 comparisons). c,d, Patterns obtained with traditional data (102 comparisons). Points represent observed values and shaded areas are 95%</figcaption></figure>

CIs obtained from Bayesian GLMMs (see Methods). Parameters with 95% CIs overlapping or not overlapping zero are represented with dashed or solid lines, respectively (see Table 1). The age differences between compared sites indicate how different the communities are, in terms of age, from each other. The mean age represents the average time since the glacier retreat of the compared communities.
<table><caption>Table 1 | Results of the Bayesian GLMMs assessing the effects of age differences between compared sites and time (mean age of compared sites) on the different β-diversity measures</caption><thead><tr><th rowspan="3">Dependent variable</th><th rowspan="3">Independent variables</th><th colspan="6">eDNA (N=771 comparisons)</th><th colspan="6">Traditional (N=102 comparisons)</th></tr><tr><th rowspan="2">B</th><th colspan="2">95% CI</th><th rowspan="2">R<sup>2</sup><sub>M</sub></th><th rowspan="2">R<sup>2</sup><sub>C</sub></th><th rowspan="2">B</th><th colspan="2">95% CI</th><th rowspan="2">R<sup>2</sup><sub>M</sub></th><th rowspan="2">R<sup>2</sup><sub>C</sub></th></tr><tr><th>Lower</th><th>Upper</th><th>Lower</th><th>Upper</th></tr></thead><tbody><tr><td rowspan="2">β-total</td><td>Age differences</td><td><strong>0.5</strong></td><td><strong>0.4</strong></td><td><strong>0.6</strong></td><td>0.1</td><td>0.3</td><td><strong>0.5</strong></td><td><strong>0.3</strong></td><td><strong>0.6</strong></td><td rowspan="2">0.3</td><td rowspan="2">0.9</td></tr><tr><td>Mean age</td><td><strong>-0.7</strong></td><td><strong>-0.8</strong></td><td><strong>-0.5</strong></td><td></td><td></td><td>-0.2</td><td>-0.5</td><td>0.03</td></tr><tr><td rowspan="2">β-richness</td><td>Age differences</td><td><strong>0.4</strong></td><td><strong>0.2</strong></td><td><strong>0.6</strong></td><td>0.05</td><td>0.6</td><td><strong>0.9</strong></td><td><strong>0.6</strong></td><td><strong>1.2</strong></td><td rowspan="2">0.3</td><td rowspan="2">0.8</td></tr><tr><td>Mean age</td><td><strong>-0.5</strong></td><td><strong>-0.7</strong></td><td><strong>-0.2</strong></td><td></td><td></td><td><strong>-0.7</strong></td><td><strong>-1.1</strong></td><td><strong>-0.4</strong></td></tr><tr><td rowspan="2">β-replacement</td><td>Age differences</td><td>-0.05</td><td>-0.2</td><td>0.1</td><td>0.01</td><td>0.6</td><td>-0.2</td><td>-0.4</td><td>0.01</td><td rowspan="2">0.2</td><td rowspan="2">0.8</td></tr><tr><td>Mean age</td><td>0.2</td><td>-0.08</td><td>0.4</td><td></td><td></td><td><strong>0.5</strong></td><td><strong>0.2</strong></td><td><strong>0.8</strong></td></tr></tbody></table>

Glacier identity and cross-site identity of sites involved in the comparisons were included as random factors. Parameters with 95% CIs not overlapping zero are highlighted in bold. R²<sub>M</sub>, marginal R; R²<sub>C</sub>, conditional R.

<figure><img src="image_9.png" /><figcaption>**Fig. 4 | Proportion of early colonizers, persisters and late colonizers in glacier forelands.** Percentages of taxa present exclusively in communities with ≤50 yr of succession ('Early colonizer'), taxa present exclusively in communities with >50 yr of succession ('Late colonizer') and taxa present in both groups of communities ('Persister'). Taxa categorization was performed within the chronosequence of each glacier foreland.</figcaption></figure>

differences in plant communities (Fig. 2). Nevertheless, their contribution to total β-diversity varied over time, supporting the hypothesis that the mechanisms driving succession after glacier retreat change over time<sup>22,47–49</sup>. Immediately after glacier retreat, richness differences contributed more to β-total than replacement, as expected in harsh environments<sup>27</sup>. This suggests an overall predominant role of taxa addition in early plant primary succession, with replacement becoming dominant after >50 yr following glacier retreat.

The more the communities differed in age, the more dissimilar they were in terms of composition. Furthermore, the dissimilarity between communities with strong age differences was mostly driven by β-richness (Fig. 3a–c and Table 1). This pattern matches the observed taxonomic accumulation from recently deglaciated terrains to late-successional stages<sup>13,47,50</sup>. Our results question studies advocating that severe environments are characterized by a constant initial floristic composition without changes in species composition over time resulting from the lack of establishment of additional species (autosuccession)<sup>12,14</sup>. In fact, only 16% of the taxa detected with eDNA and 35% of taxa detected with traditional sampling persisted after 50 yr of succession (Fig. 4). Such apparent incongruence could be explained by the ambiguity of the concept of 'severe' or harsh' environment. These terms apply to limiting conditions both linked to climate (conditions in high altitudes and/or latitudes) and to specific microhabitats within recently deglaciated terrains where soils are nutrient poor and geomorphological disturbances are frequent<sup>51</sup>. All these conditions

are typical of recently deglaciated terrains but can have distinct effects on communities. While the climatic and edaphic constraints faced by plants in Arctic and alpine environments are reported to reduce species replacement, sometimes leading to autosuccession<sup>12,14</sup>, the retreat of glaciers also results in changing physical conditions over time. The initial absence of soil and biota would constrain colonization to relatively fast-growing and opportunistic species adapted to alpine environments<sup>36,52</sup>, with the gradual addition of more dispersal-limited species<sup>32,53</sup>. Later, these species are followed by mountain specialists that are more competitive because they are shade-resistant and/or exploit nutrients efficiently<sup>24,34,36,51,54</sup>.

Total community dissimilarity was influenced not only by the differences in ages between communities but also by their mean age. In early-successional stages, β-total between communities was generally larger than between communities with similar age differences but being in late-successional stages. Thus, dissimilarity between sites decreases over time during succession. Considering that stochastically structured communities should exhibit divergent taxonomic compositions, our results suggest that community composition in early stages is strongly affected by initial conditions and/or stochastic processes (for example, priority effects, probabilistic dispersal and local extinction<sup>55</sup>), in agreement with temporal observations from studies using permanent plots<sup>9,22,56</sup>. Then, deterministic processes (for example, habitat filtering, competitive interactions) may drive more convergent community structures in late-successional stages.

When we compared communities in early-successional stages (having on average <50 yr), richness differences contributed more than replacement in determining the dissimilarity between communities (Fig. 3a–b). Immediately after glacier retreat, soils are generally nutrient poor and affected by surface instability, but early colonizers do not inhibit the establishment of new colonizing taxa¹. This may be explained either by neutral interactions (due to the predominant role of the environment<sup>57</sup> or of stochastic processes<sup>9,10</sup>) in these species-poor early stages, or by facilitative interactions¹ where the beneficiary species are not constraining the already established ones<sup>58</sup>. However, the importance of taxa addition quickly decreased over time and replacement became the dominant pattern for late-successional stages, suggesting increasing competition, as expected when resources, species richness and cover increase<sup>27</sup>. In late-successional stages, the stabilization of nutrients and terrains can allow new and more competitive alpine taxa to establish and replace pioneer species<sup>34</sup>. Such substitutions may occur either because early arriving modify the environment, making the conditions less suitable for themselves compared with other colonizers¹, or because later successional species outcompete the already established early species<sup>34,58</sup>.

Threshold dynamics have been proposed during the biotic colonization of glacier forelands, with a fast increase in alpha richness during the first 60 yr followed by a plateau and a decline in total
β-diversity<sup>22,59,60</sup>. However, in our study, the trends in β-diversity and its components during succession did not exhibit significant break-points (Supplementary Table 4), in agreement with Clements's view of successions as continuous trajectories without abrupt changes. Differences between this finding and previous results<sup>22</sup> might occur because temporal patterns obtained by focusing on one specific taxon (for example, plants) can differ from successional trajectories aggregating multiple diverse taxa<sup>35</sup> (such as plants, animals and microbes). It should be noted that our plant communities did not reach a stable point within the considered time frame. Even in our late-successional communities, β-total remained substantial and β-richness remained well above zero. Our sampling focused on terrains deglaciated since the Little Ice Age (mostly after 1850), whose ages may remain too young for stabilization of community composition. Longer time series (thousand years) would be required for a complete understanding of β-diversity changes and their drivers<sup>24,61</sup>, especially to identify whether and when β-diversity changes decelerate.

Our conclusions were highly consistent between eDNA and traditional sampling (Supplementary Fig. 1b) despite eDNA generally detecting fewer taxa per site. Similar to all sampling approaches, eDNA has its limitations. For instance, some taxa can remain undetected, it does not provide estimates of absolute biomass and taxonomic resolution is limited by the lack of complete reference databases<sup>43</sup>. Furthermore, marker amplification strongly depends on the amount of DNA present in environmental samples<sup>43</sup>, potentially favouring abundant taxa<sup>62-64</sup>. Amplifying rare taxa in soil samples can be particularly challenging, given the low diffusion rates of DNA in the soil and the presence of inhibitors<sup>44,65</sup>. Increasing the number of subplots per site, PCR replicates and sequencing depth can improve the detection of rare species. In most forelands, very recent sites without visible vegetation were not sampled by traditional methods but were sampled using eDNA. This might also contribute to the lower number of plants detected on average by eDNA. Another possible issue is that soil eDNA can detect past plant species<sup>44</sup>. However, the signal of these past species is weaker than that of current plant communities<sup>40,41,44</sup>, and eDNA inventories generally better match traditional inventories of current plant communities than those of older ones<sup>39</sup>. Despite these potential limitations, eDNA yielded temporal patterns extremely similar to those obtained by traditional methods, confirming that the former provides reliable diversity estimates, particularly for β-diversity<sup>66</sup>.

Some studies have questioned the robustness of the conclusions obtained from the analysis of chronosequences. Sites placed at similar distances from the glacier margin can differ from each other due to microhabitat conditions or the identity of first colonizers (priority effects), potentially following different trajectories<sup>38,67</sup>. Nevertheless, the analysis of temporal data from permanent plots yielded patterns of β-richness and β-replacement highly consistent with our results (see Methods and Extended data Fig. 1), confirming the robustness of conclusions drawn from the chronosequence approach. Moreover, the chronosequence approach assumes that there are no disturbances over succession and that terrain age is the community age, which is not always the case (see ref. 17). Even if we avoided sites clearly affected by geomorphological disturbances, glacier forelands are dynamic landscapes where disturbances might interact with temporal patterns according to their frequency and magnitude<sup>17</sup>. Indeed, despite retrieving general patterns of β-diversity changes, our data exhibited high variability, suggesting that other site-level environmental conditions, such as microclimate, soil properties and perturbations<sup>13,17,47</sup>, can affect the dissimilarity between sites. This calls for studies assessing how local drivers influence the contribution of taxa addition versus replacement, which will help to establish a general theory of succession that is lacking<sup>2,3</sup>.

The debate on processes shaping succession has persisted since the onset of community ecology, with both stochastic<sup>6</sup> and deterministic<sup>4</sup> processes pinpointed as key successional drivers during the

past century<sup>5</sup>. Our broad-scale study suggests that both processes play a fundamental role in community composition changes, with neutral and/or positive interactions dominating compositional variations in early-successional communities, and competition becoming more important in late-successional communities. Today, glaciers are retreating at an unprecedented rate and plant communities play a key-stone role in ecosystems developing after deglaciation<sup>13</sup>. The temporal changes in compositional drivers are expected to go beyond plant communities, affecting taxa interacting with plants through pollination, mutualism, herbivory or parasitism<sup>13,68</sup>. Understanding how β-diversity measures co-vary across different components of communities will be a key challenge in predicting the long-term consequences of climate change on ecosystems<sup>69</sup>.

## Methods

### eDNA sampling

In 46 glacier forelands (Fig. 1a), eDNA was captured from the soil in 1,255 plots during 2014 to 2020 using a previously published protocol<sup>70</sup>. In the sampled forelands, information on the dates of glacier retreat is available from the literature, remote sensing images and field surveys<sup>71</sup>. For each glacier foreland, the chronosequence approach<sup>38</sup> was used to select 3–17 sites along deglaciated terrains for which the date of each glacier retreat is known (Fig. 1b). Each site corresponded to a given age class; the number of sites depended on the number of documented positions of the glacier foreland available from the literature, and we tried to cover as much as possible the whole history of the retreat of each glacier<sup>71</sup>. We avoided sites clearly affected by geomorphological disturbance. Within each site, we sampled 2–10 plots (mean = 5, s.e. = 0.05; Supplementary Table 1). The plots within each site had similar distances to the glacier forefront and were, if possible, regularly spaced at distances of -20 m (Fig. 1b). At each plot, we collected five soil subsamples within 1 m (Fig. 1b) at a depth of 0–20 cm and pooled them to form a composite sample of -200 g per plot. We did not include soil litter and avoided roots, leaves and other large plant organs. Composite samples were homogenized; from each sample, we took 15 g of soil and desiccated it immediately in sterile boxes with 40 g of silica gel<sup>72</sup>. Before the collection of soil, all the sampling equipment underwent strict decontamination protocols (burned at >1,000 °C with a portable blow torch). In all countries, sampling was performed around the warmest season (for example, late July to early September in temperate areas of the Northern Hemisphere and February for temperate areas of the Southern Hemisphere).

Environmental DNA from the composite samples was extracted in a dedicated laboratory using the NucleoSpin soil mini kit (Macherey-Nagel), adding a preliminary step where the soil was mixed with 20 ml of phosphate buffer for 15 min<sup>73</sup>, and we eluted eDNA in 150 µl of elution buffer. To control for contamination in the extraction room, we included one extraction control every -10 samples (total: 101 extraction controls)<sup>43</sup>. We used the SperO1 primer pair<sup>74</sup> (forward: GGGCAATCCTGAGCCAA; reverse: CCATTGAGTCTCTGCACCTATC), which targets the P6 loop of the trnL intron in chloroplast DNA of Spermatophyta (seed plants). Amplicon size generally ranged from 10 to 220 bp (excluding the primers). We used reverse and forward primers that included 8-nucleotide-long tags on the 5' end. Each tag had at least five nucleotide differences from the others, thus allowing bioinformatic discrimination of PCR replicates after sequencing<sup>75</sup>. DNA extracts were randomized in 96-well plates together with extraction controls, bioinformatic blanks (that is, tagging-system controls), PCR-negative and positive controls (total across all plates: 291 blanks, 90 negative and 53 positive controls). In eDNA metabarcoding-based analyses, extraction and PCR-negative controls are pivotal for monitor contaminations, blanks allow identification of tag-jump issues and positive controls allow monitoring of potential cross-contamination of samples, as well as amplification and sequencing performance<sup>43</sup>. Positive controls consisted of a mock community composed of 16
non-tropical plant species belonging to 15 families (Taxaceae, Lamiaceae, Salicaceae, Polygonaceae, Betulaceae, Oleaceae, Pinaceae, Caprifoliaceae, Pinaceae, Aceraceae, Poaceae, Rosaceae, Brassicaceae, Geraniaceae, Ericaceae). Before amplification, we used quantitative PCR (qPCR) in a subset of samples to determine the optimal number of PCR cycles. We randomly selected 48 DNA samples and used 2 µl of undiluted or 1:10 diluted DNA and 1 µl of 1:1,000 diluted SYBR Green I nucleic acid gel stain (Invitrogen), with a real-time PCR thermal cycler set to standard mode. On the basis of qPCR results and for all samples, we performed 45 amplification cycles of 2 µl undiluted DNA in a 20 µl reaction volume with 10 µl of AmpliTaq Gold 360 Master Mix 2X (Applied Biosystems), 2 µl of primers mix (5 µM of each primer) and 0.16 µl of bovine serum albumin (Roche Diagnostic).

PCR amplifications of samples were performed in 384-well plates and consisted of an initial step of 10 min at 95 °C, followed by 45 cycles including 30 s denaturation at 95 °C, 30 s annealing at 52 °C, 60 s elongation at 72 °C and 7 min final elongation at 72 °C. All samples and controls underwent four PCR replicates⁷⁶. PCRs were performed in four distinct batches. All amplicons with a unique combination of forward and reverse tags within each batch were pooled. We used 5 µl aliquots of pooled amplicons to monitor the amplified fragment length and check for primer dimers using high-resolution capillary electrophoresis (QIAxcel Advanced System, Qiagen). Then, we purified six subsamples of the pooled amplicons using the MinElute PCR purification kit (Qiagen) following manufacturer protocol. Finally, we combined subsamples and sent them to Fasteris (Switzerland) where library preparation and sequencing were performed using the MetaFast protocol⁷⁵ and Illumina HiSeq platforms (paired-end approach, 2×150 bp), respectively.

The OBITools software suite²⁷ was used to perform the bioinformatic analyses of sequence data. First, forward and reverse reads were assembled with the 'illuminapairedend' programme, and the 'ngsfilter' programme was used to assign sequences with an alignment score >40 to the corresponding PCR replicate. Two mismatches on primers and zero mismatches on tags were allowed for this step. Then, we dereplicated sequences using the 'obiuniq' programme and filtered out those containing 'N' (representing that the nucleotide could not be determined) and/or with an unexpected sequence length (for example, <10 bp) and singletons. Subsequently, the 'obiclean' programme was used to keep sequences present in at least one PCR and that were at least twice as abundant as other related sequences differing by one base (hereafter, 'head sequences'). This step permitted removal of PCR and sequencing errors. At this point, sequences from different experiments were concatenated into one file and clustered at a threshold of 97% sequence similarity using the 'SUMACLUST' programme (https://git.metabarcoding.org/obitools/sumaclust/wikis/home). This threshold was selected on the basis of preliminary bioinformatics analyses, as it represents the threshold minimizing the risk of merging different species in the same MOTU while avoiding splitting of different sequences of the same species in different MOTUs⁷⁸. Finally, we performed a taxonomic assignment of cluster heads on the basis of the EMBL reference database (v.140). The reference database was built by carrying out an in silico PCR with the 'ecopcr' programme⁷⁹. Next, we assigned detected sequences to MOTUs using the 'ecotag' programme, following a previously described procedure⁷⁷. This programme matches each sequence in the dataset against the reference database and then uses the lowest common ancestor algorithm to identify the taxonomic level of the assignment (for example, genus, family, order)⁷⁷. To remove sequences detected at a low frequency that can be artefacts produced by PCR, contaminants and sequencing errors⁴,³,⁶⁶, we performed an additional filtering in R (v.4.3.1). Specifically, we discarded MOTUs with best identity <90% and detected less than eight times in all samples, which corresponds to the minimum number of reads that removed ≥99.99% of sequences detected in the blanks (that is, tag-jump errors). Then, we discarded MOTUs detected in only one sample, as they might represent

singletons⁸⁰; MOTUs detected in <2 PCR replicates of the same sample, as they can represent false positives⁷⁶; and MOTUs detected in more than one extraction or PCR-negative control, as they might represent contaminants⁴³. The complete codes and functions for bioinformatics and MOTU filtering are provided in Supplementary Codes 1–3. See Supplementary Table 5 for the number of sequences and MOTUs kept at each step of the procedure. Eleven of these MOTUs were removed because they were probably food contamination, as the corresponding families do not exist in the studied ecosystems and include species used as food (Supplementary Table 2).

## Traditional sampling

Even if the eDNA approach is emerging as a viable and reliable tool for sampling plant communities in soils⁴⁰,⁶²,⁶⁴,⁴⁵, some features of the method remain poorly understood (for example, the spatiotemporal scale⁴⁴). To confirm the reliability of the obtained eDNA patterns, we thus compared them with the ones obtained from morphological plant surveys. To this aim, we gathered observational inventories of plant communities from 13 of the 46 glacier forelands that were sampled with eDNA (Fig. 1a). In three cases (Carihuairazo, Pasquale and Rutor; Supplementary Table 1), floristic surveys were obtained from published studies performed in the same forelands where we collected eDNA; each of these studies sampled 4–5 dated sites per foreland. See refs. 20,54,81 for complete methodological aspects in these forelands. Literature data were complemented with inventories (Supplementary Table 3) collected along the chronosequences of 10 additional forelands (Supplementary Table 1). Along these chronosequences, 3 to 6 sampling sites corresponding to a given terrain age were sampled. Within each site, vascular plants were recorded in multiple plots (1–9, mean ± s.e.: 3.4 ± 0.4 plots per site, depending on site surface and geomorphological heterogeneity; Supplementary Table 1). Plots showed homogeneous altitudes, slopes and aspects throughout their surface. All plots within the same foreland had the same size. In most of the forelands, the plot surface was 25 m² (Supplementary Table 1), but in two forelands (Exploradores and Flaajokull), larger plots were used to better cover the geomorphological variability of the foreland (50–200 m²; Supplementary Table 1). Nevertheless, all results remain identical if the two forelands with larger plots were removed from the dataset (Supplementary Table 6). Plots were located in dated sites in the central portion of the foreland, that is, in front of the terminal part of the glacier tongue, avoiding disturbed areas (for example, those affected by glacial streams), as well as steep and unstable slopes. In each plot, we recorded every occurring vascular plant species; species that could not be identified in the field were collected and identified with the aid of identification keys for the local flora and with the aid of local experts, when necessary. Field sampling took into account all the vascular plants, that is, angiosperms, gymnosperms and pteridophytes (ferns, clubmosses and horsetails); however, to enable comparison with eDNA data, pteridophytes were not included in the analyses.

Overall, by combining literature data with original data, we gathered traditional data from forelands located in the Andes (n = 2), the Alps (n = 8), Iceland (n = 1), Nepal (n = 1) and Norway (n = 1). The overall dataset included 57 dated sites (time since glacier retreat ranging from 2 to 419 yr) where vascular plant communities were traditionally inventoried.

## β-diversity measures

For both methods, we combined the inventories obtained in all of the plots from the same dated site to recover comprehensive biodiversity inventories⁸². For eDNA and traditional data separately, we calculated the total dissimilarity (β-total) between communities within the same glacier foreland using the Jaccard's index based on presence/absence matrices. We then used the approach in ref. 46 to decompose β-total into β-replacement and β-richness, and to quantify the relative importance of these two processes. β-richness represents the richness
differences between compared communities associated with taxa
losses and gains, irrespective of nestedness<sup>46,83</sup>. We additionally parti-
tioned β-diversity into turnover (that is, replacement of some taxa by
others between communities) and nestedness (that is, richness differ-
ences where a community is a strict subset of a broader community)
following the approach in ref. 84. Despite retrieving patterns similar to
the ones obtained with the approach in ref. 46 for traditional data (Sup-
plementary Fig. 6), for eDNA data, the nestedness values were low and
the total dissimilarity was mostly driven by the turnover component
(see Supplementary Note 1 and Figs. 7 and 8). In fact, the method has
been found to underestimate richness differences and overestimate
the turnover component<sup>16,83,85</sup>, which may be specially marked when
assessing diversity with MOTUs (which is often the case for eDNA stud-
ies on a large geographic scale); the approach in ref. 46 appeared to be
less sensitive to these issues. β-diversity partitioning was performed
with the approach in ref. 46 through the ‘BAT’ R package<sup>86</sup>, while the
‘betapart’ package<sup>87</sup> was used for partitioning as proposed in ref. 84.

Seven of the 13 chronosequences sampled with traditional sam-
pling showed an important variability in the number of sampled plots
per site (Supplementary Table 1). We accounted for this unbalanced
sampling by using a subsampling procedure to calculate β-diversity
and its components between the sites of each chronosequence. First,
we determined the minimum number of sampled plots per site ‘N’.
Second, for sites having a number of plots >N, we randomly selected N
plots. Third, we calculated β-diversity between sites, only considering
the N randomly selected plots. We repeated this procedure 999 times
and calculated the mean value of β-diversity and its components for
each pair of compared communities. We ran the main Bayesian models
with these measures. Full details and the R code used are provided in
Supplementary Code 4.

Statistical analyses

First, we tested whether β-replacement and β-richness provide an
overall different contribution to the β-diversity of pairs of communi-
ties. Two-sample randomization tests for paired samples were used
to assess whether the observed β-replacement is significantly greater
or lower than β-richness. To do this, we first calculated the mean dif-
ference between β-replacement and β-richness for each site and then
compared it to values expected under randomness⁸⁸. Expected val-
ues were obtained by reshuffling the data 10,000 times across two
random groups and calculating the mean difference between groups
across permutations. For this analysis, we used the ‘EnvStats’ pack-
age, considering a two-sided alternative hypothesis. We also used this
approach to test the differences between turnover and nestedness
(Supplementary Note 1).

Subsequently, we used Bayesian GLMMs to assess how the different
β-diversity measures varied over succession. We considered two predic-
tor variables: the age differences between compared communities and
the mean age between them. We ran three GLMMs per sampling method;
each GLMM included a different measure of β-diversity as a dependent
variable (that is, β-total, β-replacement and β-richness differences).
The glacier identity and the identity of each site involved in the com-
parison were included as random intercepts (β ≥ age differences + mean
age + (1|glacier) + (1|site1) + (1|site2)). Age differences and mean ages
were log-transformed and then scaled (mean = 0, s.d. = 1) to allow
comparison of their estimated effects. Models were run assuming a
beta distribution for all the response variables. β-diversity variables
were rescaled to avoid fixed zeros and ones according to ref. 89
((value × (N − 1) + 0.5) / N; with N being the number of observations).
Three Markov chain Monte Carlo chains using 10,000 iterations and a
burn-in of 5,000 were run in the ‘brms’ R package⁹⁰. For all models,
ct was <1.01, indicating convergence. We interpreted a strong contribu-
tion of a predictor variable to the β-diversity measures if the 95% CIs of
a parameter’s posterior distribution did not overlap zero. The interac-
tion between age difference and mean age was tested, but it was not

relevant for any measure of β-diversity, so we kept the models without
interaction. In the main results, we show the random intercept models,
as models including glacier identity as both random intercepts and
slopes did not show a substantial decrease in the WAIC (widely applica-
ble information criterion) compared with models with random intercept
only (except for the β-richness and β-total models with traditional data;
Supplementary Table 7). All the effects remain identical if random slope
models were used instead of random intercept models (Supplementary
Table 8). To illustrate the main results (Fig. 3), we kept all the available
comparisons, but traditional data and eDNA data covered different
temporal ranges. eDNA data ranged from 1–475 yr and 2–350 yr for age
differences and mean ages, respectively, while traditional data ranged
from 4–397 yr and 7.5–281.5 yr for age differences and mean ages,
respectively. See Supplementary Fig. 9 (where eDNA comparisons hav-
ing age differences and mean ages outside of the ranges covered by
traditional data were removed) for a visualization of the results of the
main models within a common temporal range.

To assess threshold dynamics over time, we used segmented
regressions⁹¹,⁹². Specifically, we checked for the existence of thresh-
olds in the relationships of β-total, β-richness and β-replacement with
mean age. We used maximum likelihood to build linear mixed models
with one breakpoint ('segmented'⁹¹ package in R); glacier identity and
cross-site identity were included as random factors. Models with the
breakpoint were compared with linear mixed models without the
breakpoint on the basis of the Bayesian information criterion (BIC).
Simulations performed following the framework in ref. 92 confirmed
that this approach can detect different threshold typologies in datasets
with features similar to the ones analysed here (Supplementary Figs. 10
and 11). We used the packages ggplot2, cowplot, ggpubr and ggcorrplot
to create figures, as well as stringr, reshape2, reshape, rlist, plyr, dplyr
and reshape2 to format data.

Permanent plots to support conclusions drawn from the
chronosequence approach

Two approaches have been most commonly used to study plant succession in glacier forelands: the chronosequence approach and temporal data obtained from permanent plots. The chronosequence approach represents a space-for-time substitution to infer temporal changes in vegetation dynamics through contemporary spatial patterns. This approach is the most widely used, as it is the only one enabling us to reconstruct long-term (>100 yr) trends of ecosystem development⁶⁷,⁹³. In successions after glacier retreat, the chronosequence approach uses the assumption that communities established in terrains that are ice-free since longer times underwent a longer development time. However, the ability of chronosequences to reflect reliable temporal patterns has been questioned³⁸ (but see ref. 67). For instance, the method might mis-estimate the relative contribution of replacement and taxa addition if the landscape context strongly differs between sites along the chronosequences. Therefore, we extracted data obtained in permanent plots from ref. 94 to validate the conclusions obtained through our analysis of chronosequences. Ref. 94 sampled 12 permanent plots during the first decade after the deglaciation of two glaciers in the Alps: Goldbergkes (Austria) and Lenksteinferner (Italy). Vascular plants in the plots were sampled with traditional methods every 2 yr from 2005 to 2015. We used this compositional data from permanent plots to calculate and decompose β-diversity between communities using the approach in ref. 46. On average, taxa addition contributed to 85% of total β-diversity, whereas replacement contributed only to 15% (Extended Data Fig. 1a). Randomization tests for paired samples showed that β-richness and β-replacement were significantly different (*P* < 0.001), with a pattern similar to the one observed with chronosequences. Furthermore, the importance of addition decreases over time, even if the trend is less marked than for our main results (Extended Data Fig. 1b). Overall, the analysis of permanent plots showed findings highly consistent with the analysis of chronosequences.
**Reporting summary**

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

**Data availability**

Raw sequence data (Sper01 marker) generated using the protocols described in the 'Methods' section are deposited in the 'Sper01_raw_sequences.zip' folder available in Zenodo (https://zenodo.org/record/6620359#.Y8E10P6Z05d)⁹⁵. The data that support the findings of this study are provided as Supplementary Tables 1–3.

**Code availability**

Scripts for reproducing the results in this study are available as Supplementary codes.

Supplementary Code 1. Code reproducing bioinformatics steps.

Supplementary Code 2. Code reproducing taxonomic assignation.

Supplementary Code 3. R code for the MOTU filtering after bioinformatic analyses to remove sequences with best identity <90% and detected at a low frequency that can be artefacts produced by PCR, contaminants and sequencing errors.

Supplementary Code 4. R code to calculate beta-diversity and its components, run the main models and illustrate results.

Supplementary Code 5. R code to test the ability of our sampling design to detect breakpoints in segmented regressions.

**References**

1. Connell, J. H. & Slatyer, R. O. Mechanisms of succession in natural communities and their role in community stability and organization. *Am. Nat.* **111**, 1119–1144 (1977).

2. Prach, K. & Walker, L. R. *Comparative Plant Succession among Terrestrial Biomes of the World* (Cambridge Univ. Press, 2020).

3. Walker, L. R. & del Moral, R. *Primary Succession and Ecosystem Rehabilitation* (Cambridge Univ. Press, 2003).

4. Clements, F. E. *Plant Succession: An Analysis of the Development of Vegetation* (Carnegie Institution of Washington, 1916).

5. Pulsford, S. A., Lindenmayer, D. B. & Driscoll, D. A. A succession of theories: purging redundancy from disturbance theory. *Biol. Rev.* **91**, 148–167 (2016).

6. Gleason, H. A. The individualistic concept of the plant association. *Bull. Torrey Bot. Soc.* **53**, 7–26 (1926).

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

8. Bayle, A. et al. Local environmental context drives heterogeneity of early succession dynamics in alpine glacier forefields. *Biogeosciences* **20**, 1649–1669 (2023).

9. Li, S. et al. Convergence and divergence in a long-term old-field succession: the importance of spatial scale and species abundance. *Ecol. Lett.* **19**, 1101–1109 (2016).

10. Fukami, T., Martijn Bezemer, T., Mortimer, S. R. & Putten, W. H. Species divergence and trait convergence in experimental plant community assembly. *Ecol. Lett.* **8**, 1283–1290 (2005).

11. Marteinsdóttir, B., Thórhallsdóttir, T. E. & Svavarsdóttir, K. An experimental test of the relationship between small-scale topography and seedling establishment in primary succession. *Plant Ecol.* **214**, 1007–1015 (2013).

12. Matthews, J. A., Hill, J. L., Winkler, S., Owen, G. & Vater, A. E. Autosuccession in alpine vegetation: testing the concept on an altitudinal bioclimatic gradient, Jotunheimen, southern Norway. *CATENA* **170**, 169–182 (2018).

13. Ficetola, G. F. et al. Dynamics of ecological communities following current retreat of glaciers. *Annu. Rev. Ecol. Evol. Syst.* **52**, 405–426 (2021).

14. Svoboda, J. & Henry, G. H. R. Succession in marginal Arctic environments. *Arct. Alp. Res.* **19**, 373–384 (1987).

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

16. Bosson, J. B. et al. Future emergence of new ecosystems caused by glacial retreat. *Nature* **620**, 562–569 (2023).

17. Wojcik, R., Eichel, J., Bradley, J. A. & Benning, L. G. How allogenic factors affect succession in glacier forefields. *Earth Sci. Rev.* **218**, 103642 (2021).

18. Fischer, A., Fickert, T., Schwaizer, G., Patzelt, G. & Groß, G. Vegetation dynamics in Alpine glacier forelands tackled from space. *Sci. Rep.* **9**, 13918 (2019).

19. Körner, C. *Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems* (Springer, 2021).

20. Rosero, P. et al. Multi-taxa colonisation along the foreland of a vanishing equatorial glacier. *Ecography* **44**, 1010–1021 (2021).

21. Llambí, L. D. et al. Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier. *Front. Ecol. Evol.* **9**, 657755 (2021).

22. Hanusch, M., He, X., Ruiz-Hernández, V. & Junker, R. R. Succession comprises a sequence of threshold-induced community assembly processes towards multidiversity. *Commun. Biol.* **5**, 424 (2022).

23. Erschbamer, B., Niederfriniger Schlag, R. & Winkler, E. Colonization processes on a central Alpine glacier foreland. *J. Veg. Sci.* **19**, 855–862 (2008).

24. Losapio, G. et al. The consequences of glacier retreat are uneven between plant species. *Front. Ecol. Evol.* **8**, 616562 (2021).

25. Koffel, T., Boudsocq, S., Loeuille, N. & Daufresne, T. Facilitation- vs. competition-driven succession: the key role of resource-ratio. *Ecol. Lett.* **21**, 1010–1021 (2018).

26. Callaway, R. M. et al. Positive interactions among alpine plants increase with stress. *Nature* **417**, 844–848 (2002).

27. Bertness, M. D. & Callaway, R. Positive interactions in communities. *Trends Ecol. Evol.* **9**, 191–193 (1994).

28. Walker, L. R., Clarkson, B. D., Silvester, W. B. & Clarkson, B. R. Colonization dynamics and facilitative impacts of a nitrogen-fixing shrub in primary succession. *J. Veg. Sci.* **14**, 277–290 (2003).

29. Chapin, F. S., Walker, L. R., Fastie, C. L. & Sharman, L. C. Mechanisms of primary succession following deglaciation at Glacier Bay, Alaska. *Ecol. Monogr.* **64**, 149–175 (1994).

30. Gerla, D. J., Mooij, W. M. & Huisman, J. Photoinhibition and the assembly of light-limited phytoplankton communities. *Oikos* **120**, 359–368 (2011).

31. Losapio, G. et al. Network motifs involving both competition and facilitation predict biodiversity in alpine plant communities. *Proc. Natl Acad. Sci. USA* **118**, e2005759118 (2021).

32. Erschbamer, B. & Caccianiga, M. S. In *Progress in Botany* Vol. 78 (eds Cánovas, F. M. et al.) 259–284 (Springer, 2016).

33. Grime, J. P. *Plant Strategies, Vegetation Processes, and Ecosystem Properties* (Wiley, 2001).

34. Pérez, C. A. et al. Ecosystem development in short-term postglacial chronosequences: N and P limitation in glacier forelands from Santa Inés Island, Magellan Strait. *Austral Ecol.* **39**, 288–303 (2014).

35. Pothula, S. K. & Adams, B. J. Community assembly in the wake of glacial retreat: a meta-analysis. *Glob. Change Biol.* **28**, 6973–6991 (2022).

36. Anthelme, F., Carrasquer, I., Ceballos, J. L. & Peyre, G. Novel plant communities after glacial retreat in Colombia: (many) losses and (few) gains. *Alp. Bot.* **132**, 211–222 (2022).

37. Whittaker, R. H. Vegetation of the Siskiyou Mountains, Oregon and California. *Ecol. Monogr.* **30**, 279–338 (1960).

38. Johnson, E. A. & Miyanishi, K. Testing the assumptions of chronosequences in succession. *Ecol. Lett.* **11**, 419–431 (2008).
39. Ariza, M. et al. Plant biodiversity assessment through soil eDNA reflects temporal and local diversity. *Methods Ecol. Evol.* **14**, 415–430 (2023).

40. Yoccoz, N. G. et al. DNA from soil mirrors plant taxonomic and growth form diversity. *Mol. Ecol.* **21**, 3647–3655 (2012).

41. Foucher, A. et al. Persistence of environmental DNA in cultivated soils: implication of this memory effect for reconstructing the dynamics of land use and cover changes. *Sci. Rep.* **10**, 10502 (2020).

42. Taberlet, P., Coissac, E., Hajibabaei, M. & Rieseberg, L. H. Environmental DNA. *Mol. Ecol.* **21**, 1789–1793 (2012).

43. Zinger, L. et al. DNA metabarcoding—need for robust experimental designs to draw sound ecological conclusions. *Mol. Ecol.* **28**, 1857–1862 (2019).

44. Johnson, M. D. et al. Environmental DNA as an emerging tool in botanical research. *Am. J. Bot.* **110**, e16120 (2023).

45. Banerjee, P. et al. Environmental DNA analysis as an emerging non-destructive method for plant biodiversity monitoring: a review. *AoB Plants* **14**, plac031 (2022).

46. Carvalho, J. C., Cardoso, P. & Gomes, P. Determining the relative roles of species replacement and species richness differences in generating beta-diversity patterns: partitioning beta diversity. *Glob. Ecol. Biogeogr.* **21**, 760–771 (2012).

47. Raffl, C., Mallaun, M., Mayer, R. & Erschbamer, B. Vegetation succession pattern and diversity changes in a glacier valley, Central Alps, Austria. *Arct. Antarct. Alp. Res.* **38**, 421–428 (2006).

48. Tscherko, D., Hammesfahr, U., Zeltner, G., Kandeler, E. & Böcker, R. Plant succession and rhizosphere microbial communities in a recently deglaciated alpine terrain. *Basic Appl. Ecol.* **6**, 367–383 (2005).

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

50. Cauvy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat. *Nat. Ecol. Evol.* **3**, 1675–1685 (2019).

51. Zanzottera, M., Dalle Fratte, M., Caccianiga, M., Pierce, S. & Cerabolini, B. E. L. Community-level variation in plant functional traits and ecological strategies shapes habitat structure along succession gradients in alpine environment. *Community Ecol.* **21**, 55–65 (2020).

52. Anthelme, F., Cauvy-Fraunié, S., Francou, B., Cáceres, B. & Dangles, O. Living at the edge: increasing stress for plants 2–13 years after the retreat of a tropical glacier. *Front. Ecol. Evol.* **9**, 584872 (2021).

53. Erschbamer, B. & Mayer, R. Can successional species groups be discriminated based on their life history traits? A study from a glacier foreland in the Central Alps. *Plant Ecol. Divers.* **4**, 341–351 (2011).

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

55. Chase, J. M. & Myers, J. A. Disentangling the importance of ecological niches from stochastic processes across scales. *Phil. Trans. R. Soc. B* **366**, 2351–2363 (2011).

56. del Moral, R. Increasing deterministic control of primary succession on Mount St. Helens, Washington. *J. Veg. Sci.* **20**, 1145–1154 (2009).

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

58. Paterno, G. B., Siqueira Filho, J. A. & Ganade, G. Species-specific facilitation, ontogenetic shifts and consequences for plant community succession. *J. Veg. Sci.* **27**, 606–615 (2016).

59. Brambilla, M. & Gobbi, M. A century of chasing the ice: delayed colonisation of ice-free sites by ground beetles along glacier forelands in the Alps. *Ecography* **37**, 33–42 (2014).

60. Gobbi, M., Fontaneto, D. & De Bernardi, F. Influence of climate changes on animal communities in space and time: the case of spider assemblages along an alpine glacier foreland. *Glob. Change Biol.* **12**, 1985–1992 (2006).

61. Delgado-Baquerizo, M. et al. Changes in belowground biodiversity during ecosystem development. *Proc. Natl Acad. Sci. USA* **116**, 6891–6896 (2019).

62. Carrasco-Puga, G. et al. Revealing hidden plant diversity in arid environments. *Ecography* **44**, 98–111 (2021).

63. Hartvig, I., Kosawang, C., Kjær, E. D. & Nielsen, L. R. Detecting rare terrestrial orchids and associated plant communities from soil samples with eDNA methods. *Biodivers. Conserv.* **30**, 3879–3901 (2021).

64. Edwards, M. E. et al. Metabarcoding of modern soil DNA gives a highly local vegetation signal in Svalbard tundra. *Holocene* **28**, 2006–2016 (2018).

65. Wang, H., Qi, J., Xiao, D., Wang, Z. & Tian, K. A re-evaluation of dilution for eliminating PCR inhibition in soil DNA samples. *Soil Biol. Biochem.* **106**, 109–118 (2017).

66. Calderón-Sanou, I., Münkemüller, T., Boyer, F., Zinger, L. & Thuiller, W. From environmental DNA sequences to ecological conclusions: how strong is the influence of methodological choices? *J. Biogeogr.* **47**, 193–206 (2019).

67. Walker, L. R., Wardle, D. A., Bardgett, R. D. & Clarkson, B. D. The use of chronosequences in studies of ecological succession and soil development: chronosequences, succession and soil development. *J. Ecol.* **98**, 725–736 (2010).

68. Albrecht, M., Riesen, M. & Schmid, B. Plant-pollinator network assembly along the chronosequence of a glacier foreland. *Oikos* **119**, 1610–1624 (2010).

69. Urban, M. C. et al. Improving the forecast for biodiversity under climate change. *Science* **353**, aad8466 (2016).

70. Guerrieri, A. et al. Local climate modulates the development of soil nematode communities after glacier retreat. *Glob. Chang. Biol.* **30**, e17057 (2024)

71. Marta, S. et al. The retreat of mountain glaciers since the little ice age: a spatially explicit database. *Data* **6**, 107 (2021).

72. Guerrieri, A. et al. Effects of soil preservation for biodiversity monitoring using environmental DNA. *Mol. Ecol.* **30**, 3313–3325 (2021).

73. Taberlet, P. et al. Soil sampling and isolation of extracellular DNA from large amount of starting material suitable for metabarcoding studies. *Mol. Ecol.* **21**, 1816–1820 (2012).

74. Taberlet, P. et al. Power and limitations of the chloroplast trnL (UAA) intron for plant DNA barcoding. *Nucleic Acids Res.* **35**, e14 (2007).

75. Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. *Environmental DNA for Biodiversity Research and Monitoring* (Oxford Univ. Press, 2018).

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

77. Boyer, F. et al. OBITools: a unix-inspired software package for DNA metabarcoding. *Mol. Ecol. Resour.* **16**, 176–182 (2016).

78. Bonin, A., Guerrieri, A. & Ficetola, G. F. Optimal sequence similarity thresholds for clustering of molecular operational taxonomic units in DNA metabarcoding studies. *Mol. Ecol. Resour.* **23**, 368–381 (2023).

79. Ficetola, G. F. et al. An in silico approach for the evaluation of DNA barcodes. *BMC Genomics* **11**, 434 (2010).

80. Bálint, M. et al. Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes. *FEMS Microbiol. Rev.* **40**, 686–700 (2016).

81. Caccianiga, M., Luzzaro, A., Pierce, S., Ceriani, R. M. & Cerabolini, B. The functional basis of a primary succession resolved by CSR classification. *Oikos* **112**, 10–20 (2006).
82. Dickie, I. A. et al. Towards robust and repeatable sampling methods in eDNA based studies. *Mol. Ecol. Resour.* **18**, 940–952 (2018).

83. Baselga, A. & Leprieur, F. Comparing methods to separate components of beta diversity. *Methods Ecol. Evol.* **6**, 1069–1079 (2015).

84. Baselga, A. Partitioning the turnover and nestedness components of beta diversity: partitioning beta diversity. *Glob. Ecol. Biogeogr.* **19**, 134–143 (2010).

85. Legendre, P. & De Cáceres, M. Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. *Ecol. Lett.* **16**, 951–963 (2013).

86. Cardoso, P., Rigal, F. & Carvalho, J. C. BAT – Biodiversity Assessment Tools, an R package for the measurement and estimation of alpha and beta taxon, phylogenetic and functional diversity. *Methods Ecol. Evol.* **6**, 232–236 (2015).

87. Baselga, A. & Orme, C. D. L. betapart: an R package for the study of beta diversity: Betapart package. *Methods Ecol. Evol.* **3**, 808–812 (2012).

88. Manly, B. F. *J. Randomization, Bootstrap and Monte Carlo Methods in Biology* (Chapman and Hall/CRC, 2017).

89. Smithson, M. & Verkuilen, J. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. *Psychol. Methods* **11**, 54–71 (2006).

90. Bürkner, P.-C. brms: an R package for Bayesian multilevel models using Stan. *J. Stat. Softw.* **80**, 1–28 (2017).

91. Muggeo, V. M. R. segmented: an R package to fit regression models with broken-line relationships. *R News* **8**, 20–25 (2008).

92. Ficetola, G. F. & Denoël, M. Ecological thresholds: an assessment of methods to identify abrupt changes in species–habitat relationships. *Ecography* **32**, 1075–1084 (2009).

93. Poorter. L. et al. Multidimensional tropical forest recovery. *Science* **374**, 1370–1376 (2021).

94. Fickert, T. & Grüninger, F. High-speed colonization of bare ground – permanent plot studies on primary succession of plants in recently deglaciated glacier forelands. *Land Degrad. Dev.* **29**, 2668–2680 (2018).

95. Guerrieri, A., Bonin, A., Gielly, L. & Ficetola, F. Raw sequencing data for studying the colonization of soil communities after glacier retreat (2022). https://doi.org/10.5281/zenodo.7628236

## Acknowledgements

This study was funded by the European Research Council under the European Community's Horizon 2020 Programme, Grant Agreement no. 772284 (IceCommunities) to I.C., A.C., A.G., S.M., A.B., R.A., R.S.A., F.G., L.G., N.K., G.A.D., J.P., W.T., M.C. and G.F.F. This research was also funded by Biodiversa+, the European Biodiversity Partnership under the 2021–2022 BiodivProtect joint call for research proposals, co-funded by the European Commission (grant agreement no. 101052342 'PrioritlIce-Vanishing habitats: conservation priorities for

glacier-related biodiversity threatened by climate change') to I.C., R.A., W.T., M.C., M.G. and G.F.F. and with the funding organizations MUR and ANR. We thank R. Kaufmann, A. Guisan, K. Sieron and M. A. Morales-Martínez for help and discussions at various phases of this project.

## Author contributions

I.C., M.C. and G.F.F. conceived, developed and wrote the paper, with input from A.C., R. A., F.A., S.C.-F., M.G., A.R., A. Zerboni, P.T., J.P. and W.T.; I.C. performed the statistical analyses; A.G., S.M., A.B., F.G. and G.F.F. contributed to data preparation and curation; A.G., A.B. and L.G. performed laboratory analyses; A.G., S.M., A.B., R.A., F.A., R.S.A., P.A., P.A.G., S.C.-F., J.L.C.L., P.C., M.C.S., J.C., J.A.C.R., C.C., R.C.E., O.D., A.E., S.E., A.F., L.G., F.G., M.G, S.H., N.K., R.I.M., G.P., F.P., A.R., N.U., Y.Y., V.Z., A. Zerboni, A. Zimmer, G.A.D., J.P. M.C. and G.F.F. participated in sampling and the initial development of the study. All authors reviewed and/or provided input on the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Extended data is available for this paper at
https://doi.org/10.1038/s41477-023-01609-4.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41477-023-01609-4.

**Correspondence and requests for materials** should be addressed to Isabel Cantera.

**Peer review information** *Nature Plants* thanks Ingolf Kühn, Robert Junker and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at
www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2024

Isabel Cantera Ⓘ¹✉, Alexis Carteron Ⓘ¹², Alessia Guerrieri Ⓘ¹³, Silvio Marta Ⓘ¹⁴, Aurélie Bonin Ⓘ¹³, Roberto Ambrosini Ⓘ¹, Fabien Anthelme⁵, Roberto Sergio Azzoni¹,⁶, Peter Almond⁷, Pablo Alviz Gazitúa⁸, Sophie Cauvy-Fraunié⁹, Jorge Luis Ceballos Lievano¹⁰, Pritam Chand Ⓘ¹¹, Milap Chand Sharma¹², John Clague¹³, Justiniano Alejo Cochachín Rapre¹⁴,³³, Chiara Compostella⁶, Rolando Cruz Encarnación¹⁴, Olivier Dangles Ⓘ¹⁵, Andre Eger¹⁶, Sergey Erokhin¹⁷, Andrea Franzetti¹⁸, Ludovic Gielly¹⁹, Fabrizio Gili Ⓘ¹,²⁰, Mauro Gobbi²¹, Sigmund Hågvar²²,²³, Norine Khedim²⁴, Rosa Isela Meneses²⁵,²⁶, Gwendolyn Peyre²⁷, Francesca Pittino¹⁸, Antoine Rabatel Ⓘ²⁸, Nurai Urseitova¹⁷, Yan Yang²⁹, Vitalii Zaginaev³⁰, Andrea Zerboni Ⓘ⁶, Anaïs Zimmer Ⓘ³¹, Pierre Taberlet¹⁹,²³, Guglielmina Adele Diolaiuti¹, Jerome Poulenard¹⁵, Wilfried Thuiller¹⁹, Marco Caccianiga Ⓘ³² & Gentile Francesco Ficetola¹,¹⁹
¹Dipartimento di Scienze e Politiche Ambientali, Università degli Studi di Milano, Milano, Italy. ²Université de Toulouse, Ecole d'Ingénieurs de Purpan, UMR INRAE-INPT DYNAFOR, Toulouse, France. ³Argaly, Bâtiment CleanSpace, Sainte-Hélène-du-Lac, France. ⁴Institute of Geosciences and Earth Resources, CNR, Pisa, Italy. ⁵Laboratory AMAP, IRD, University of Montpellier, CIRAD, CNRS, INRA, Montpellier, France. ⁶Dipartimento di Scienze della Terra 'Ardito Desio', Milano, Italy. ⁷Department of Soil and Physical Sciences, Lincoln University, Lincoln, New Zealand.

⁸Departamento de Ciencias Biológicas y Biodiversidad, Universidad de Los Lagos, Osorno, Chile. ⁹INRAE, UR RIVERLY, Centre de Lyon-Villeurbanne, Villeurbanne, France. ¹⁰Instituto de Hidrología, Meteorología y Estudios Ambientales IDEAM, Bogotá, Colombia. ¹¹Department of Geography, School of Environment and Earth Sciences, Central University of Punjab, VPO- Ghudda, Bathinda, Punjab, India. ¹²Centre for the Study of Regional Development - School of Social Sciences, Jawaharlal Nehru University, New Delhi, India. ¹³Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada. ¹⁴Área de Evaluación de Glaciares y Lfagunas, Autoridad Nacional del Agua, Huaraz, Peru. ¹⁵CEFE, Univ Montpellier, CNRS, EPHE, IRD, Univ Paul Valéry Montpellier 3, Montpellier, France. ¹⁶Mannaki Whenua - Landcare Research, Soils and Landscapes, Lincoln, New Zealand. ¹⁷Institute of Water Problems and Hydro-Energy, Kyrgyz National Academy of Sciences, Bishkek, Kyrgyzstan. ¹⁸Department of Earth and Environmental Sciences (DISAT), - University of Milano-Bicocca, Milano, Italy. ¹⁹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, Grenoble, France. ²⁰Department of Life Sciences and Systems Biology, University of Turin, Turin, Italy.

²¹Research and Museum Collections Office, Climate and Ecology Unit, MUSE-Science Museum, Trento, Italy. ²²Faculty of Environmental Sciences and Natural Resource Management (INA), Norwegian University of Life Sciences, Ås, Norway. ²³UiT - The Arctic University of Norway, Tromsø Museum, Tromsø, Norway. ²⁴Univ. Savoie Mont Blanc, Univ. Grenoble Alpes, EDYTEM, Chambéry, France. ²⁵Herbario Nacional de Bolivia: La Paz, La Paz, Bolivia. ²⁶Universidad Católica del Norte, Antofagasta, Chile. ²⁷Department of Civil and Environmental Engineering, University of the Andes, Bogotá, Colombia. ²⁸Univ. Grenoble Alpes, CNRS, IRD, Grenoble-INP, Institut des Géosciences de l'Environnement (IGE, UMR 5001), Grenoble, France. ²⁹Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu, China. ³⁰Mountain Societies Research Institute, University of Central Asia, Bishkek, Kyrgyzstan. ³¹Department of Geography and the Environment, University of Texas at Austin, Austin, TX, USA. ³²Dipartimento di Bioscienze, Università degli Studi di Milano, Milano, Italy. ³³Deceased: Justiniano Alejo Cochachín Rapre. e-mail: isa_cantera@hotmail.com
a

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

b

<table><thead><tr><th rowspan="2">Dependent variable</th><th rowspan="2">Independent variables</th><th rowspan="2">B</th><th colspan="2">95C.I.</th><th colspan="2">85C.I.</th><th rowspan="2">R<sup>2</sup><sub>Conditional</sub></th><th rowspan="2">R<sup>2</sup><sub>Marginal</sub></th></tr><tr><th>lower</th><th>upper</th><th>lower</th><th>upper</th></tr></thead><tbody><tr><td rowspan="2">β-total</td><td>Age differences</td><td>0.8</td><td><strong>0.4</strong></td><td><strong>1</strong></td><td><strong>0.5</strong></td><td><strong>0.8</strong></td><td rowspan="2">0.9</td><td rowspan="2">0.8</td></tr><tr><td>Mean age</td><td>-0.4</td><td>-0.8</td><td>0.04</td><td><strong>-0.6</strong></td><td><strong>-0.1</strong></td></tr><tr><td rowspan="2">β-richness</td><td>Age differences</td><td>0.8</td><td><strong>0.5</strong></td><td><strong>1.1</strong></td><td><strong>0.6</strong></td><td><strong>1</strong></td><td rowspan="2">0.9</td><td rowspan="2">0.8</td></tr><tr><td>Mean age</td><td>-0.6</td><td><strong>-1.1</strong></td><td><strong>-0.05</strong></td><td><strong>-0.9</strong></td><td><strong>-0.3</strong></td></tr><tr><td rowspan="2">β-replacement</td><td>Age differences</td><td>-0.1</td><td>-0.5</td><td>0.4</td><td>-0.4</td><td>0.2</td><td rowspan="2">0.6</td><td rowspan="2">0.4</td></tr><tr><td>Mean age</td><td>0.5</td><td>-0.2</td><td>1.2</td><td><strong>0.05</strong></td><td><strong>1</strong></td></tr></tbody></table>

**Extended Data Fig. 1** | Patterns of β-diversity and its components over time measured with temporal data to support the conclusions based on the chronosequence approach. Data was obtained from Fickert & Grüninger⁹⁴, which sampled vascular plants with traditional methods in permanent plots during the first decade after the deglaciation of two glaciers in the Alps (N = 30 comparisons). a. β-diversity components. Boxplots indicate median (middle

line), 25th, 75th percentiles (box), as well as 1.5* interquartile range (whiskers) and outliers (dots). Diamonds indicate the mean values. b. Results of the Bayesian generalized mixed models assessing the effects of mean age and age differences between compared sites on the different β-diversity measures. Glacier identity and identity of sites involved in the comparisons were included as random factors. Parameters with 95% CI non-overlapping zero are highlighted in bold.
## Reporting Summary

Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our *Editorial Policies* and the *Editorial Policy Checklist*.

## Statistics

For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a | Confirmed

☐ ☑ The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement

☐ ☑ A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

☐ ☑ The statistical test(s) used AND whether they are one- or two-sided
*Only common tests should be described solely by name; describe more complex techniques in the Methods section.*

☐ ☑ A description of all covariates tested

☐ ☑ A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons

☐ ☑ A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

☐ ☑ For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and *P* value noted
*Give *P* values as exact values whenever suitable.*

☐ ☑ For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

☐ ☑ For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes

☑ ☐ Estimates of effect sizes (e.g. Cohen's *d*, Pearson's *r*), indicating how they were calculated

*Our web collection on statistics for biologists contains articles on many of the points above.*

## Software and code

Policy information about availability of computer code

Data collection | No software was used for data collection.

Data analysis
Sequence data were analyzed using the OBITools software suite (Boyer et al., 2016. doi:DOI: 10.1111/1755-0998.12428) and the SUMACLUST program (https://git.metabarcoding.org/obitools/sumaclust/wikis/home).
As reference database for genetic assignment, we used EMBL release 140 (https://www.ebi.ac.uk/about/news/service-news/release-140-enas-assembledannotated-sequences-now-available-0)
Statistical analyses were performed using R version 4.0.5 (2021-03-31)z.

R packages:

-BAT (2.8.1) and betapart (1.5.4) to measure and decompose beta-diversity.

-brms (2.15.0) for Bayesian generalized linear mixed models and performance (0.8.0) to obtain R² of the models.

-segmented (1.4.1) to build linear mixed models with breakpoints.

-ggplot2 (3.4.0); cowplot (1.1.1); ggpubr (0.4.0); and ggcorrplot (0.1.3) for figures.

-stringr (1.4.0); reshape2 (1.4.4); reshape (0.8.8); rlist (0.4.6.2); plyr (1.8.8); dplyr (1.1.3) to format data.

-EnvStats (2.8.1) test the differences on the contribution of β-richness and β-replacement, as well as between beturnover and nestedness

Detailed description of methods is available in the methods sections: "eDNA sampling", "β-diversity measures" and "Statistical analyses".

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.
# Data

Policy information about availability of data

All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

- Accession codes, unique identifiers, or web links for publicly available datasets

- A description of any restrictions on data availability

- For clinical datasets or third party data, please ensure that the statement adheres to our policy

Raw sequence data (Sper01 marker) generated using the protocols described in the "Methods" section is deposited in the "Sper01_raws.zip" folder available at the Zenodo link (https://zenodo.org/record/6620359#.Y8E1OP6ZO5d). The data that support the findings of this study are provided as supplementary tables (Tables S1-S3).

# Research involving human participants, their data, or biological material

Policy information about studies with human participants or human data. See also policy information about sex, gender (identity/presentation), and sexual orientation and race, ethnicity and racism.

Reporting on sex and gender not applicable

Reporting on race, ethnicity, or other socially relevant groupings not applicable

Population characteristics not applicable

Recruitment not applicable

Ethics oversight not applicable

Note that full information on the approval of the study protocol must also be provided in the manuscript.

# Field-specific reporting

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

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

Life sciences

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

Behavioural & social sciences

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

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see [nature.com/documents/nr-reporting-summary-flat.pdf](https://nature.com/documents/nr-reporting-summary-flat.pdf)

# Ecological, evolutionary & environmental sciences study design

All studies must disclose on these points even when the disclosure is negative.

**Study description** We used both environmental DNA and traditional surveys to quantify the relative importance of species addition versus replacement in total beta-diversity of vascular plant communities colonizing glaciers forelands.

**Research sample** We obtained data from 323 dated sites (eDNA sampling: 266 sites; Traditional sampling: 57) covering 1-483 years since glacier retreat across 46 glaciers around the world. 2-9 sites per glacier forelands (average: 5) were selected along the chronosequences of glacier retreat (i.e., the chronological sequence of specific geomorphological features along deglaciated terrains for which the date of glacier retreat is known). A site represented one plant community and correspond to an age class indicating the past position of the glacier at known dates.

**Sampling strategy** At each site, we established 1-20 plots regularly spaced for a more repeatable biodiversity estimate of the plant community. To characterize communities, 1255 plots and 258 plots were used for the eDNA and traditional analyses, respectively.

**Data collection**

**eDNA data:** At each plot, we collected five soil subsamples within one meter and at a depth of 0-20 cm and pooled them together to form a composite sample of ~200 g per plot. Composite samples were homogenized; from each sample we took 15 g of soil and desiccated it immediately in sterile boxes with 40 g of silica gel and stored until molecular analyses. Environmental DNA from the soil samples was extracted in a dedicated room using the NucleoSpin Soil Mini Kit (Macherey-Nagel), adding a preliminary step where the soil was mixed with 20 ml of phosphate buffer for 15 min and we eluted eDNA in 150 µl of elution buffer. To control for contamination in the extraction room, we included one extraction control every ~10 samples (total: 101 extraction controls). We used the Sper01 primer pair, which targets the P6 loop of the trnL intron in chloroplast DNA of Spermatophyta (seed plants).

**Traditional data:**

We gathered observational inventories of plant communities from 13 of the 46 glacier forelands that were sampled with eDNA. In three cases (Carihuairazo, Pasquale, Rutor) floristic surveys were obtained from published studies performed in the same forelands
where we collected eDNA; each of these studies sampled 4-5 dated sites per foreland. Literature data were complemented with unpublished data collected along the chronosequences of nine additional forelands. Within each site, vascular plants were recorded in multiple homogeneous plots within each site (1-9; mean ± SE: 3.2 ± 0.4 plots per site), depending on site surface and geomorphological heterogeneity. In most of the forelands, the plot surface was 25 m2 (Table S1). In two forelands (Exploradores and Flaajokull), larger plots were used to better cover the geomorphological variability of sites (50-200 m2; Table S1). Plots were located in dated sites, in the central portion of the foreland, i.e. in front of the terminal part of the glacier tongue, avoiding disturbed areas (e.g. those affected by glacial streams) as well as steep and unstable slopes. Plots showed homogeneous altitudes, slopes, and aspects throughout their surface. In each plot, we recorded every occurring vascular plant species; species that could not be identified on the field were collected and identified with the aid of identification keys for the local flora and with the aid of local experts when necessary. Overall, by combining literature data with original data we gathered traditional data from forelands located in the Andes (n=2), the Alps (n=8), Iceland (n=1), Nepal (n=1), and Norway (n=1). The overall dataset included 57 dated sites (time since glacier retreat ranging 1-419 years) where vascular plant communities were traditionally inventoried.

### Timing and spatial scale

From 2000 to 2020, we sampled 1255 soil samples in 46 glacier forelands representative of some of the main mountain chains of Europe, Asia, the Americas and Oceania and 258 plant vascular plant inventories located across 13 glacier forelands in the Andes (n=2), the Alps (n=8), Iceland (n=1), Nepal (n=1), and Norway (n=1).

### Data exclusions

We excluded sites impacted by recent landslides / fluvial erosion.

The eDNA from each plot was analysed in four replicated PCRs. We then removed MOTUs with best identity < 90%, detected in only one sample (as they represent singletons), MOTUs that were detected in <2 PCR replicates of the same sample, as they can represent false positives and MOTUs detected in less than 20 reads in all samples, which corresponds to the minimum number of reads that removed ≥99.99% of sequences detected in the blanks (i.e., tag-jump errors). The complete code and functions are provided in supplementary information (Supplementary Script 1). Finally, 11 MOTUs were removed from eDNA inventories because they were probably food contamination, as the families do not exist in the studied ecosystems and include species used as food. See Table S7 for the number of sequences and MOTUs kept at each step of the procedure.

To enable comparison with eDNA data, Pteridophytes were not included in the traditional data analyses.

### Reproducibility

The coordinates of sampled sites are available with the dataset, thus sampling is reproducible. The results presented in the manuscript were obtained by analysing the output of statistical models; all data and script are available and the results are thus fully reproducible.

### Data availability

Raw sequence data generated using the protocols described in the "Material and Methods" section is deposited on Zenodo https://zenodo.org/record/6620359#.Y8E1OP6ZO5d. The data that support the findings of this study are provided as supplementary materials (Tables S1-S3).

Scripts for reproducing the results in this study are available as Supplementary codes.

Supplementary code 1: Code reproducing bioinformatics steps.

Supplementary code 2: Code reproducing taxonomic assignation.

Supplementary code 3: R code for the MOTU filtering after bio-informatic analyses to remove sequences with best identity < 90% and detected at a low frequency that can be artefacts produced by PCR, contaminants, and sequencing errors.

Supplementary code 4: R code to calculate beta-diversity and its components, run the main models, and illustrate results.

Supplementary code 5: R code to test the ability of our sampling design to detect breakpoints in segmented regressions.

### Randomization

DNA extracts were randomized in 96-well plates together with extraction controls, bioinformatic blanks (i.e. tagging-system controls), PCR negative and positive controls (total across all plates: 291 blanks, 90 negative and 53 positive controls).

### Blinding

Site location was replaced by a code for molecular analyses. Therefore, plant inventories derived from environmental DNA data were obtained without prior knowledge of site location, dating and environmental variables.

Did the study involve field work?

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

Yes

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

No

## Field work, collection and transport

### Field conditions

In all countries, sampling was performed during the warmest season (e.g., late July-early September in temperate areas of the Northern hemisphere and February for temperate areas of the Southern hemisphere).

### Location

Samples were collected in 46 glacier forelands, covering five continents from the equator to polar regions. The country were they were sampled and the precise location (latitude, longitude) are available in Table S1.

### Access & import/export

The access to all sites, eDNA collection and samples export was legally permitted. When applicable, the study complies with access and benefit permits authorizing collection, transport, analysis and divulgation of DNA data used in this study.

### Disturbance

No disturbance was caused by the field work.
We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

Materials & experimental systems

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☑</td><td>Antibodies</td></tr><tr><td>☑</td><td>Eukaryotic cell lines</td></tr><tr><td>☑</td><td>Palaeontology and archaeology</td></tr><tr><td>☑</td><td>Animals and other organisms</td></tr><tr><td>☑</td><td>Clinical data</td></tr><tr><td>☑</td><td>Dual use research of concern</td></tr><tr><td>☐</td><td>Plants</td></tr></tbody></table>

Methods

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☑</td><td>ChIP-seq</td></tr><tr><td>☑</td><td>Flow cytometry</td></tr><tr><td>☑</td><td>MRI-based neuroimaging</td></tr></tbody></table>