Biodiversity stabilizes plant communities through statistical-averaging effects rather than compensatory dynamics

Empirical data

We applied our theory to two datasets (Table 1): the plant survey dataset and the biodiversity-manipulated experiment dataset. The plant survey dataset contains nine sites of long-term grassland experiments across the United States (see also Hallett et al.¹⁰, and Zhao et al.²³). Five of nine sites are from the Long Term Ecological Research (LTER) network (see Table 1). Plant abundances were measured either as biomass or as percent cover. In percent-cover cases, summed values can exceed 100% due to vertically overlapping canopies. All sites were sampled annually and were spatially replicated. We only used data of the plant survey dataset from unmanipulated control plots. Methods for data collection were constant over time.

The biodiversity-manipulated experimental dataset comprises two long-term grassland experiments, BigBio and BioCON, at the Cedar Creek Ecosystem Science Reserve. Both experiments directly manipulated plant species number (1, 2, 4, 8, 16 for BigBio; and 1, 4, 9, 16 for BioCON). BioCON also contains different treatment levels for nitrogen and atmospheric CO₂, but here only data from the ambient CO₂ and ambient N treatments were used. We excluded plots with only one species. BigBio comprises 125 plots over 17 years, and BioCON comprises 59 plots over 22 years (Table 1).

Theory

Let x_i(t) denote the biomass of species i = 1, …, S at time t = 1, …, t and let μ_i = mean (x_i (t)), σ_i = ({{mbox{sd}}})(x_i (t)), and ({v}_{i}={sigma }_{i}^{2}) be the mean, standard deviation and variance of species i, computed through time. Let v_ij = cov (({x}_{i}left(tright),, {x}_{j}left(tright))) be the covariance, through time, of the dynamics of species i and j. Let x_tot (left(tright)={sum }_{i}{x}_{i}(t)), ({mu }_{{{mbox{tot}}}}={sum }_{i}{mu }_{i}), ({v}_{{{mbox{tot}}}}={sum }_{i,j}{v}_{{ij}}), and ({{{{{{rm{sigma }}}}}}}_{{{{{{rm{tot}}}}}}}=sqrt{{v}_{{{{{{rm{tot}}}}}}}}). When population time series are uncorrelated, ({v}_{{{{{{rm{tot}}}}}}}={sum }_{i}{v}_{i}).

As defined previously^10,15, community stability is the inverse coefficient of variation of ({x}_{{{mbox{tot}}}}left(tright)), ({S}_{{{{{{rm{com}}}}}}}={mu }_{{{{{{rm{tot}}}}}}}/{sigma }_{{{{{{rm{tot}}}}}}}). Population stability is the inverse of weighted-average population variability⁹, ({sum }_{i}frac{{mu }_{i}}{{mu }_{{{{{{rm{tot}}}}}}}}{{CV}}_{i}={sum }_{i}frac{{mu }_{i}}{{mu }_{{{{{{rm{tot}}}}}}}}frac{{sigma }_{i}}{{mu }_{i}}={sum }_{i}frac{{sigma }_{i}}{{mu }_{{{{{{rm{tot}}}}}}}}), i.e, ({S}_{{pop}}={mu }_{{{{{{rm{tot}}}}}}}/{sum }_{i}{sigma }_{i}). The ratio of community stability over population stability is the Loreau-de Mazancourt asynchrony index¹⁴, Φ = ({sum }_{i}{sigma }_{i}/{sigma }_{{{{{{rm{tot}}}}}}}), so that

Now we suppose a hypothetical community with the same species-level variances and means as the original community but with species covariances equal to zero. Then, (1) becomes S_{com_ip} = (SAE)S_pop, where ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{ip}}}}}}}=frac{{mu }_{{{{{{rm{tot}}}}}}}}{sqrt{{sum }_{i}{v}_{i}}}=frac{{mu }_{{{{{{rm{tot}}}}}}}}{sqrt{{sum }_{i}{sigma }_{i}^{2}}}) is the value of community stability in the case of uncorrelated or independent populations and SAE is the component of Φ due to statistical averaging (here, “ip” stands for “independent populations”). The equation S_{com_ip} = (SAE)S_pop can be interpreted as a definition of SAE. We then have

The compensatory effect is then the rest of Φ, i.e.,

Considering the classic variance ratio ({{{{{rm{varphi }}}}}}=frac{{V}_{{{{{{rm{tot}}}}}}}}{{sum }_{i}{V}_{i}}=frac{{sigma }_{{{{{{rm{tot}}}}}}}^{2}}{{sum }_{i}{sigma }_{i}^{2}}), our CPE is (1/sqrt{varphi }). Values CPE > 1 (respectively, <1) correspond to greater (resp., lesser) community stability than would be expected if dynamics of different taxa were independent, reflecting compensatory (resp., synchronous) dynamics. We also have

where ({p}_{i}={sigma }_{i}/{sum }_{i}{sigma }_{i}).

In the Introduction and Results, we characterized SAE as relating to statistical mechanisms, but Eq. (4) shows that the strength of SAE depends on the evenness of species variances, which may be influenced by species biology. In fact, ({SAE}=frac{{S}_{{{{{{rm{com}}}}}}_{{{{{rm{ip}}}}}}}}{{S}_{{{{{{rm{pop}}}}}}}}) is a comparison between the two values that community stability would take in the scenarios of independent populations (({S}_{{{{{{rm{com}}}}}}_{{{{{rm{ip}}}}}}})) and entirely synchronous populations (({S}_{{{{{{rm{pop}}}}}}};) see elsewhere¹⁴ for a proof that ({S}_{{{{{{rm{pop}}}}}}}) is the value community stability would take if populations were perfectly synchronous), and so may be best characterized as a statistical effect, the strength of which depends on evenness of variances (and therefore on biology). We define ({{SAE}}_{{{{{{rm{even}}}}}}}) to be the value SAE would take in the case of perfect evenness of species variances, i.e., if all species variances were identical (that means σ₁ = σ₂ = … = σ). We can then define the evenness effect, ({EVN}=frac{{SAE}}{{{SAE}}_{{{{{{rm{even}}}}}}}}), so that the equation ({SAE}={EVN}times {{SAE}}_{{{{{{rm{even}}}}}}}) separates out biological/evenness effects from purely statistical effects. ({{SAE}}_{{{{{{rm{even}}}}}}}) is a purely statistical version of statistical averaging because it only captures the stabilizing effects of averaging independent random variables, being uninfluenced by evenness. It can straightforwardly be shown that ({{SAE}}_{{{{{{rm{even}}}}}}}=sqrt{S}): applying Eq. (4) to the scenario of equal population variances, one sees that ({{SAE}}_{{{{{{rm{even}}}}}}}=sqrt{frac{1}{{{sum }_{i}left(1/Sright)}^{2}}}=sqrt{S}). This is the maximum possible value of SAE for a given species richness, S.

It is possible to explicitly quantify the effects of the evenness of variability. Using the formula ({{log }}left({S}_{{{{{{rm{com}}}}}}}right)={{log }}left({SAE}right)+{{log }}left({CPE}right)+{{log }}({S}_{{{{{{rm{pop}}}}}}}),) the main observation of our study has been that regression slopes of ({{log }}left({SAE}right)) against ({{log }}left(Sright)) tend to have substantially positive slopes, whereas regression slopes of ({{log }}left({CPE}right)) and ({{log }}({S}_{{{{{{rm{pop}}}}}}})) against ({{log }}left(Sright)) tend to have negative or only slightly positive slopes; and hence regression slopes of ({{log }}left({S}_{{com}}right)), which are sums of the other slopes, tend to be positive primarily because of the positive relationship between log(SAE) and log(S). But, in another decomposition, we can also write ({{log }}left({S}_{{com}}right)={{log }}left({{SAE}}_{{{{{{rm{even}}}}}}}right)+{{log }}left({ENV}right)+{{log }}left({CPE}right)+{{log }}({S}_{{{{{{rm{pop}}}}}}})=frac{1}{2}{{log }}left(Sright)+{{log }}left({ENV}right)+{{log }}left({CPE}right)+{{log }}({S}_{{{{{{rm{pop}}}}}}})). We can observe that, because all our regressions of ({{log }}left({SAE}right)) against ({{log }}left(Sright)) had slope less than 1/2 (Figs. 2e and 4e), regressions of ({{log }}left({EVN}right)) versus ({{log }}left(Sright)) will have negative slope for our empirical systems. Therefore ({{log }}left({S}_{{com}}right)) tends to depend on ({{log }}left(Sright)) with a “baseline” slope of 1/2, for purely statistical reasons (the ({{log }}left({{SAE}}_{{{{{{rm{even}}}}}}}right)) effect), with effects of evenness (({{log }}left({ENV}right))), compensatory dynamics (log(CPE)) and population variability (log(S_pop)) tending to reduce this slope. This again supports the conclusion that statistical averaging is the main reason for positive diversity-stability relationships in our data, now using a definition of statistical averaging that is purely statistical. We used SAE instead of ({{SAE}}_{{{{{{rm{even}}}}}}}) in Results to represent the concept of statistical averaging because it is a quantification of ideas which were previously present in the literature²⁰, although it can be viewed instead as a combination of statistical averaging effects and effects of evenness of species variances.

Note that Loreau and de Mazancourt¹⁴ got synchrony as 1/S for the special case of identical ({sigma }_{i}) and zero-correlation, a value which seems to contrast with our result, (sqrt{S}), described above. But Loreau and de Mazancourt¹⁴ measured community-wide synchrony via a statistic (theta=frac{{sigma }_{T}^{2}}{{left({sum }_{i}{sigma }_{i}right)}^{2}}). We can make further transformation: (theta=frac{{({sigma }_{T}/{mu }_{{tot}})}^{2}}{{left({sum }_{i}{sigma }_{i}/{mu }_{{tot}}right)}^{2}}=frac{{{CV}}_{{com}}^{2}}{{left({sum }_{i}frac{{mu }_{i}}{{mu }_{{tot}}}bullet frac{{sigma }_{i}}{{mu }_{i}}right)}^{2}}=frac{{{CV}}_{{com}}^{2}}{{{CV}}_{{pop}}^{2}}). Here CV_pop² means the square of the weighted average of CV_i. That is, θ measures a ratio of community-level CV² to population level CV². But in our system, our asynchrony measures a ratio of community-level stability (1/CV) to population-level stability. These differences in notational and terminological choices explain why we got (sqrt{S}) but Loreau and de Mazancourt¹⁴ got 1/S in this special case – our results are actually consistent with those of Loreau and de Mazancourt when notational choices are accounted for.

Further decomposition of CPE

In previous sections, we decomposed asynchrony into a statistical averaging effect (SAE) and a compensatory effect (CPE) via constructing a surrogate or hypothetical community stability value S_{com_ip} by quantifying what community stability would be if all species were to fluctuate independently and thus species covariances equal zero. By comparing the stability of the real community, S_com, to the hypothetical community stability under the assumption of species independence, S_{com_ip}, we then computed CPE as S_com/S_{com_ip}. But compensatory dynamics among species in a community can arise from at least two distinct mechanisms: (1) differential response of species to environmental perturbations, and (2) direct or indirect interactions among species¹⁹. Correspondingly, we can decompose CPE into two parts: compensatory effects due to differential responses to environmental perturbations (CPE_env) and compensatory effects due to species interactions (CPE_int). We here initiate that decomposition.

To separate CPE_env from CPE_int, we construct another surrogate/hypothetical community stability value, called ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}), for which the effects on community stability of components of species covariances due to species interactions are eliminated, keeping only components due to common species responses to environmental perturbations. In the following section “Surrogate community stability from different plots”, we establish surrogate community stability values for the systems of the plant survey dataset by means of time series from other plots from the same system. Time series from distinct plots may covary, but if they do, it is likely because of common environmental influences rather than species interactions. In the section “Surrogate community stability from monocultures”, we establish surrogate community stability values for the systems of the biodiversity-manipulated experimental dataset by means of monocultures which were established alongside the multi-species plots. Monocultures in separate plots may covary, but if they do, it is likely again because of common environmental influences rather than species interactions. Both approaches lead to approximations of the value community stability would take if species interactions were eliminated, and therefore ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}={{{CPE}}_{{{{{{rm{env}}}}}}}times {SAE}times S}_{{{{{{rm{pop}}}}}}}={{{CPE}}_{{{{{{rm{env}}}}}}}times S}_{{{com}}_{{{{rm{ip}}}}}}). Thus, we calculate CPE_env by

The compensatory effect from species interactions is then the rest of CPE, i.e.

Both CPE_env and CPE_int can be either greater than or less than 1, corresponding to the fact that joint species responses to environmental perturbations (CPE_env) and species interactions (CPE_int) can be either stabilizing or destabilizing. Thus, we have decomposed community stability, S_com, into four parts (Supplementary Fig. 5): population stability, S_pop; the statistical averaging effect, SAE; compensatory effects due to joint species responses to environmental perturbations, CPE_env; and compensatory effects due to species interactions, CPE_int.

In the Introduction and Results, the term “portfolio effect” was used for the stabilizing effect of biodiversity from independent fluctuations (zero correlation) of species abundances through time²⁰, measured by SAE. However, other studies have adopted an alternative definition of portfolio effects, allowing for non-zero correlations between species induced by common responses to the environment, but excluding components of correlations resulting from species interactions^12,19,21. In our framework, this alternative definition of portfolio effects could be represented by SAE × CPE_env (Supplementary Fig. 5). We prefer to group CPE_env together with CPE_int rather than with SAE, because SAE corresponds to primarily statistical mechanisms, whereas CPE_env involves species responses to the environment, which we regard as containing biological information. But other interpretations and classifications can certainly be reasonable, as demonstrated by earlier studies which defined portfolio effects in a way that corresponds to our product SAE × CPE_env^12,19,21. In the following sections, we decompose diversity-stability relationships into components that come from the relationship of each of the four terms of Supplementary Fig. 5 to diversity.

Surrogate community stability from different plots

In this section, we explain how we construct surrogate communities that apply to observational data without monoculture data, such as our plant survey dataset. In doing so, we calculate an approximate hypothetical/surrogate community stability index, S_{com_sur}, using populations from different mixture plots. The general idea is to construct a surrogate covariance matrix,

which will be used to calculate S_{com_sur}. Again using x_i(t) to denote the biomass of species i = 1, …, S at time t = 1, …, t in a focal plot, the ith diagonal element of the surrogate covariance matrix, V_ii, is the variance of x_i(t). The off-diagonal element C_ij equals the covariance of x_i and Z_j, where ({Z}_{j}left(tright)=frac{{z}_{j}left(tright)-{bar{z}}_{j}}{{{{{{rm{sd}}}}}}left({z}_{j}right)}times {{{{{rm{sd}}}}}}left({x}_{j}right)+{bar{x}}_{j}) is a surrogate time series, with z_j representing a time series of species j taken from another plot. When more than one other plot contained species j, we randomly chose one. The time series ({Z}_{j}left(tright)) was constructed to have the same mean and variance as x_j(t), but its correlations with x_j(t) reflects correlations between two non-interacting populations of species i and j from distinct plots. We then calculated ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}=frac{{bar{X}}_{{{{{{rm{tot}}}}}}}}{sqrt{mathop{sum}limits_{{ij}}{C}_{{ij}}}}), and we calculated CPE_env and CPE_int from Eqs. (5) and (6). This method was applied only to the plant survey dataset, where plots at a site within the plant survey dataset were considered replicates. In the biodiversity-manipulated experimental dataset, different plots of the same diversity had different artificially maintained species composition.

Surrogate community stability from monocultures

In this section, we explain how we construct surrogate communities that apply to a biodiversity experiment with monoculture data. Given a community consisting of species i = 1, …, S, recall that x_i(t) denotes the biomass of species i = 1, …, S at time t = 1, …, t. Suppose monocultures, y_i(t), were maintained for each species over the same time period. When more than one monoculture was maintained for a species, we randomly selected one. The hypothetical/surrogate community stability value ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}) was computed as ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}=frac{{bar{Y}}_{{{{{{rm{tot}}}}}}}}{{{{{{rm{sd}}}}}}left({Y}_{{{{{{rm{tot}}}}}}}right)}), where ({Y}_{{{{{rm{tot}}}}}}left(tright)={sum }_{i}{Y}_{i}(t)) and ({Y}_{i}left(tright)=frac{{y}_{i}left(tright)-{bar{y}}_{i}}{{{{{{rm{sd}}}}}}left({y}_{i}right)}times {{{{{rm{sd}}}}}}left({x}_{i}right)+{bar{x}}_{i}). The time series ({Y}_{i}left(tright)) were constructed to have the same mean and variance as the time series x_i(t), but with correlations equal to those between the corresponding monocultures. CPE_env and CPE_int were then calculated from Eqs. (5) and (6). We repeated this process 100 times, randomly selecting from among available monocultures for each species on each repeat calculation, and taking means to get CPE_env and CPE_int. This method was only used for the biodiversity-manipulated experimental dataset because monocultures were not available for the plant survey dataset.

It is worth noting that our construction of the surrogate community assumes that in the absence of species interactions, species in the mixture should exhibit the same correlation as those in monocultures. While this assumption may hold if population dynamics were solely influenced by environmental stochasticity, the influence of demographic stochasticity may lead to a lower correlation in the mixture than in monocultures due to the relatively lower population size in the mixture. In other words, our surrogate community might overestimate the correlation between species and thus underestimate ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}) (i.e., the expected stability of a community without interactions). Moreover, the possible underestimation of ({S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}) could be more pronounced for more diverse communities (which have substantially lower population sizes in mixtures than monocultures), thus possibly creating an artificially negative correlation between species richness and CPE_env ((={S}_{{{{{{rm{com}}}}}}_{{{{{rm{sur}}}}}}}/{S}_{{{{{{rm{com}}}}}}_{{{{{rm{ip}}}}}}})), as well as a positive correlation between species richness and CPE_int ((={S}_{{{{{{rm{com}}}}}}}/{S}_{{com}_{sur}})). So, the patterns in Supplementary Fig. 7a, b might be partially explained by such artificial effects. However, how to separate out such artificial effects from other processes was by no means obvious³⁰. We note, importantly, that such effects, even if they exist, should not affect our main conclusion. Indeed, the positive correlation between species richness and SAE × CPE_env (Supplementary Fig. 7d) would be even stronger if such effects existed and were accounted for.

Statistical analysis

For the plant survey dataset, the values of each of our theoretical quantities were averaged across all plots for each site and then relationships between site averages were assessed with regression. For the biodiversity-manipulated experimental dataset, theoretical quantities were calculated for individual plots because species richness was manipulated at the plot level. All theoretical quantities were log₁₀ transformed prior to generating plots and preforming regressions, except for the Shannon index, because the Shannon index has a log transformation already embedded in its definition. All analyses were programmed in R 4.0.5³⁷.

To compute the confidence interval of CPE, we constructed bootstrapped species time series. For each plot, separately, the AAFT procedure³⁸, implemented in the surrog function in the R package wsyn, was used to generate bootstrapped species time series for the plot, for which time series were rendered uncorrelated while retaining the autocorrelation and marginal-distribution properties of the original data. Such a bootstrapped dataset instantiates the null hypothesis that species were unrelated while retaining other statistical properties of the data not related to species relationships. For each plot, 1000 bootstrapped datasets were computed, giving rise to 1000 bootstrapped quantities S_{com_ip}⁽ⁱ⁾ for i = 1, …, 1000. We then computed the quantities SAE⁽ⁱ⁾ = S_{com_ip}⁽ⁱ⁾/S_pop and CPE⁽ⁱ⁾ = S_com/(SAE⁽ⁱ⁾S_pop). Confidence intervals for CPE were computed using quantiles of the distribution CPE⁽ⁱ⁾.

Reporting summary

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

Source: Ecology - nature.com