Introduction

Plant functional traits are morphological, physiological or life-history characteristics that directly or indirectly influence plant fitness via their effects on survival, growth and reproduction (Violle et al. 2012). Since Hutchison (1957), the multidimensional space spanned by these traits served to define the ecological niche of a species and the degree of niche division (the functional diversity, FD) amongst the species in animal and plant communities (Violle and Jiang 2009). FD might serve as an indicator of functional variability during the course of community assembly (e.g. Lavorel and Garnier 2002; Díaz et al. 2007; Garnier et al. 2016).

Three basic processes govern plant community assembly: habitat filtering (Kraft et al. 2015), species interactions (Aschehoug et al. 2016) and neutral colonization and extinction dynamics (Hubbell 2001). The first two processes select for certain species traits and should result in a non-random community composition compared to neutral assembly (Ulrich et al. 2016). In this context, Kahmen and Poschlod (2004) and Ulrich et al. (2014) reported a covariance with species richness and abundance during plant succession for several important traits, especially for plant height and seed mass. However, most of the traits analysed by Kahmen and Poschlod (2004) did not show predictable successional trends, indicating trait assembly to be independent of the specific processes of community assembly and to be in accordance with neutral processes. Similarly, Schleicher et al. (2011) failed to confirm that community assembly during plant succession can be predicted from trait distributions (being trait driven in the sense of Keddy 1992). It is possible that filtering and competition processes counterbalance one another, making it difficult to trace the drivers of FD (Cornwell et al. 2006; Kumordzi et al. 2015). For example, increased interspecific competition and strong habitat filter regimes during succession might lead to limited FD (MacArthur and Levins 1967) even in cases of overall increasing richness. Such a trade-off between richness and interaction effects might be manifest in stable functional diversity during succession (Mason et al. 2008).

In a long-term study of primary plant succession, Zaplata et al. (2013) found that total plant cover tended to increase over time, whilst at times being associated with high small-scale spatial variability in cover and in species composition. Such a spatial variability might mask the respective successional trends in trait expression and trait diversity. Additionally, environmental variability, particularly soil variability, is linked to the respective variability in community composition, making it difficult to assess whether community composition or soil variability triggers changes in FD (Robertson et al. 1988).

To disentangle the factors that influence FD during plant succession we would need to partition the major components that influence trait expression, particularly richness and abundance effects and their environmental correlates. One tool for such a task is the Price equation that had been primarily used in ecology to infer the influences of species richness and community composition on ecosystem functioning (Loreau and Hector 2001; Fox and Kerr 2012; Isbell et al. 2018). The Price approach in ecology applies the multiplication law of discrete calculus to partition a multiplicative term into two components, each containing the changes in the two factors involved (Fox 2006; Bourat et al. 2023). Applied to FD, let SA and SB denote the species richness in two communities A and B and \(\overline{FD }\) the mean expression of functional diversity. The change in FD (∆FD = ∆\((\overline{FD }S\)) can be partitioned by

$$\Delta FD=\overline{{FD }_{B}}{S}_{B}-\overline{{FD }_{A}}{S}_{A}=\overline{{FD }_{A}}\Delta S+{S}_{A}\Delta \overline{FD }+\Delta \overline{FD }\Delta S=\overline{{FD }_{A}}\Delta S+{S}_{B}\Delta \overline{FD }$$
(1)

Fox (2006) interpreted the three partition form in terms of richness, habitat context (including random effects) and community composition (covariance), respectively. Importantly, the partitions of Eq. 1 itself contain terms to which the multiplication rule can be applied. It can be extended into a theoretically unlimited number of further partitions that cover additional aspects of the change in community characteristics. For instance, Ulrich et al. (2022) proposed a five partition version that explicitly contains terms for changes in richness, average trait expression, relative and total species abundances and covariation in traits and relative abundance. Such a partitioning offers an elegant way to disentangle different drivers of plant primary succession within a trait-based approach.

Here, we use data from long-term monitoring of initial primary plant succession carried out between 2005 and 2011 in a German artificial catchment created to study the initial hydrological processes and the development of soil and vegetation (Frouz et al. 2020). Our major goal was to study the relative importance of major drivers behind the change in functional diversity. In earlier work we reported high small-scale variability in the pattern of species occurrences and phylogenetic community composition (Ulrich et al. 2017). We also found a steady increase in the diversity of important functional traits (Ulrich et al. 2014). FD and multifunctionality, defined as the manifold of ecological ecosystem services realized in a community (Hector and Bagchi 2007; Gamfeldt and Roger 2017), were comparatively low, although increasing in time (Winter et al. 2019). Based on these studies and the above discussion we ask three basic questions. (1) How do three major factors (increase in richness, change in abundance, and the dynamics in community assembly) influence FD during early plant succession? (2) Are there trade-offs between these factors that stabilize the dynamics in functional diversity in space and time? (3) Does the importance of each factor predictably change in time and across local patches during plant succession?

Methods

From 2005 to 2011, we studied the early vegetation succession of a six-ha area in the artificial catchment Chicken Creek (German: Hühnerwasser) within the partially decarbonised opencast lignite mine Welzow Süd in NE Germany (details in Gerwin et al. 2009) (supplementary material Appendices A and B, Fig. B1). For the present study, we used annual data on quantitative plant surveys from 107 single plots of 25 m2 (excluding a number of unvegetated plots in 2005 and 2006) (Fig. A1, details in Zaplata et al. 2013). For each plot, at the beginning of the study, initial topsoil parameters (pH, organic carbon (OC) and the percentage of sand in the soil) were determined. We combined these soil variables using principle coordinates analysis and used the dominant eigenvector (EVS), which explained 98% of variance, in further analyses. The complete set of raw data and the Pearson correlation coefficients between the soil variables are contained in supplementary material Table A1 of Appendix A and Table B1 of Appendix B.

For each species, we estimated the abundance of each species by the cover degree (C) according to a modified Londo scale (Londo 1976; 0.1: ≤ 0.1%; 0.5: > 0.1–0.5%; 1: > 0.5–1%; 2: > 1–2%, in 1% steps up to 10; 15: > 10–15%, in 5% steps up to 30; 40: > 30–40%, in 10% steps up to 100). As species can overlap in space, the sum of all Londo values might exceed 100%. Bryophyta and Marchantiophyta were not identified to lower taxonomic levels. In total, we found 168 morphospecies during the seven study years. The complete data of species identities and cover degrees of all study years used in this study are contained in Table A1 and in Ulrich et al. (2014, 2016).

We used the Leda (Kleyer et al. 2008) and BioFlor (Klotz et al. 2002) databases and compiled one plant morphological trait (specific leaf area, SLA), one reproductive trait (seed number) and one dispersal trait (seed weight). Previous analysis had shown that the total expression of these traits showed a high small-scale spatial and temporal variability, although the underlying drivers of this variability remained obscure (Ulrich et al. 2014, 2017). Therefore, we have chosen these traits to answer our starting questions. In addition, we compiled data on three environmental demand traits using Ellenberg indicator values (Ellenberg et al. 1992) (light, soil moisture, soil fertility). The traits served as proxies to the changing environmental conditions during succession. We combined these traits using principle coordinates analysis and used the dominant eigenvector (EVD), which explained 94% of variance, in further analyses. Raw data and pairwise correlations of these variables are contained in Tables A3 and B1.

For the present analysis, we constructed three data matrices, a 168 species × 3 traits matrix T containing for all 168 species the raw values of the three traits, a 168 species × 642 plot/year matrix M containing the cover data for all plot × consecutive study year combinations (6 pairs of consecutive years × 107 plots = 642), and a 6 soil variables × 642 plot/year matrix V, containing the soil parameters. These three matrices are provided in Appendix A.

To quantify functional diversity, we used the functional attribute diversity (FAD) of Walker et al. (1999) being the normalized sum of all pairwise functional trait dissimilarities between the species of a community. FAD is therefore a community-wide measure of trait space obtained from the trait expression of each species. To obtain FAD, we first calculated for each of the three traits the 168 × 168 species × species Euclidean dissimilarity matrix D. The inner product A = DM then provides the FAD values at each of the 642 plot/year combinations. As the matrix A accounts for differences in species cover our approach is equivalent to the abundance-based approach of Schmera et al. (2009).

We used two approaches for subsequent data analyses. First, we compared observed FAD values with those obtained from a statistical null standard. This was obtained from a reshuffling of species names across trait values breaking of the cover-trait combinations. The associated null distributions were obtained from 200 randomizations. For comparisons we used standardized effect sizes \({SES}_{j}= \frac{{y}_{j}-\mu }{\sigma },\) where yj denotes the observed parameter and μ and σ are the respective mean and standard deviation of the null model distribution.

Second, Ulrich et al. (2022) used the Price approach (Eq. 1) and proposed to partition the difference in the expression of a focal variable (here FAD) between two communities A and B into five parts. In the present case, ‟community” refers to the plants occurring in two plots separated in space or time. The partition then takes the form

$$\Delta FAD=E({FAD}_{A})\Delta S+{C}_{B}\mathrm{E}(\Delta t\Delta p)+{C}_{B}\mathrm{E}(\Delta t)+{C}_{B}\mathrm{E}(\Delta p)+f\left(\Delta C\right)E({FAD}_{A})$$
(2)

containing terms quantifying the (1) Difference in species richness (∆S). (2) Average difference in trait expression (∆t). (3) Difference in relative cover (∆p). (4) Combined effect of differences in relative cover and trait expression (∆t∆p) and (5) Effect of differences in cover (∆C).

In Eq. 2, E denotes the statistical expectation (mean) in community A of the average value of FADA, the combined and pure changes in trait expression (∆t), and the relative cover p (∆p). CB symbolizes the total cover in community B. The proof of Eq. 2 is given in Ulrich et al. 2022. We note that these partitions, although being additive, are not always linearly independent. Particularly, the present richness partitions ∆S for seed number and seed mass were moderately correlated with the relative (∆p) and absolute (∆C) cover partitions (Table B2). Therefore, these co-variances must be taken into account when interpreting these partitions.

Here, we used FAD as the variable of interest and consequently, ∆FAD refers to the total change in FAD in a given plot between two consecutive years. We calculated ∆FAD and its partitions for all combinations of consecutive study years and plots (6 × 107 = 642). This approach allowed a full separation of the influences of plot characteristics at each point in time (richness, cover, environment) and the temporal change during succession.

To answer to our first and second starting questions, we used general linear mixed modelling with robust error estimation as implemented in SPSS 28. The changes in FAD and its Price partitions from six pairs of consecutive study years in the 107 plots (6 × 107 = 642 individual changes) served as response vectors. Predictors were species richness (SA), total cover (CA), total trait expression (TA), and the dominant eigenvectors for the soil variables (EVS) and plant demands (EVD) for community A. These predictors entered the resulting models as fixed variables. We repeated this for each trait variable leading to a total of 3 traits × 6 Price partition variables = 18 single models each having 7 predictors. Study year served as random qualitative predictor because each year can be seen as a sample containing a random climatic component. Predictor and response variables were Z transformed prior to analysis. The environmental predictor variables used in the models were at most moderately corrected, and multicollinearity should not affect the parameter estimates (Appendix B, Table B1). We note that the plots were temporarily not independent and significance levels may be inflated. However, this non-independence, expressed by the changes in trait space within each plot, is exactly what we wanted to study with this approach. Consequently, we focus on general trends and effect sizes.

To answer our third starting question, we used three approaches. First, Eq. 2 allows for a comparison of the proximate factors that affect temporal changes in FAD, as well as for separate analyses of the ultimate covariates. To assess the relative impact of each Price partition in space, we used average values and respective variance—mean ratios (VMR = σ2/μ, where σ and μ are the respective standard deviations and mean values) calculated from all 107 single plots. VMR = 1 indicates a Poisson random distribution. Finally, we calculated the third and fourth statistical moments (skewness and kurtosis) to assess the shape of the distributions of ∆FAD and its partitions across study plots.

Results

The importance of single partitions differed between traits, although the change in trait expression ∆t had in all three cases the highest impact on the change in FAD as quantified by ∆FAD (Table 1). With respect to seed numbers and SLA, differences in relative abundance (∆p) and total cover (∆C) always had a comparatively low impact on ∆FAD, whilst the change in seed mass was strongly influenced by the change in richness (∆S, Table 1).

Table 1 Moments of Price partitions (Eq. 1) of 107 plot × 6 study year combinations (N = 642) for trait diversity: arithmetic mean (μ), standard deviation (σ), skewness (γ) and kurtosis (δ)
Full size table

The average species richness of the plots increased constantly during the seven years of succession, whilst the average plant cover of the plots reached a plateau after four years (Fig. 1a, b). The variance mean ratio (VMR) of plot species richness scattered in all study years around unity indicating a Poisson random distribution of richness across the study plots (Fig. 1c). In turn, VMR of plant cover was from 2006 onwards well above unity and reached a plateau in 2009 (Fig. 1c).

Fig. 1
figure 1

Species richness (a), total cover (b), the average variance/mean ratios (σ2/μ) of species richness (dark grey dots) and total plant cover (light grey) (c) and average FAD values on 107 plots of early plant succession from 2005 to 2011 (d). Boxplots in a and b: median, quartiles, maxima and minima and outliers. Colours in d: red: no. of seeds, blue: seed size, green: specific leaf area (SLA)

Full size image

Average FAD steadily increased during the first 4 years of succession and reached a plateau in 2009 (Fig. 1d). We did not find respective significant annual trends in the annual changes of FAD (∆FAD) and its partitions (Fig. B2). However, two other patterns were obvious. First, variability in ∆FAD and consequently in the partitions was small during the first 3 years (Fig. B2). There was a clear step in variability from 2007 to 2008 and from 2008 onwards the variability in the annual change of FAD amongst the study plots was comparatively high (Fig. B2).

Fig. 2
figure 2

Skewness values (γ∆x, where ∆x is one of the partitions) amongst 107 plots of the SES transformed Price partitions for ∆FAD (dark blue), ∆S (yellow), ∆t∆p (grey), ∆t (red), ∆p (blue), ∆C (green). a Seed numbers, b seed sizes, c specific leaf area (SLA). Error bars denote bootstrapped 95% confidence limits

Full size image

We did not find plot-specific patterns in the change of FAD (Table B3). Pearson correlations of the plot ∆FAD values between study years showed only low correlations, explaining less than 10% of variance (Table B3). There was one notable exception. For all three traits, the change in FAD between 2007 and 2008 appeared to be synchronized amongst plots, as indicated by the higher plot correlations of r ≥ 0.4.

The corresponding analysis of higher moments of the distributions of the Price partitions across all study years revealed a high spatial variance in the average change of the FAD and its partitions between the study plots (Table 1). The coefficient of variation was for all three traits and all of the partitions well above 1 or below − 1 (Table 1). The high skewness values (γ >  > 0) indicate that this variability in ∆FAD stem from a small number of plots with comparatively high change (Table 1), whereas plot identity changed from year to year (Table B3). The detailed analysis of the spatial distribution of ∆FAD at its partitions revealed temporal trends and trait-specific pattern (Fig. 2). The respective distributions deviated significantly from zero for most partitions and were mostly strongly positively skewed, consistent with a small number of plots with very high annual changes in functional diversity (Fig. 2). In terms of seed number (Fig. 2a) and SLA (Fig. 2c), this pattern was most pronounced in the change of FAD from 2007 to 2008, where skewness values peaked. In line with these findings, kurtosis values were in all cases well above δ = 3, indicating a more flat distribution than expected from a normal distributions of FAD change across plots (Table 1).

The average interannual changes in FAD (∆FAD) of seed numbers and SLA, but not of seed size, varied significantly between the study years (Fig. 3). The Price partitioning revealed that the weak and insignificant influences of the model predictors on ∆FAD resulted from respective contrasting effects on these partitions (Fig. 3). Particularly, species richness had insignificant effects on ∆FAD for all three traits (Fig. 3). The reason for this was a trade-off between changes in richness relative to cover (Fig. 3). Similar trade-offs occurred for TA (seed number, seed size and SLA) and CA (seed number). Initial soil characteristics and spatial plot distance did not significantly influence the temporal changes in FAD (Fig. 3).

Fig. 3
figure 3

Parameters values and respective standard errors of 18 single general linear mixed models (study year as random effect) with ∆FAD (a, g, m), ∆S (b, h, n), ∆t∆p (c, i, o), ∆t (d, j, p), ∆p (e, k, q) and ∆C (f, l, r) as trait diversity response variables, and study year, the dominant eigenvector of the Euclidean plot distance matrix (Space), total cover (CA) and species richness (SA), total trait expression (TA) and the initial soil conditions (Soil) and the Ellenberg plant demand traits (Demand) as predictors. a–f: Average seed numbers, g–l: average seed size, m–r: SLA. Dark grey bars: p < 0.01, light grey: p > 0.01. N = 642

Full size image

Discussion

We have confirmed that an additive partitioning of the change in an aggregate community trait, here FAD, into major components is able to reveal the effects of abiotic and biotic drivers behind this change. In ecology, Price partitioning has first been applied to study the richness effect for ecosystem functioning (Loreau and Hector 2001), whilst Fox (2006), Isbell et al. (2018) and Ulrich et al. (2022) extended the approach to frameworks for the analysis of changes in community properties. Critics have argued that Price partitions are not statistically independent (Pillai and Gouhier 2019) and lack a simple, unequivocal interpretation, or even become meaningless (Van Veelen 2020). However, the present work as well as Ulrich et al. (2022) and Loreau and Hector (2019) have shown that non-independence (cf. Table B2 for the present data) does not exclude application and that partitions have meaningful interpretation. Ladouceur et al. (2022) demonstrated how careful use can shed light on the basic processes behind plant biomass dynamics, particularly the effects of species gains and losses. In this respect, Price partitions resemble variance partitioning methods and non-orthogonal principle coordinates analyses whose components can be interpreted despite of being correlated.

Our first starting question asked about the major drivers of functional diversity during early plant succession. This was for all three traits, the change in the trait space of individual species followed by the richness effect (Table 1). Changes in relative and absolute cover were of comparatively minor importance (Table1). If community assembly were temporarily neutral according to the ecological drift model (Hubbell 2001), the species functional equivalence (identical trait spaces) should have caused major effects of richness and cover but not of trait space. Our results, therefore, indicate that the change in functional diversity during early succession (as measured by FAD) depends mainly on species richness and the associated variability in trait expression due to species turnover. The associated changes in community composition and cover seem to be of minor importance. Therefore, we argue that functional diversity in plant communities is mainly species driven.

In this respect, much work has been devoted to the dynamics of trait expression during plant succession. It is generally accepted that local habitat filters in combination with the initial habitat conditions select for certain traits in a temporally ordered sequence leading to predictable changes in community composition at each successional stage (Weiher and Keddy 1999; Zaplata et al. 2013; Chang and Turner 2019, but see Götzenberger et al. 2012). In this respect, Douma et al. (2012) have shown that these changes can be related to changing light availability and initial abiotic soil conditions that shift the expression of major plant functional traits. Here, we extend this knowledge to the change in functional diversity and therefore to the total trait space spanned by a community at a certain successional state.

We found significant differences in FAD change (∆FAD) amongst study years, although these were only weakly linked to initial soil conditions and only for SLA significantly related to plant habitat demands (Fig. 3). Unfortunately, we lack small-scale light conditions at the plot level. The partitioning confirmed that this lack was not due to trade-offs between richness (∆S), composition (∆t∆p, ∆p), and cover (∆C) partitions (Fig. 3). Our results instead point solely to the levels of richness, cover, and functional diversity that determine the degree of annual change in the course of succession (Fig. 3). Under this hypothesis, the change in functional diversity of the three traits studied here would resemble a Markov chain process applied to successional processes and community assembly (Horn 1975; Usher 1981). In classical Markov chain models, each subsequent state is solely determined by the previous state (in the present case, richness, cover and functional diversity) in combination with a specific transition rule. In this respect, these models resemble the present Price analysis where the change in trait value (here FAD) of two consecutive temporal states were compared and the partitions contained (except for the fixed multiplicator CB) the initial state values only. Classical Markov chain models with fixed transition probabilities have been found to be too rigid in most cases (e.g. Usher 1981), whilst extended models with adjustable transition probabilities, particularly those with unobserved (hidden) states (Rabiner and Juang 1986), appeared to be suited to describe patterns of temporal and spatial species replacements (e.g. Logofet et al. 1997; McClintock et al. 2020). Forthcoming work has to show whether and how Price partitions can be reformulated within a hidden Markov framework. The associated transition matrices and their covariates would then inform about the quantitative parameters behind the succession process.

Classical Markov models are consistent with neutral community assembly (Steiner and Tuljapurkar 2012), whereas extended and hidden models include various degrees of niche and filter mechanisms (McClintock et al. 2020). The fact that SA, CA, and TA explained generally less than 50% of the variance in ∆FAD indicates that other factors, presumably niche and filter effects, are important. This result does not point to simple neutral mechanisms of community assembly and successional change. Our results thus corroborate Ulrich et al. (2014), who already reported that filter processes gained importance for community assembly and FAD, particularly during the first four years of succession. This filter effect acted mainly via a reordering of species dominances (the ∆p and ∆t∆p partitions, Table 1).

Answering to our second starting question, the present Price partitions gave insight into the trade-offs between the drivers of the change in FAD (Fig. 3). For seed numbers and SLA, we found contrasting effects of the richness (∆S) and relative cover (∆p) components and to a lesser degree between ∆p and the covariance between trait and relative cover (∆t∆p) (Fig. 3). This result was unexpected as the temporal increase in richness (Fig. 1) should mainly affect rare species, whilst the total FAD should be largely unaffected. Apparently, this was not the case and newly colonizing species immediately changed the dominance orders with major influence on FAD (Fig. 3). This was most apparent in the FAD change from the third to the fourth year of succession (2007/2008, cf. Zaplata et al. 2013). Incoming species tended to increase FAD, whilst the reordering of the cover (abundance) distributions mediated this trend leading to a stabilizing effect on the annual changes in FAD (Fig. 3). Indeed, classical and current work (Bazzaz 1975; Yin et al. 2018) demonstrated that community evenness increases during plant succession. The higher proportions of species with intermediate cover should stabilize trait spaces as the probability of rearrangements within this group of intermediate species increases in comparison to rearrangements within the groups of the most abundant and the rare species. Soliveres et al. (2015) called this the proper middle class effect, because these species might take over dominance in cases where dominants face population breakdown. This effect should particularly hold for FAD being calculated from the product of trait expression and cover. Consequently, our findings support the hypothesis that the temporal variability in relative cover mediates the increase in total FAD due to increasing species richness in the course of plant succession.

Finally, we asked whether the importance of each partition predictably changes in time during plant succession. This was not the case (Fig. B2). We did not see consistent annual trends for each partition. This lack of trend indicates that the annual changes within the study system are spatially and temporarily point patterns triggered solely by the current state of the environmental and the composition of the local biotic community. Of course, in prior work we found a steady increase in plant species richness and cover (cf. Fig. 1), as well as a trend towards increasing spatial species segregation (increasing β-diversity amongst plots; Zaplata et al. 2013), but also a high small-scale variability in richness and phylogenetic and functional diversity (Ulrich et al. 2014). The fact that environmental and especially soil conditions were similar after catchment construction suggests that initial small local stochastic processes of plant colonization or establishment from the soil seed bank (Zaplata et al. 2010, 2011) may trigger this large small-scale variability in plant community composition. Prior, Wiegleb and Felinks (2001) already found that plant succession in comparable post-mining landscapes did not follow predictable chronosequences. Here, we have shown that this result also holds for the temporal change in functional diversity.

This result is strongly corroborated by the high spatial variability in ∆FAD and its partitions (Tables 1, B3, Figs. 3, B2). Ulrich et al. (2014) already pointed to the high small-scale spatial variability of phylogenetic diversity in this system. The high positive skewness and kurtosis values in the ∆FAD distribution across plots particularly point to a few plots in each year with extraordinary high annual temporal variability in FAD. Importantly, these plots varied in time as indicated by the low correlations of ∆FAD between plots, particularly amongst consecutive years (Table B3). Therefore, we could not identify highly variable and more predictable plots, so that the successional process was even more determined by hidden factors. Consequently, the annual change in functional diversity again appears to be a process without memory, as is typical for Markov processes with variable transition probabilities. However, our study design could not clarify whether accounting for plot-specific annual characteristics manifested in species composition could elucidate the underlying patterns. Subsequent plot-specific analyses should uncover predictable successional patterns in this system.