The importance of species addition versus replacement varies over succession in plant communities after glacial retreat

doi:10.21203/rs.3.rs-2482972/v1

Download PDF

Article

The importance of species addition versus replacement varies over succession in plant communities after glacial retreat

https://doi.org/10.21203/rs.3.rs-2482972/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 17 Jan, 2024

Read the published version in Nature Plants →

Version 1

posted

You are reading this latest preprint version

Mechanisms underlying plant succession remain highly debated. A global quantification of the relative importance of species addition versus replacement is lacking due to the local scope of most studies. We quantified their role in the variation of plant communities colonizing the forelands of 46 retreating glaciers distributed worldwide, using both environmental DNA and traditional surveys. Both mechanisms concur in determining community changes over time but their relative importance varied over time along successions. Taxa addition predominated immediately after glacier retreat, as expected in harsh environments, while replacement became more important for late-successional communities. Those changes were aligned with total beta-diversity changes, which were larger between early successional communities than between late-successional communities (>50 years since glacier retreat). Despite the complexity of community assembly over plant succession, our 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.

Biological sciences/Ecology/Community ecology

Biological sciences/Ecology/Biodiversity

Ecological successions - how communities change or replace one another over time - have been a cornerstone of ecology since its inception^1,2. Primary successions, i.e. the development of ecosystems where a severe disturbance opens up large areas lacking most of 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 over a century since Clements’s⁴ work. Still, despite decades of work, the mechanisms that drive plant primary succession remain not fully understood^2,5. In the pioneering Clementsian deterministic view, succession occurs in a progressive, directional, and homogenous manner, with a relatively stable and predictable community structure (the climax stage) reached after some time without punctuated changes^4,5. This deterministic predictability was rapidly questioned and debated⁶. Over the last century, some studies showed homogenous responses among plant successions^7,8 while others showed that successions are not always convergent^9,10, barely reach an equilibrium³, and are largely determined by stochastic processes¹¹.

While there is no clear concession about succession being stochastic or deterministic, it has been demonstrated that over time during succession, the number of plant species increases over time and it is generally accompanied by changes in community composition^2,3 (but see¹²). Those 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 (i.e. pioneer species) and the addition of new species. In contrast, replacement involves the substitution over time of early colonizers by other species^1,13. The relative importance of these two key mechanisms on plant succession is difficult to quantify as it may vary over the succession time, but also in function of the physical context, like the local topographic and environmental characteristics^2,8,12,14. Thus, addressing the importance of each mechanism requires the assessment of plant successions in multiple environmental settings.

Ongoing climate change is dramatically accelerating the retreat of glaciers worldwide¹⁵, exposing new ground surfaces to the development of plant succession¹³. As such, glacial forelands provide ideal conditions for studying plant primary succession. Recently deglaciated terrains are harsh environments as vegetation is sparse, they are nutrient poor, and they often are isolated and difficult to reach. It has been demonstrated that succession in those harsh environments is primarily dominated by a gradual addition of species during ecosystem development^16,17, while replacement is often hypothesized to be weak or absent in harsh environments^2,12,14. Still, it is possible that the relative importance of species addition versus replacement changes through time, and the short time window covered by the studies could have overemphasized the importance of addition. In principle, species accumulation could result from the joint effect of dispersal-related processes (i.e. species require time to disperse in a new area) and facilitative interactions^{7,12, 18–20}, which have a key role in physically harsh environments^21–23. Plant species that initially colonize recently deglaciated terrains can modify the environment through the accumulation of nutrients and can create new micro-environmental conditions^{21,22, 24–26} that are more suitable to the establishment and development of subsequent colonizers¹. Nevertheless, competition and facilitation were found to jointly affect biodiversity in alpine plant communities²⁷. Their relative contributions vary along environmental gradients¹² and are thus expected to vary over time along the glacier foreland. As environmental conditions become less extreme along the succession sequence, productivity and nutrient availability increase, and thus the importance of competition is expected to increase²⁸. This may happen because, beyond a given level of complexity, species richness and/or density, nutrient availability, space, and light become limiting^3,20. Therefore, species addition is expected to dominate plant succession soon after glacier retreat, with decreasing importance along succession and replacement becoming more important through time. 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 beta diversity (β-diversity, i.e. compositional variation between assemblages²⁹) of plant communities colonizing 46 glacier forelands distributed worldwide (Fig. 1). Specifically, we decomposed β-diversity into a β-richness (i.e. the species richness difference between communities due to species gain/loss) and a β-replacement part (i.e. the true turnover in species composition, without affecting species richness) to test the relative importance of taxa addition/loss versus replacement, respectively. We tested whether β-diversity and the relative contributions of its components change in a predictable way across the globe and over time and whether those changes are progressive or show breakpoints¹⁸. Glacier forelands include chronosequences of progressively older landscape units the further away from the glacier front³⁰, thus we used this space-for-time substitution to encompass a wide range of temporal changes in communities (sites deglaciated from 1 to more than 400 years). Within the chronosequences, we collected environmental DNA (eDNA) from soil to obtain vascular plant inventories in 266 communities; i.e. assemblage of taxa inhabiting a site with a specific age since glacier retreat along the chronosequence. The method produces inventories matching well recent botanical inventories³¹ and enables the rapid assessment of communities over broad geographic scales and from remote areas³². Additionally, the plant patterns obtained with eDNA were compared with the ones obtained for a subset of forelands (Fig. 1) using traditional sampling (i.e. based on the morphological identification of species), as knowledge gaps in the correspondence of the patterns obtained with the different approaches still exist.

After data filtering and clustering, 519 molecular taxonomic units (MOTUs) of vascular plants were detected with eDNA metabarcoding across the 266 dated sites (Table S1), with 0–60 detected MOTUs per site (mean = 10.3, SE = 0.64) and 5-150 MOTUs per foreland (mean = 36, SE = 4.4). No vascular plant MOTUs were detected at 17 sites, all of which were less than 32 years old. Using traditional floristic surveys, we detected 365 vascular plant taxa across 57 dated sites from 13 glacier forelands (Table S2), with 0-133 (mean = 23, SE = 3) taxa per site and 15–162 taxa (mean = 48, SE = 10) per foreland. Traditional sampling did not observe any vascular plant species only in one site, which was deglaciated just for one year (Carihuairazo glacier). To assess the variation in community composition over time, as well as the relative contribution of taxa addition versus taxa replacement on plant successions in glacier forelands, we quantified and partitioned the total β-diversity between pairs of communities within the same chronosequence of a foreland (i.e. comparisons among communities from different glacier forelands were not assessed). This resulted in 771 and 102 comparisons of communities for the eDNA and the traditional data respectively. Total β-diversity (β-total) was partitioned into taxa replacement (β-replacement) and richness differences (β-richness) following Carvalho et al.’s approach³³, which allows the quantification of the taxa addition process irrespective of whether communities are nested or not^33,34 (see Methods section for more details). Overall, β-richness and β-replacement values 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 higher than β-replacement values for eDNA data (Wilcoxon paired test: V = 112464 and P < 0.01), but not for traditional data (Wilcoxon paired test: V = 2238, P = 0.9).

Bayesian generalized mixed models were used to assess how the β-diversity components varied over succession, considering two independent variables: “age differences” between communities, and their “mean age”. The age differences between compared sites 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 (e.g., early versus early or late versus late), while high values represent comparisons between communities with different successional stages. The mean age is the time since glacier retreat, averaged 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 β-diversity increased with age differences between compared communities (Fig. 3a, c; Table S3). Furthermore, comparisons between old communities generally showed a lower β-diversity; 95% credible intervals of this estimate were consistently negative for eDNA data, while slightly overlapped zero for traditional data (Fig. 3b, d; Table S3). β-richness increased with age differences between compared communities (Fig. 3a,c; Table S3) but, after taking into account age differences, β-richness was lower for comparisons involving old communities than for comparisons involving early successional communities, with a highly consistent pattern between eDNA and traditional data (Fig. 3b, d; Table S3). β-replacement was not strongly related to age differences (Fig. 3a, c; Table S3) and tended to increase with the mean age of communities; for this relationship credible intervals were negative for traditional data and slightly overlapped zero for eDNA data (Fig. 3b, d; Table S3). Repeating the analyses using Sørensen’s index instead of Jaccard’s index to calculate dissimilarities yielded consistent results (Fig. S1), confirming the robustness of our conclusions regardless of the used index. Re-running models after excluding the 80 MOTUs found in a single site (rare MOTUs) yielded patterns consistent with the ones obtained with all the detected MOTUs (Fig. S2).

Total β-diversity and its components of β-richness and β-replacement changed progressively and linearly over time, given that segmented regressions did not reveal significant breaking points for the relationships between those variables and mean age (all P > 0.05). Additionally, the Bayesian information criterion (BIC) values of the linear models were lower than the BIC of the models with breaking points (Table S4).

Predicting ecosystem responses to disturbance events and environmental changes requires understanding the mechanisms that govern community assembly during primary succession 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 400 years after deglaciation, both mechanisms showed a similar contribution to compositional differences in plant communities (Fig. 2). Nevertheless, their contribution to total β-diversity varied over time. Immediately after glacier retreat, richness differences contributed more to β-diversity than replacement, as expected in harsh environments²³. This suggests an overall predominant role of taxa addition in plant primary succession, but replacement becomes dominant after more than 50 years 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 rather than β-replacement (Fig. 3b and d, Table S3). This matches the observed pattern of plant taxonomic accumulation from recently deglaciated terrains to late-successional stages^13,35,36. 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 the establishment of additional species (autosuccession)^12,14. In fact, only 16% of the taxa detected with eDNA and 35% of taxa detected with traditional sampling persisted after 50 years of succession (Fig. S3). Such apparent incongruity could be explained by the ambiguity of the concept of “severe” or harsh” environment. These terms apply to limiting conditions both linked to regional climate (high altitude and/or latitude, and related climatic conditions) and occurring specifically on recently-deglaciated terrains, characterized by nutrient-poor soils and frequent geomorphological disturbances³⁸. While the first conditions were reported to reduce species replacement ultimately leading to autosuccession^12,14, the second ones can favor fast-growing species (ruderals²⁸), which often include ubiquitous species^39,40, but are subsequently replaced by species exhibiting more competitive strategies^20,38,41. “Severe” environmental conditions thus include a wide array of ecological factors and relative responses with multifaceted effects along succession.

Moreover, total community dissimilarity is not only affected by the differences in ages between communities, but also by their mean age. The total β-diversity between communities in early successional stages was generally larger than between communities having similar age differences but being in late successional stages. Thus, our results point out that dissimilarity among 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 (e.g., priority effects, probabilistic dispersal, and local extinction⁴²), while deterministic processes (e.g., habitat filtering, competitive interactions) may drive more convergent community structures in late successional stages, converging with temporal observations from studies using permanent plots^9,18,43.

When we compared communities in early successional stages (having on average less than ~ 50 years), richness differences contribute much more than replacement in determining the dissimilarity between communities (the credible intervals of the importance of the two components did not overlap for any method, Fig. 3a-d), while replacement becomes the dominant pattern for late-successional stages. This supports the hypothesis that the mechanisms driving succession after glacial retreat can change through time^18,35,44,45. Immediately after glacial retreat, soils are generally nutrient-poor and affected by physical disturbances 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⁴⁶ or of stochastic processes^9,10) in these species-poor early stages or by facilitative interactions¹, where the beneficiary species are not constraining the already established species⁴⁷. However, the importance of taxa addition quickly decreased over the succession sequence, suggesting an increasing competition, as expected when environmental harshness decreases and species richness and cover increases²³. In late successional stages, the stabilization of resources and terrains can reduce the strength of environmental filtering, and the most competitive taxa may replace the least competitive ones. Such substitutions are likely driven by biotic interactions, either because early arrivers modify the environment making the conditions suitable for other colonizers but less suitable for them¹ or because later successional species exploit similar resources and outcompete the already established early species⁴⁷.

We found a consistent temporal trend of beta-diversity and its components across 46 chronosequences around the world, suggesting common features in the changes in the assembly rules underlying plant successions, despite existing variability across forelands. Local and regional environmental conditions are known to strongly influence the trajectory and rate of plant successions^{12,14,35,46,48}. Forthcoming studies using a global dataset covering a wide geographical scale, such as the one used here, will allow assessing the local and regional drivers that shape successional trajectories. The eDNA metabarcoding method enables rapid, cost-efficient, and standardized biodiversity assessment over broad geographical and taxonomic scales, yielding data that would have been challenging to assemble with traditional methods. Like all sampling approaches, eDNA has some limitations such as missing the rarest taxa, not providing estimates of absolute biomass, and limited taxonomic resolution⁴⁹. Nevertheless, eDNA yielded temporal patterns of β-diversity strongly consistent with traditional methods, suggesting that our results are robust to these methodological limitations.

Threshold dynamics have been proposed during the biotic colonization of glacier forelands, with a fast increase in community alpha richness during the first 60 years followed by a plateau and a decline in β-diversity^18,41. However, in our study, the global trends of β-diversity and its components during succession were characterized by linear relationships over time without significant breakpoints (Table S4), in agreement with Clements's views of successions as progressive, linear changes. It should be noted that our plant communities did not reach a stable point in the considered time frame. Even in our late successional communities, the total β-diversity remained substantial and β-richness remained well above zero. Our sampling was focused on terrains deglaciated since the Little Ice Age (LIA, after ca. 1850 for the majority of forelands), whose ages may remain too young for observing a stabilization of the composition of plant communities. Longer time series (thousand years) would be required for a complete understanding of β-diversity changes and their drivers^20,50, especially to identify whether and when β-diversity changes stop.

The robustness of conclusions about successions obtained from the analysis of chronosequences has been questioned, given that sites placed at different distances from the glacier margin can differ from each other due to microhabitat conditions or the identity of first colonizers (priority effects) and might thus potentially follow different trajectories^30,51. Nevertheless, the analysis of temporal data from permanent plots yielded temporal patterns of β-richness and β-replacement highly consistent with our results (Supplementary Note1, Fig. S4, Table S5), confirming the robustness of conclusions from the chronosequence approach. Importantly, the combination of the chronosequence approach and eDNA sampling allowed readdressing old ecological questions left behind due to the lack of properly standardized and comprehensive datasets. Even though we find common successional trends among different regions of the world, the temporal patterns of total β-diversity and its components highlight the complexity of plant successions after glacial retreat, where the drivers of community assembly changed over time. The debate on processes shaping succession continues since the onset of community ecology, with both stochastic⁶ and deterministic⁴ processes pinpointed as key drivers during the last century⁵. Our broad-scale study suggests that both processes have a fundamental role in community composition changes, with neutral and/or positive interactions dominating compositional variations in early successional communities, while competition becomes more important in late successional communities.

Today, glaciers are retreating at an unprecedented rate, and plant communities have an important role in ecosystems developing after deglaciation¹³. The temporal changes of compositional drivers are expected to go beyond plant communities, thus affecting taxa interacting with plants through pollination, mutualism, herbivory, or parasitism^13,52. 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⁵³.

eDNA sampling

In 46 glacier forelands (Fig. 1), 1255 soil samples were collected from 2014 to 2020 to capture eDNA. In the sampled forelands, information on the dates of glacier retreat is available from the literature, remote sensing images, and field surveys⁵⁴. For each glacier foreland, the chronosequence approach³⁰ was used to select 3–17 sites along deglaciated terrains for which the date of each glacier retreat is known. 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 ⁵⁴. Each site corresponded generally to the line delimiting the position of the glacier’s front on a given date. Along each line, 2–10 (mean = 5, SE = 0.8, Table S6) regularly spaced sampling plots (if possible, regularly spaced at distances of 20 m) were located; in most cases plots belonging to the same site were located approximately at the same distance from the glacier front. 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. 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⁵⁵. Before the collection of each sample, all the sampling equipment underwent strict decontamination protocols (burned at > 1,000°C with a portable blow torch). 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).

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). 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⁵⁸. 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). In eDNA metabarcoding-based analyses, extraction and PCR negative controls are pivotal to monitor contaminations, blanks allow identification of tag-jump issues, and positive controls allow monitoring of potential cross-contamination of samples, amplification and sequencing performance⁴⁹. 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). Prior to 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. Based on 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 the manufacturer's 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, 2x150 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 program and the ngsfilter program 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 program and filtered out those containing “N” and/or with an unexpected sequence length (e.g. <10 bp) and singletons. Subsequently, the obiclean program 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 to remove 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 program (https://git.metabarcoding.org/obitools/sumaclust/wikis/home). This threshold was selected based on preliminary bioinformatics analyses as it represents the threshold minimizing the risk of merging different species in the same MOTU while avoiding to split different sequences of the same species in different MOTUs⁶¹. Finally, we performed a taxonomic assignment of cluster heads based on the EMBL reference database (version 140). The reference database was built by carrying out an in-silico PCR with the ecopcr program⁶². Next, we assigned detected sequences to molecular operational taxonomic units (MOTUs) using the ecotag program, following the procedure described in Boyer et al.⁶⁰. This program 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 (e.g. genus, family, order)⁶⁰. We then performed additional filtering in R (version 4.0) to remove contaminants and tag-jump errors based on sequences detected in controls and blanks⁶³. Specifically, we discarded MOTUs with best identity < 90%, detected less than eight times in the overall data set, or detected in at least two extractions or PCR negative controls⁴⁹. We also removed MOTUs detected in only one sample, as they represent singletons, and MOTUs that were never found in at least two replicated PCRs per sample, as they are possible contaminants⁵⁹. See Table S7 for the number of sequences kept at each step in the analyses. Eleven MOTUs were removed because they were probably food contamination, as the families do not exist in the studied ecosystems and include species used as food.

Traditional sampling

We gathered plant vascular plant inventories from 13 glacier forelands (Fig. 1) in the Andes (n = 2), the Alps (n = 8), Iceland (n = 1), Nepal (n = 1), and Norway (n = 1). Along the chronosequences of those glaciers, 3 to 6 sampling sites corresponding to a given terrain age were sampled (one site per terrain age per glacier), resulting in 57 dated sites (time since glacier retreat ranging from 1-419 years) where vascular plant communities were inventoried. For the Carihuairazo glacier (Ecuadorian Andes) we used the inventories available from Rosero et al.¹⁶’s work. For the other forelands, we recorded vascular plants in 1–20 homogeneous plots within each site (Table S6), depending on site surface and geomorphological heterogeneity. Plots were sampled from 2000 to 2019 and their surface ranged from 100 to 200 m². Plot size was assessed for each glacier foreland by sampling increasing surfaces and plotting the number of species detected against the cumulative area until a plateau is reached. The resulting plot sizes may vary as a function of community structure and diversity as well as the geomorphological setting of the site. Plot surface was constant within each chronosequence unless geomorphological constraints forced otherwise. Plots were located in the central portion of the line delimiting the position of the glacier’s front on a given date, 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 seven chronosequences, we had high variability in the number of sampled plots per site (Table S6). We accounted for this unbalanced sampling by performing a subsampling procedure to calculate β-diversity and its components among the sites of each chronosequence. For each heterogeneously sampled chronosequence, we selected the minimum number of sampled plots obtained among the different sites. For the sites having a higher number of plots, we randomly selected the minimum number of plots per chronosequence to compare sites sampled with the same number of plots. We repeated this procedure 999 times and then we calculated the mean value of β-diversity and its components for each pair of compared communities.

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 local flora and with the aid of local experts when necessary. The ground cover value was visually estimated for each species with a resolution of 5% (i.e. assigning cover values of < 5%, 5%, 10%, 15%, and so on), together with the overall cover of vascular plants, rock outcrops, and debris (i.e. loose rocky material of any grain size). Field sampling took into account all the vascular plants, i.e. Angiosperms, Gymnosperms, and Pteridophytes (ferns, clubmosses, and horsetails) but, to enable comparison with eDNA data, Pteridophytes were not included in the further analyses.

β-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 β-diversity between communities within the same glacier foreland using the Jaccard’s index based on presence/absence matrices. The comparisons between communities from different glacier forelands were not included in our analyses. Total β-diversity (or dissimilarity) quantifies the variation in community composition between two communities²⁹. Two ecological mechanisms may explain compositional dissimilarities among communities during succession: the replacement of taxa (substitution of taxa in one community by different taxa in another community) and richness differences due to taxa gain and loss [a community has more (or less) richness than the other, e.g. because of the addition (or removal) of a taxon]. These processes contribute jointly to the compositional differentiation of communities but their relative importance may vary along environmental gradients. The Carvalho et al.³³’s approach was used through the BAT R package⁶⁵ to decompose dissimilarity into β-replacement and β-richness and to quantify the relative importance of those two processes. β-richness represents the richness differences between compared communities associated with taxa losses and gains, irrespective of nestedness^33,34. We additionally partitioned β-diversity into turnover (i.e. replacement of some taxa by others between communities) and nestedness (i.e. richness differences where a community is a strict subset of a broader community) following Baselga et al.⁶⁶’s approach. Despite retrieving similar patterns than the ones obtained with Carvalho et al.³³’s approach for traditional data, for eDNA data the nestedness values were low and the total dissimilarity was mostly driven by the turnover component (see Supplementary Note 2 and Fig. S5-S7). In fact, the method has been found to underestimate nestedness and overestimate the turnover component^33,34,67, which may be specially marked when assessing diversity with MOTUs (which is often the case for eDNA studies on a large geographic scale) and Carvalho’s approach appeared to be less sensitive to these issues.

Statistical analyses

We used Bayesian generalized linear mixed models (GLMMs) to assess how the different β-diversity measures varied through succession. We considered two 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 (i.e. Total β-diversity, β-replacement, and β-richness differences). The glacier identity, and the identity of each site involved in the comparison, were included as random intercepts to take into account the non-independence of sites within the same glacier, and multiple comparisons involving the same site [β ~ age differences + mean age + ( 1|glacier) + (1|site1) + (1|site2)]. Age differences and mean ages were log-transformed and then scaled (mean = 0, SD = 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 Smithson & Verkuilen⁶⁸[(value * (N-1) + 0.5 )/ N; with N = 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, \(\widehat{c}\) was < 1.01, indicating convergence. We interpreted a significant contribution of a predictor variable to the β-diversity measures if the 95% credible interval of a parameter’s posterior distribution did not include zero. The interaction 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. We also tested preliminary mixed models including glacier identity as random slopes. However, these models did not always show a substantial decrease in the WAIC (Widely Applicable Information Criterion) than the models with random intercept (Table S8), thus we retained the random intercept models.

To assess threshold dynamics, we used segmented regressions^70,71. Specifically, we checked the existence of thresholds in the relationships of total β-diversity, β-richness, and β-replacement with mean age. We used maximum likelihood to build linear mixed models with one breakpoint (segmented⁷⁰ package in R) and glacier identity and cross-site identity as random factors. Models with the breakpoint were compared with linear mixed models without the breakpoint based on the Bayesian Information Criterion (BIC).

Data availability

Raw sequence data generated using 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.

Connell, J. H. & Slatyer, R. O. Mechanisms of Succession in Natural Communities and Their Role in Community Stability and Organization. The American Naturalist 111, 1119–1144 (1977).
Prach, K. & Walker, L. R. Comparative Plant Succession among Terrestrial Biomes of the World. (Cambridge University Press, 2020). doi:10.1017/9781108561167.
Walker, L. R. & del Moral, R. Primary Succession and Ecosystem Rehabilitation. (Cambridge University Press, 2003). doi:10.1017/CBO9780511615078.
Clements, F. E. Plant succession: an analysis of the development of vegetation. (Carnegie Institution of Washington Publication Sciences, 1916).
Pulsford, S. A., Lindenmayer, D. B. & Driscoll, D. A. A succession of theories: purging redundancy from disturbance theory: Purging redundancy from disturbance theory. Biol Rev 91, 148–167 (2016).
Gleason, H. A. The Individualistic Concept of the Plant Association. Bulletin of the Torrey Botanical Club 53, 7–26 (1926).
Zimmer, A. et al. Time lag between glacial retreat and upward migration alters tropical alpine communities. Perspectives in Plant Ecology, Evolution and Systematics 30, 89–102 (2018).
Bayle, A. et al. Local environmental context drives heterogeneity of early succession dynamics in alpine glacier forefields. EGUsphere 2022, 1–33 (2022).
Li, S. et al. Convergence and divergence in a long-term old-field succession: the importance of spatial scale and species abundance. Ecology Letters 19, 1101–1109 (2016).
Fukami, T., Martijn Bezemer, T., Mortimer, S. R. & Putten, W. H. Species divergence and trait convergence in experimental plant community assembly. Ecol Letters 8, 1283–1290 (2005).
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).
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).
Ficetola, G. F. et al. Dynamics of Ecological Communities Following Current Retreat of Glaciers. Annu. Rev. Ecol. Evol. Syst. 52, 405–426 (2021).
Svoboda, J. & Henry, G. H. R. Succession in Marginal Arctic Environments. Arctic and Alpine Research 19, 373 (1987).
Zemp, M. et al. Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. Nature 568, 382–386 (2019).
Rosero, P. et al. Multi-taxa colonisation along the foreland of a vanishing equatorial glacier. Ecography 44, 1010–1021 (2021).
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).
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).
Erschbamer, B., Niederfriniger Schlag, R. & Winkler, E. Colonization processes on a central Alpine glacier foreland. Journal of Vegetation Science 19, 855–862 (2008).
Losapio, G. et al. The Consequences of Glacier Retreat Are Uneven Between Plant Species. Front. Ecol. Evol. 8, 616562 (2021).
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).
Callaway, R. M. et al. Positive interactions among alpine plants increase with stress. Nature 417, 844–848 (2002).
Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends in Ecology & Evolution 9, 191–193 (1994).
Walker, L. R., Clarkson, B. D., Silvester, W. B. & Clarkson, R., Beverley. Colonization dynamics and facilitative impacts of a nitrogen-fixing shrub in primary succession. Journal of Vegetation Science 14, 277–290 (2003).
Chapin, F. S., Walker, L. R., Fastie, C. L. & Sharman, L. C. Mechanisms of Primary Succession Following Deglaciation at Glacier Bay, Alaska. Ecological Monographs 64, 149–175 (1994).
Gerla, D. J., Mooij, W. M. & Huisman, J. Photoinhibition and the assembly of light-limited phytoplankton communities. Oikos 120, 359–368 (2011).
Losapio, G. et al. Network motifs involving both competition and facilitation predict biodiversity in alpine plant communities. Proc. Natl. Acad. Sci. U.S.A. 118, e2005759118 (2021).
Grime, J. P. Plant strategies, vegetation processes, and ecosystem properties. (Wiley, 2001).
Whittaker, R. H. Vegetation of the Siskiyou Mountains, Oregon and California. Ecological Monographs 30, 279–338 (1960).
Johnson, E. A. & Miyanishi, K. Testing the assumptions of chronosequences in succession. Ecol Letters 11, 419–431 (2008).
Ariza, M. et al. Plant biodiversity assessment through soil eDNA reflects temporal and local diversity. Methods Ecol Evol 2041–210X.13865 (2022) doi:10.1111/2041-210X.13865.
Taberlet, P., Coissac, E., Hajibabaei, M. & Rieseberg, L. H. Environmental DNA. Molecular ecology 21, 1789–1793 (2012).
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. Global Ecology and Biogeography 21, 760–771 (2012).
Baselga, A. & Leprieur, F. Comparing methods to separate components of beta diversity. Methods Ecol Evol 6, 1069–1079 (2015).
Raffl, C., Mallaun, M., Mayer, R. & Erschbamer, B. Vegetation Succession Pattern and Diversity Changes in a Glacier Valley, Central Alps, Austria. Arctic, Antarctic, and Alpine Research 38, 421–428 (2006).
Cauvy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat. Nature Ecology & Evolution 3, 1675–1685 (2019).
Egler, F. E. Vegetation science concepts I. Initial floristic composition, a factor in old-field vegetation development with 2 figs. Vegetatio 4, 412–417 (2004).
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 Ecology 21, 55–65 (2020).
Anthelme, F., Carrasquer, I., Ceballos, J. L. & Peyre, G. Novel plant communities after glacial retreat in Colombia: (many) losses and (few) gains. Alp Botany 132, 211–222 (2022).
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).
Gobbi, M. et al. Plant adaptive responses during primary succession are associated with functional adaptations in ground beetles on deglaciated terrain. Community Ecology 11, 223–231 (2010).
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).
del Moral, R. Increasing deterministic control of primary succession on Mount St. Helens, Washington. Journal of Vegetation Science 20, 1145–1154 (2009).
Tscherko, D., Hammesfahr, U., Zeltner, G., Kandeler, E. & Böcker, R. Plant succession and rhizosphere microbial communities in a recently deglaciated alpine terrain. Basic and Applied Ecology 6, 367–383 (2005).
Kaufmann, R. Invertebrate succession on an alpine glacier foreland. Ecology 82, 2261–2278 (2001).
Matthews, J. A. The ecology of recently-deglaciated terrain. A geoecological approach to glacier forelands and primary succession. (Cambridge University Press, 1992).
Paterno, G. B., Siqueira Filho, J. A. & Ganade, G. Species-specific facilitation, ontogenetic shifts and consequences for plant community succession. Journal of Vegetation Science 27, 606–615 (2016).
Poorter, L. et al. Wet and dry tropical forests show opposite successional pathways in wood density but converge over time. Nat Ecol Evol 3, 928–934 (2019).
Zinger, L. et al. DNA metabarcoding—Need for robust experimental designs to draw sound ecological conclusions. Molecular Ecology 28, 1857–1862 (2019).
Delgado-Baquerizo, M. et al. Changes in belowground biodiversity during ecosystem development. Proc. Natl. Acad. Sci. U.S.A. 116, 6891–6896 (2019).
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. Journal of Ecology 98, 725–736 (2010).
Albrecht, M., Riesen, M. & Schmid, B. Plant-pollinator network assembly along the chronosequence of a glacier foreland. Oikos 119, 1610–1624 (2010).
Urban, M. C. et al. Improving the forecast for biodiversity under climate change. Science 353, aad8466 (2016).
Marta, S. et al. The Retreat of Mountain Glaciers since the Little Ice Age: A Spatially Explicit Database. Data 6, (2021).
Guerrieri, A. et al. Effects of soil preservation for biodiversity monitoring using environmental DNA. Molecular Ecology 30, 3313–3325 (2021).
Taberlet, P. et al. Soil sampling and isolation of extracellular DNA from large amount of starting material suitable for metabarcoding studies. Molecular Ecology 21, 1816–1820 (2012).
Taberlet, P. et al. Power and limitations of the chloroplast trnL (UAA) intron for plant DNA barcoding. Nucleic Acids Research 35, e14–e14 (2007).
Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. Environmental DNA for biodiversity research and monitoring. (Oxford University Press, 2018).
Ficetola, G. F. et al. Replication levels, false presences and the estimation of the presence/absence from eDNA metabarcoding data. Molecular Ecology Resources 15, 543–556 (2015).
Boyer, F. et al. OBITools: a unix-inspired software package for DNA metabarcoding. Molecular ecology resources 16, 176–182 (2016).
Bonin, A., Guerrieri, A. & Ficetola, G. F. Optimal sequence similarity thresholds for clustering of molecular operational taxonomic units in DNA metabarcoding studies. Molecular Ecology Resources 00, 1–14 (2022).
Ficetola, G. F. et al. An In silico approach for the evaluation of DNA barcodes. BMC Genomics 11, 434 (2010).
Zinger, L. et al. DNA metabarcoding—Need for robust experimental designs to draw sound ecological conclusions. Molecular Ecology 28, 1857–1862 (2019).
Dickie, I. A. et al. Towards robust and repeatable sampling methods in eDNA based studies. Molecular Ecology Resources (2018) doi:10.1111/1755-0998.12907.
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 in Ecology and Evolution 6, 232–236 (2015).
Baselga, A. Partitioning the turnover and nestedness components of beta diversity: Partitioning beta diversity. Global Ecology and Biogeography 19, 134–143 (2010).
Legendre, P. & De Cáceres, M. Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecology Letters 16, 951–963 (2013).
Smithson, M. & Verkuilen, J. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychological methods 11, 54–71 (2006).
Bürkner, P.-C. brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software 80, 1–28 (2017).
Muggeo, V. M. R. segmented: An R Package to Fit Regression Models with Broken-Line Relationships. 8, 7 (2008).
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).

There is NO Competing Interest.

Download PDF

Journal Publication

published 17 Jan, 2024

Read the published version in Nature Plants →

Version 1

posted

You are reading this latest preprint version

The importance of species addition versus replacement varies over succession in plant communities after glacial retreat

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

eDNA sampling

Traditional sampling

β-diversity measures

Statistical analyses

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1