Succession comprises a sequence of threshold-induced community assembly processes towards multidiversity

Research on successions and community assembly both address the same processes such as dispersal, species sorting, and biotic interactions but lack unifying concepts. Recent theoretical advances integrated both research lines proposing a sequence of stochastic and deterministic processes along successional gradients. Shifts in ecosystem states along successional gradients are predicted to occur abruptly once abiotic and biotic factors dominate over dispersal as main driver. Considering the multidiversity composed of five organismal groups including plants, animals, and microbes, our results imply that stochastic, likely dispersal-dominated, processes are replaced by rather deterministic processes such as environmental filtering and biotic interactions after around 60 years of succession in a glacier forefield. The niche-based character of later successional processes is further supported by a decline in multi-beta-diversity. Our results may update concepts of community assembly by considering multiple taxa, help to bridge the gap between research on successions and community assembly, and provide insights into the emergence of multidiverse and complex ecosystems.

S uccession 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 ideas 1 . 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 succession 1 . 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 assembly 2,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 shift 1,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 taxa 5,6 . Multidiversity, an aggregate measure of biodiversity that integrates the standardized diversities of multiple taxa 7 , has been shown to reflect ecosystem functionality 8 , species composition 9 , and multitrophic interactions 10 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 drivers [11][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 communities 16,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 only 16 .
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 diversity 1,17 . The framework by Chang and HilleRisLambers (2016) 1 predicts an increase in functional and phylogenetic diversity of local singletaxon 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 patterns 18 , 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 age 19 and regulate plant-animal interactions 20 , 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 succession 21,22 , and studies on successional convergence have mostly been focused on singular taxa neglecting potential interactions and co-occurrences across organismal groups [23][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 multibetadiversity 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 decline 7,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 Alps 26 . 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 latesuccessional communities.

Results
Breaking point in multidiversity. Multidiversity follows a nonlinear 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 inference 28 . 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 (nonsignificant relationship with time) for the late successional stage (early: t 42 = 5.77, p < 0.001, r 2 = 0.44; late: t 64 = −0.55, p = 0.58, r 2 = 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).
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 growthseason 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 environments 29,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 1-4). 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 events 31,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 interactions 33,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).
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; mbD early = 0.75 ± 0.09; mbD late = 0.70 ± 0.08; Fig. 4).

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   Single taxa occurrence along the successional gradient for each taxon that occurs on at least three plots along the gradient. 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. ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03372-2 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 change 4,6 or for restoration efforts of anthropogenically altered landscapes 2,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 America 11,[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 community 33,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 negative 40 . For instance, in host-pathogen relationships, more diverse host communities harbor a higher diversity of obligate pathogens 41 and below-ground fungal diversity is promoted by diverse bacterial assemblages through complex feedback loops 42 . 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 sources 43 . 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 environments 44,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 environment 46 . Some fungal-feeding taxa show preferences for certain types of fungal hyphae that frequently have distinct associations with selected plant species as pathogens 47 , mycorrhizal partners 48 , or litter decomposers 49 , 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 retreat 51 , 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 mature 52,53 . Although co-occurrence patterns must not necessarily reflect biotic interactions 54 , 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 assembly 55,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-webs 1,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 reoccurring 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. Soilmicrobiome 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 N 2 ) 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 groups 7 . 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 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 dynamics 1 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 regimeshifts of the ecosystem and can be interpreted as ecological thresholds 27 . 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 thresholds 58 . 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.9 28 . 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 datasets [59][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 CSSnormalized 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 levels 29 . 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 crosschecking 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 community 33,34 . All path models were estimated using the R-package "lavaan" v0.6-7 62 .
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 distances 63,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-8 65 .
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.

Data availability
Raw sequences of next-generation 16  Code availability