Introduction

Succession and community assembly research are neighboring fields, each with a long history of constructive debates and a large body of ecological literature. Both research lines are founded on overlapping concepts considering processes such as dispersal, environmental filtering, biotic interactions, and stochasticity as structuring elements of local communities and ecosystems. Nevertheless, succession and community assembly research focus on different spatial and temporal scales, which may have hindered a mutual exchange of ideas1. Succession research relates to the initial development of ecosystems and communities over time. Community assembly studies, on the other hand, try to elucidate the past drivers of local and recent diversity patterns regardless of temporal components. Attempts to synthesize the two fields led to the development of an integrated conceptual framework of succession and community assembly dynamics. Here, drivers of community structure are considered from an explicitly temporal perspective, assuming a temporal sequence of assembly processes shifting from stochastic (likely dispersal) to rather deterministic (niche-based and interaction-mediated) processes along the course of succession1. The proposed framework also emphasizes the importance of threshold dynamics, predicting that gradual changes in the abiotic and biotic environment are followed by rapid shifts in the predominant processes that determine community assembly2,3. These shifts in processes are expected to be associated with various community-level responses such as changes in diversity, species composition, and functionality which ultimately leads to a community state shift1,2,4. Recent research suggests that sudden ecosystem shifts are not restricted to changes in species composition of single organismal groups but that local diversity patterns emerge from on-site interactions between multiple taxa5,6. Multidiversity, an aggregate measure of biodiversity that integrates the standardized diversities of multiple taxa7, has been shown to reflect ecosystem functionality8, species composition9, and multitrophic interactions10 more accurately than diversity measures considering only one taxon. Additionally, the responses of different organismal groups to successional gradients may differ either due to stochastic or deterministic drivers11,12,13,14,15. For instance, the assembly of belowground bacterial and fungal communities follows individual trajectories during primary succession that differ from the assembly of vascular plant communities16,17. Thus, the concept of multidiversity is a well-suited approach that can reflect various community-level responses and allows more concise conclusions about different ecosystem states along successions, as well as underlying assembly processes than considering a single or few taxa only16.

The assumed changes of assembly processes along successions, a decrease in the importance of stochasticity, and a simultaneous increase in niche-based and interaction-based processes, have consequences for the taxonomic, functional, and phylogenetic composition of communities and thus beta diversity1,17. The framework by Chang and HilleRisLambers (2016)1 predicts an increase in functional and phylogenetic diversity of local single-taxon communities in later successions due to niche differentiation and intensified biotic interactions. In a multidiversity context where multiple taxa are considered, increased interactions and thus potentially strong dependencies between taxa may result in aggregated co-occurrence patterns18, which may have consequences for the beta-diversity between local communities—and not necessarily for the functional and phylogenetic diversity within a community. For instance, interactions between plants and soil-inhabiting microorganisms intensify with successional age19 and regulate plant-animal interactions20, which is reflected in the co-occurrence of the partners involved in these interactions. Such interactions may thus result in reduced species turnover across local assemblages in older successional stages. Yet, it is an ongoing debate whether ecological communities generally converge towards a core set of species over the course of succession21,22, and studies on successional convergence have mostly been focused on singular taxa neglecting potential interactions and co-occurrences across organismal groups23,24,25. To specifically test the extent of species turnover across several taxa, we introduce multi-betadiversity, an aggregated measure of taxonomic community dissimilarity. This index allows to test whether the expected increase in biotic interactions is reflected in a reduced compositional variation (i.e., reduced multi-betadiversity regarding several taxa) of the multidiverse community, i.e., a community that comprises a number of taxonomic groups, once niche-based and interaction-mediated assembly processes dominate over stochastic dispersal as main mechanisms shaping multidiversity.

Assuming that multidiversity more precisely reflects ecosystem responses and functions, we tested predictions deduced from contemporary succession hypotheses using a dataset on multidiversity comprising inventories of plants, animals, and microorganisms along a successional gradient. Most studies on multidiversity have been undertaken in well-developed ecosystems that have undergone severe anthropogenic alterations in the past and focused on the drivers of biodiversity decline7,9,10. The mechanisms behind the successional emergence of multidiversity and ecosystem complexity under natural conditions, however, are yet poorly understood. Thus, considering multidiversity within a successional framework will help to gain insights into the temporal dynamics of the processes that shape the initial emergence of biodiversity in natural ecosystems that comprise multiple organismal groups with unique and complementary ecological functions. We adopted the concept of multidiversity to empirically test predictions of the integrated framework of succession and community assembly dynamics on an ecological gradient of primary succession following glacial retreat in the Austrian Alps26. We assessed the multidiversity and multi-betadiversity of five organismal groups, including vascular plants, bryophytes, invertebrates, and soil-inhabiting fungi and bacteria on 110 plots spanning 170 years of primary succession. We estimated multidiversity by the mean of ranked normalized Shannon-diversities and calculated multi-betadiversity as the mean Bray-Curtis dissimilarities between plots regarding the composition of individual taxa. According to the conceptual framework by Chang and HilleRisLambers (2016)1, we predicted that the multidiverse community will undergo a rapid shift that is induced by a sudden change of the biotic and abiotic environment along the successional gradient. We expect this shift to be associated with increased biotic interactions and stronger environmental influences that result in a reduced compositional variation between plots of the multidiverse community.

Accordingly, we identified a threshold in multidiversity development (i.e., a shift from an increase of multidiversity with time to a rather stationary phase) and estimated the relative importance of stochasticity, as well as environmental and biotic drivers on the multidiverse community before and after the threshold using path analysis. We further highlighted the more structured character of community assembly after the shift by comparing multi-betadiversity estimates of the early and late-successional communities.

Results

Breaking point in multidiversity

Multidiversity follows a non-linear trajectory along the successional gradient. We screened the gradient for the predicted shift from a developing (characterized by an increase in multidiversity over time) to a more stationary ecosystem (characterized by non-monotonic variation in multidiversity over time) using breaking point analysis, which detects changes in the slope of associations. According to Groffman (2006)27 breaking points represent ecological thresholds, which is characterized by a rapid change of the slope in the association between time since deglaciation as explanatory variable and multidiversity as dependent variable in our study. We found such a breaking point by cross-validating a broken-stick model to alternative models that either did not include a breaking point or included a disjunct breaking point using Bayesian inference28. We then used Bayes Factors analysis to validate the exact location of the breaking point after about 60 years of succession, i.e., plot index 44 along the gradient (lower-upper bounds: 40–74 years after deglaciation; see “Methods”, Fig. 1, Supplementary Table 1). Resulting piecewise linear models showed a significant linear increase in multidiversity during the early successional stage and a stationary multidiversity (non-significant relationship with time) for the late successional stage (early: t42 = 5.77, p < 0.001, r2 = 0.44; late: t64 = −0.55, p = 0.58, r2 = 0.004; Fig. 1). To assess whether the delineated threshold is also reflected in the composition of the multidiverse community, we calculated the relative abundance of each taxon along the successional gradient and visualized it in a temporally ordered multi-taxa community table of species occurrences (see “Methods”, Fig. 2).

Fig. 1: Multidiversity along the successional gradient. Multidiversity increases during the early successional stage and remains stationary as succession continues.
figure 1

The vertical dashed line indicates the estimated threshold in multidiversity development, which is the threshold in community assembly as predicted by the framework of Chang and HilleRisLambers (2016)1. The grey regression line indicates the linear correlation between multidiversity and time since deglaciation using a broken-stick model. Blue dashed lines indicate the 2.5% and 97.5% quantiles of fitted values. The yellow curve on the x-axis resembles the posterior distribution for the estimate of the breaking point, i.e., higher density of the frequency indicates a higher probability for the breaking point to be located at a given value.

Fig. 2: Single taxa occurrence along the successional gradient for each taxon that occurs on at least three plots along the gradient.
figure 2

Tiles indicate the relative abundance of each taxon at a given plot relative to the mean abundance of the taxon along the gradient. Taxa are ordered by the weighted average plot of their occurrence, i.e., taxa that reach their abundance optima early in succession appear in the top left and late successional specialist in the bottom right. Relative taxon abundance is color-coded and reflects z-scores with the mean abundance as 0 and one standard deviation as +1 or −1, rescaled between −1 and 1. Multidiversity is indicated as light pink density plot on top of the community composition plot. Grey density plot on the left side indicates the range size of the taxa, i.e., whether their occurrence is restricted to a narrow range of the successional gradient or whether taxa occur along the whole gradient. Black bars on the right represent the organismal group of the taxa in each row, symbols of organismal groups are ordered by decreasing frequency; white dashed horizontal line marks the threshold.

Drivers of multidiversity

To quantify the putative relative importance of stochastic dispersal, species interactions, and environmental drivers before and after the threshold, we specified a path model in which time since deglaciation and mean growth-season soil temperature act as exogenous variables and directly affect diversities of single organismal groups. Soil temperature is an environmental factor that is important for the diversity of several trophic levels in alpine environments29,30. We also tested the effects of other environmental variables, such as soil nutrient content and soil-pH value but could not detect any significant effects and thus removed these variables from the final model (see “Methods” and Supplementary Data 14). After accounting for environmental variation, we consider the effect of time since deglaciation as a proxy for stochastic events, such as successful dispersal and colonization events31,32, leading to cumulating multidiversity over time. The total effects of time since deglaciation (proxy for stochastic dispersal) and environmental drivers on multidiversity were estimated as indirect effects mediated through the organismal groups. Biotic interactions were modeled as the residual covariance among the diversity values among all organismal groups that remained after accounting for stochasticity and environmental variation. If corrected for such underlying common causes, strong positive covariances within communities can be interpreted as biotic interactions33,34. We further tested whether the residuals of the estimated models show significant spatial autocorrelation but could not detect any unexplained variance that could be attributed to geographic variation (see Methods and Supplementary Table 2).

The path analysis provided good model fit (Model fit: pchi-square > 0.05; CFI > 0.95; TLI ≥ 0.9; SRMR < 0.09; RMSEAearly < 0.05; RMSEAlate = 0.066 with lower bound of confidence interval = 0) and revealed different influences of the exogenous variables between the early and late successional environments (Fig. 3, Supplementary Data 4). During early succession, time since deglaciation had a strong positive effect on the diversity of bryophytes (β = 0.66, p < 0.01), arthropods (β = 0.62, p < 0.01), and bacteria (β = 0.32, p = 0.04), whereas there was no significant effect of temperature on any organismal group. Time since deglaciation also had a significant positive indirect effect on total multidiversity (β = 0.74, p < 0.01) during early succession, whereas temperature did not affect multidiversity (β = 0.01, p = 0.93, Fig. 3a). Furthermore, diversities of the five organismal groups varied independently of each other suggesting little or no mutual influences. In late succession, mean growth-season soil temperature had a significant positive direct effect on the diversity of bryophytes (β = 0.35, p = 0.01) and a clear positive indirect effect on multidiversity mediated through all organismal groups (β = 0.30, p = 0.01). Plot age had a negative effect on fungal diversity (β = −0.35, p < 0.01) and a positive effect on bryophyte diversity (β = 0.30, p = 0.01) but no effect on multidiversity (β = −0.01, p = 0.89). We further detected positive residual covariances between vascular plant (β = 0.26, p = 0.04), invertebrate (β = 0.25, p = 0.05) and bacterial (β = 0.29, p = 0.03) diversities with fungal diversity, suggesting mutual influences between taxa (Fig. 3b).

Fig. 3: Path analysis regarding the effects of stochastic and deterministic drivers on the multidiverse community for the early and late successional stages.
figure 3

Time since deglaciation and mean growth-season soil temperature act as exogenous variables with directed effects on the diversities of vascular plants, bryophytes, invertebrates, and soil-inhabiting fungi and bacteria. Indirect effects of exogenous variables on multidiversity were calculated by mediation analysis through the direct effects of the organismal groups on multidiversity. Biotic interactions are modeled as covariances between organismal groups. Numbers represent standardized path coefficients and are given for significant paths only. The path model was calculated for (a) the early and (b) the late successional stage separately. Silhouette images were obtained from PhyloPic (http://phylopic.org) under a public domain license (CC0 1.0 license).

Multi-betadiversity

We predicted that the shift from dispersal dominated stochastic events in early successional stages to more deterministic niche-based processes that are induced through environmental filtering and biotic interactions in late successional stages is associated with lower beta-diversity among communities after the threshold. To test this prediction, we calculated a measure of total community dissimilarity, multi-betadiversity (mbD) for both successional stages by the mean Bray-Curtis dissimilarities of the multidiverse community (see “Methods”). Of n = 5 organismal groups surveyed, n = 4 showed a decrease in betadiversity and in total, multi-betadiversity was significantly lower in the late successional stage (mean ± SD; mbDearly = 0.75 ± 0.09; mbDlate = 0.70 ± 0.08; Fig. 4).

Fig. 4: Betadiversity values for each organismal group and community-wide multi-betadiversity values for the early and late successional stages.
figure 4

Mean Bray-Curtis dissimilarities were used to calculate betadiversity and multi-betadiversity values. Green violin plots resemble the dissimilarity values of the early successional stage (n = 44 plots), rose violin plots indicate the late successional stage (n = 66 plots). Boxplots indicate the upper and lower quartile of the estimates, and the size of the notch represents the confidence interval of the median. a mean of all organismal groups, i.e., multi-betadiversity (b) bacteria, c fungi, d vascular plants, e bryopyhtes, and (f) invertebrates. Silhouette images were obtained from PhyloPic (http://phylopic.org) under a public domain license (CC0 1.0 license).

Discussion

Our study is an empirical test of the conceptual framework integrating succession and community assembly dynamics proposed by Chang and HilleRisLambers (2016)1 and modified here in the context of multidiverse communities. We showed that threshold dynamics play an important role in the generation of multidiversity under natural conditions and that succession comprises a sequence of different ecosystem states that can be detected using multidiversity. These threshold-mediated shifts in ecosystem states are associated with substantial changes in the community assembly processes. So far, the importance of threshold dynamics has been recognized for anthropogenically altered ecosystems either as catastrophic shifts with respect to climate change4,6 or for restoration efforts of anthropogenically altered landscapes2,3. Here, we revealed a threshold after about 60 years of ecosystem development when multidiversity became stationary. Significant thresholds that are marked by a steep increase followed by stationarity after 40–60 years of succession appear to be a generalizable pattern that occurs in various aspects of ecosystem development, such as plant, invertebrate, and microbial diversities and functionality in glacier forefields in Europe and Northern America11,35,36,37,38. The threshold became also evident in the change of community composition with different sets of taxa before and after the breaking point. Prior to the threshold, communities mainly consist of species that reach their abundance optima early in succession. These pioneering species were soon accompanied by taxa with no clear preference of early or late successional stages leading to an increase in multidiversity over time. After the threshold, specialists for early successions are replaced by specialists for late successional stages, consequently multidiversity remained stationary over time (Fig. 2). Thus, acknowledging the effects of threshold dynamics during the development of natural multidiverse ecosystems, and their universality in ecological systems provides valuable insights into the patterns and processes of initial ecosystem development.

Our results indicate an ecosystem state shift associated with different assembly processes: the first 60 years of succession are characterized by stochastic species additions, afterwards rather deterministic processes such as niche-filtering and biotic interactions dominate. Although different organismal groups follow individual trajectories during primary succession, multidiversity is mainly promoted by stochastic drivers during the initial phase of ecosystem development, most assumably by heterogenous dispersal events. Over the course of succession, stochasticity is replaced by environmental filtering and biotic interactions as the structuring mechanism of multidiversity.

Our data do not allow to directly test for mutual influences between organismal groups on the taxon level, but we find a strong pattern of positive covariances of the diversities of various taxonomic groups after accounting for environmental variation that is indicative of biotic dependencies within the community. Community covariances have initially been applied to evaluate the relative importance of biotic interactions in multitrophic communities by estimating synchronized fluctuations in population abundances where positive covariances have been interpreted as a result of facilitative interactions within a community33,34,39. On the diversity level as analyzed in this study, however, positive covariances among organismal groups do not necessarily reflect only facilitative interactions, because diversity in one group is likely to be positively linked to diversity of another trophic group through specialized interactions, either positive or negative40. For instance, in host-pathogen relationships, more diverse host communities harbor a higher diversity of obligate pathogens41 and below-ground fungal diversity is promoted by diverse bacterial assemblages through complex feedback loops42. A previous study in the same glacier forefield showed that after more than 50 years of ecosystem development, the soil microbial community is mainly supported by carbon from recent plant production, whereas during the initial stage of ecosystem development, the microbial community is mainly sustained by ancient carbon sources43. The time point of the shift in the main carbon source largely coincides with our delineated threshold in multidiversity development. The increased interdependence of the soil microbial and plant community during the later successional stage is also reflected in the positive covariances of vascular plant, fungal and bacterial diversities in our path model. Previous studies have also shown strong mutual dependencies of invertebrate and fungal communities in developing alpine environments44,45. Invertebrates increase the diversity of soil-inhabiting microbial taxa by propagule dispersal and create a variability of suitable habitats through modifications of the physical environment46. Some fungal-feeding taxa show preferences for certain types of fungal hyphae that frequently have distinct associations with selected plant species as pathogens47, mycorrhizal partners48, or litter decomposers49, which well explains the pattern of covariance among these groups in our analysis. Further, the importance of fungi-animal interactions for the development of ecological complexity was predicted in foundational works on succession theory by Connell and Slatyer (1977)50, as heterotrophic fungi were expected to feed on the carcasses and dung of animals and in turn, certain arthropods were hypothesized to be reliant on decomposed substrate for the developmental stages of their life cycles. Although bryophytes have been shown to closely interact with epiphytic microorganisms after glacial retreat51, potential interactions of bryophytes with microorganisms may not be detectable in our path model as we sampled soil microbial communities and not all substrate types colonized by bryophytes in our study site (e.g., rocks, scree, and litter).

The threshold-mediated character of successional processes is further supported by a pronounced decline in multi-betadiversity that is indicative of an increased dependence structure among the organismal groups. In line with that assumption, we detected a lower beta-diversity in all those organismal groups that also show positive community covariances in the path model. Out of the five organismal groups we studied, only bryophytes were characterized by a higher betadiversity in the late successional stage that might be attributable to increasing variability in microhabitats as the communities mature52,53. Although co-occurrence patterns must not necessarily reflect biotic interactions54, reduced beta-diversity between plots in later successional stages still may be indicative for stronger biotic interactions, which finds support in studies on plant community assembly55,56. The literature suggests the increasing importance of interactions in community development over time. One explanation put forward is that with increasing niche differentiation, energy and nutrient flows in an ecosystem follow more complex biochemical pathways that necessitate biotic interactions to retain nutrients in the ecosystem through closed mineral cycles and complex food-webs1,17. Our results support the increasing importance of biotic interactions as in the late successional stage, our path model revealed positive community covariances between the soil microbial, vascular plant, and arthropod communities that may represent such complex nutritional cycles comprising primary producers, heterotrophic organisms, and decomposers For the organismal groups forming these putative cycles, we find a reduced betadiversity between the plots that indicates a strong interrelatedness of individual taxa after the threshold, which may lead to re-occurring core assemblages of species of several taxa. These assumptions and the threshold suggested by our data after 60 years are further supported by the study of Bardgett et al. (2007)43 that revealed a high dependence of heterotrophic microorganisms on plant communities in the Ödenwinkel forefield after a similar number of years. We thus recommend extending the integrative framework of succession and community assembly for multiple interacting taxa that mutually shape their diversity and composition and thus cause a reduction of betadiversity.

Understanding the relative importance and temporal dynamics of deterministic and stochastic processes is a key challenge in community ecology, especially in natural systems and has been in the center of broad debates among ecologists. Using the multidiversity and multi-betadiversity approach allows us to comprehensively understand the processes leading to stable, resilient, and complex ecosystems, which may remain vague in single-taxon approaches. Thus, our study contributes to a synthesis of community ecological theories into succession research, acknowledging the fundamental importance of abrupt state shifts in natural ecosystems. These results are not only a proof of concept but also (re-)emphasize that succession is a multi-faceted rather than a linear process that comprises a sequence of assembly processes towards multidiversity and ecosystem complexity.

Methods

Study design

The study was conducted in the long-term ecological research platform Ödenwinkel which was established in 2019 in the Hohe Tauern National Park, Austria (Dynamic Ecological Information Management System—site and dataset registry: https://deims.org/activity/fefd07db-2f16-46eb-8883-f10fbc9d13a3, last access: March 2021). A total of n = 135 permanent plots was established within the glacier forefield of the Ödenwinkelkees which was covered by ice at the latest glacial maximum in the Little Ice Age around 1850. For this study, we used a subset of n = 110 plots with complete datasets with all biotic and abiotic variables available. The plots represent a successional gradient spanning over 1.7 km in length and were distributed within the glacier forefield to reflect the gradient of glacial retreat at an altitude ranging from 2070 to 2170 m a.s.l (distance between neighboring plots, median ± SD: 19 m ± 6.3 m). Each plot was given a unique index number according to the location along the gradient: the plot with index 1 is located closest to the glacier, plot index 135 is the farthest away from the glacier. Plot age was estimated based on data of historical glacial extent and as glacial loss did not occur at a constant rate, plot locations were pre-selected to represent a linear gradient of time since deglaciation rather than a spatial gradient with equally spaced intervals between plots. To ensure independence of plots in areas of slow glacial retreat, we chose a minimum distance of 5 m between two neighboring plots. Plots were defined as squares with an area of 1 m² and were all oriented in the same cardinal direction. Further details on the design of the research platform, exact plot positions and details on the surrounding environment, as well as on the historical glacial extent can be found in Junker et al. (2020)26.

Sampling

In 2019, we identified all vascular plant and bryophyte species present on the plots and estimated their cover with a resolution of 0.1%. We sampled above-ground arthropod diversity by installing two pitfall traps on each plot. Traps were set active for a total of n = 7 days. The abundance of all arthropods, excluding Collembola and Acari, larger than 3 mm was counted. The abundance of Collembola and Acari and of animals smaller than 3 mm was estimated based on random samples of aliquots of the total sample. All arthropods and other animals are identified to the order level. Soil-inhabiting bacteria and fungi were sampled from soil cores from an approximate depth of 3 cm, as soil development in the proglacial study area was limited to the top layers of the pedosphere. Soil-microbiome samples were analyzed by next-generation sequencing, and microbiome profiling of isolated DNA was performed by Eurofins Genomics (Ebersberg, Germany). Prior to the statistical analysis of microbial communities, we performed a cumulative sum scaling (CSS) normalization (R-package “metagenomeSeq” v1.28.2)57 on the count data to account for differences in sequencing depth among samples. Detailed information on the sampling strategies of all organismal groups can be found in Junker et al. (2020)26. Soil temperature measurements were done by installing temperature loggers (MF1921G iButton, Fuchs Elektronik, Weinheim, Germany) 10 cm north of each plot center, at the same depth of 3 cm at which the microbial samples were taken. Mean growth-season soil temperature was calculated based on the recordings ranging from 26th of June to 16th of September representing the period in which the plots were free of permanent snow cover before and after the winter 2019/2020. In 2020, soil samples were taken and soil nutrients (Ca, P, K, Mg, and total N2) as well as soil pH were measured on all plots by AGROLAB Agrar und Umwelt GmbH (Sarstedt, Germany). Soil nutrient analysis was performed according to the Ö-Norm Boden: L 1087: 2012-12 (K and P—mg/1000 g), L 1093: 2010-12 (Mg—mg/1000 g), and L 1083: 2006-04 (pH). Total N2 (%) was determined according to the DIN EN 16168: 2012-11.

Statistics and reproducibility

Calculation of multidiversity

Multidiversity is defined as the cumulative diversity of a number of taxonomic groups7. The multidiversity of the n = 110 plots composed of the diversities of vascular plants, bryophytes, invertebrates, fungi, and bacteria present in each plot and was defined as follows: First, we calculated the Shannon diversities of each of the taxonomic groups in each plot. Second, we ranked the plots by increasing diversity for each of the five taxonomic groups individually, i.e., for each plot we received n = 5 ranks. The mean rank of each plot was defined as multidiversity mD. For better interpretability of the index, we then normalized the mean ranks as

$$x{{{\prime} }}=(x-{\min }(x))/({\max }(x)-{\min }(x))$$
(1)

to scale them between zero and one. Ranks were used to give the same weighting to each taxonomic group despite deviations in the absolute values of Shannon diversity and to reduce the impact of outliers in mD values. To account for differences in the sequencing depth of the microbial raw dataset and sampling effort of the invertebrate trapping, we performed multiple rarefactioning prior to the calculation of Shannon diversity by averaging the results of n = 999 iterations (R-package “rtk” v0.2.5.7) and used original read numbers instead of using the CSS-normalized dataset for the microbial dataset.

Breaking point analysis

According to the hypothesis of community assembly dynamics1 and by visual inspection of the relationship between multidiversity and successional age of the plots, we expected two different stages of community establishment along the successional gradient. These stages are separatable by the transition from a developing (characterized by an increase in multidiversity over time) to a more stationary ecosystem (characterized by non-monotonic variation in multidiversity over time). Such non-linear associations are suggestive of regime-shifts of the ecosystem and can be interpreted as ecological thresholds27. Piecewise regression models have been shown to be the most suitable tool for the detection of ecological thresholds in natural systems, as they are able to correctly estimate the probability, as well as the number and position of ecological thresholds58. We screened the successional gradient for the existence of a threshold by comparing four different models with mD as the dependent and plot age as the independent variable using the R-package “mcp” v0.3.0.928. The mcp-method fits piecewise regression models with a pre-defined number of breaking points and is based on Bayesian inference. We specified a base model m1 (mD-values remain constant along the gradient), and compared it to three alternative models (m2 = constant linear increase of mD-values without a breaking point, m3 = one breaking point in mD-values with a segregated slope (abrupt-threshold model), m4 = one breaking point in mD-values with a joined slope (broken-stick or smooth-threshold model)). For each model, we separately ran three Markov Chain Monte Carlo estimators with a uniform prior for a total of n = 11,000 generations while discarding the first n = 1000 generations as burn-in. Model convergence was estimated by visually inspecting the trace plots and checking that all model parameters reached stationarity. We then compared the predictive performance of the four models using leave-one-out-cross-validation and confirmed the exact position of the breaking point using Bayes Factors. Both validation methods can be applied for a robust and accurate testing of competing hypotheses in ecological datasets59,60,61.

Community taxa occurrence

To visualize the distribution of individual taxa and thus the community-wide taxonomic turnover along the successional gradient, we estimated the abundance optimum of each taxon that occurred on at least three plots. Abundance optima were estimated by calculating the weighted mean plot of occurrence with abundance of each taxon as weighting factor (i.e., cover for vascular plants and bryophytes, individual count for invertebrates, and CSS-normalized read number for microorganisms). We further calculated the relative abundance of each taxon per plot by using z-scores of abundances. The mean abundance was set as 0 with one standard deviation as +1 or −1. To allow a comparison across taxa, z-scores were then rescaled between −1 and 1. Range size was estimated by calculating the variance of occurrence plots (i.e., the span of plots on which a taxon was found along the gradient) weighted by the abundance of the taxon on the respective plots.

Path analysis

We used path analysis (i) to model the influence of the abiotic environment on the diversities of all taxonomic groups individually, (ii) to estimate the strength of covariance between the diversities of those groups, (iii) to infer the effect sizes of the diversity of each group on total mD, and (iv) estimate the strength of indirect effects of the exogenous variables that are mediated through the organismal groups on mD. We built separate models with identical structure for the early (n = 44 plots) and late (n = 66 plots) successional stages. The stages were delineated by the threshold identified in the breakpoint analysis. Exogenous variables in the model were time since deglaciation and mean growth-season soil temperature. Time since deglaciation reflects the plot age and can be seen as a proxy for the increasing chance of heterogeneous dispersal events that occur over time. The mean temperature of the growing season is an estimate for environmental heterogeneity and has been shown to affect diversity directly and indirectly on various trophic levels29. All variables were scaled by subtracting the mean and dividing through the standard deviation. We first ran a full model that included soil pH and soil nutrient content as additional exogenous variables. None of the two variables showed a significant direct effect on any organismal group or a significant indirect effect on multidiversity during either the young or late successional stage. Accordingly, we stepwise removed the two variables from the model while cross-checking whether the effect strength of one variable became significant in the absence of the other variable. As no significant effects occurred or the model fit significantly decreased, we decided to remove both variables from the final model. Within the final model, the strength of residual covariance between the diversities of all taxonomic groups was estimated while accounting for influences of time since deglaciation and temperature. If corrected for an underlying common cause, such as environmental autocorrelation, strong covariances between members of a community can be interpreted as biotic interactions and especially positive community covariances are indicative of facilitative effects within the community33,34. All path models were estimated using the R-package “lavaan” v0.6-762.

Test for spatial autocorrelation

Studies taking advantage of space-for-time substitution generally are prone to erroneous conclusions caused by spatial autocorrelation, especially when sampling points are located at varying distances63,64. We tested for the presence of residual spatial autocorrelation in the estimated path models by calculating spatial neighbor matrices for distance classes of 5 m and 10 m between plots (i.e., the potential to include a maximum of three or five neighboring plots). For both distance classes, we estimated Moran´s I for the residuals of the piece wise linear models of the breaking point analysis and the case-wise residuals of the exogenous variables of the path models using the R-Package “spdep” v1.1-865.

Multi-betadiversity

We defined Multi-betadiversity (mbD) as the cumulative averaged pairwise-dissimilarity across a number of organismal groups. First, we split the total community composition data table (sites x species, n = 110 plots) for each organismal group into two tables containing the plots before (n = 44 plots) and after (n = 66 plots) the threshold. Second, we calculated pairwise Bray-Curtis dissimilarities for each of the organismal groups individually. Third, we calculated mbD as the mean of the n = 5 dissimilarity estimates of each plot pair.

Reporting summary

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