Species richness stabilizes productivity via asynchrony and drought-tolerance diversity in a large-scale tree biodiversity experiment

Species richness promotes forest stability via enhanced asynchrony, which is positively related to drought-tolerance diversity.


Fig. S1. Principal component analysis (PCA) biplot for the trait-based definition of two drought-tolerance trait gradients.
(1) A gradient of stomata control that characterizes species as 'water spenders' if they down-regulate their stomatal conductance only at high levels of measured water pressure deficit (high VPDPOI and VPDMAXFIT), while 'water savers' are species that already down-regulate stomatal conductance at low water pressure deficits and have leaves characterized by a high a stomata density (high STODENS, STOIND; see Table S1 for a trait description). (2) A gradient of drought resistance based on the water potential at which 50% of xylem conductivity is lost (50) as key physiological trait that expresses a species resistance to water stress (2), with less negative values of 50 indicating a higher susceptibility to drought induced cavitation. The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few versus abundant stomata) and drought resistance (high versus low cavitation resistance). Note, that classic traits of the leaf economics spectrum (LES) (35) are associated with cavitation resistance in that species with high cavitation resistance also have traits commonly used to ascribe a conservative resource use strategy (tough leaves (LEAFT) and a high C/N ration (CN)) while low cavitation resistance is associated to a high resource acquisitiveness (high SLA) (33,36,43). We refer to this trait gradient therefore as 'resistanceacquisition' gradient. All traits (Table S1) were measured on site and used to calculate species level mean trait values, see refs. (33,34) for details. Species identity is shown as species code; see Table S3 for a detailed species list. Deciduous and evergreen species are coloured as black dots and green triangles, respectively. PCA was performed with the rda function in the vegan package version 2.5-6 (64). The principal components explained 31% (PC1) and 23% (PC2) of variation. Varimax rotated principal components (base R version 3.6. 2,70) were used in the analysis to achieve a good alignment of the two orthogonal trait gradients with the first and second PCA axis.

Fig. S3. Effects of population stability on community stability with increasing asynchrony.
The line is a linear mixed-effects model fit that shows a significant increase in community stability with population stability in mixtures (P<0.001). The thin 1-1 line shows where community stability is equal to average species-level population stability. Points are coloured according to their value with deeper red indicating increasing asynchrony. Grey bands represent a 95% confidence interval. See Table S2 for details on the fitted model. The model is the same as the one shown in Fig. 3D.   Fig. S4. Relationships between drought-tolerance diversity and asynchrony. Lines are linear mixed-effects model fits that show (a) significant increases in asynchrony with functional diversity of stomatal control (FD stomatal control; P<0.001) and (b) significant increases in asynchrony with functional diversity of resistance-acquisition (FD resistance-acquisition; P<0.001) in mixtures. Asynchrony ranges from 0 to 1. 0 represents complete synchrony and 1 represents complete asynchrony. Functional diversity calculated as abundance-weighted functional dispersion. Grey bands represent a 95% confidence interval. See Table S2 for details on the fitted models.

Fig. S5
. Relationships between drought-tolerance diversity and community stability. Lines are linear mixed-effects model fits that show (a) a marginally significant increase in community stability with functional diversity of stomatal control (FD stomatal control; P=0.058) and (b) a non-significant relationship of community stability with functional diversity of resistanceacquisition (FD resistance-acquisition; P=0.27) in mixtures. Functional diversity calculated as abundance-weighted functional dispersion. Grey bands represent a 95% confidence interval. See Table S2 for details on the fitted models.

Fig. S6. Relationships between CWMs of drought-tolerance traits and community stability.
Lines are linear mixed-effects model fits that show (a) a non-significant relationship of community stability with the CWM of stomatal control (CWM stomatal control; P=0.65) and (b) a non-significant relationship of community stability with the CWM of resistance-acquisition (CWM resistance-acquisition; P=0.88) in mixtures. Grey bands represent a 95% confidence interval. See Table S2 for details on the fitted models.

Fig. S7. Relationships between CWMs of drought-tolerance traits and population stability.
Lines are linear mixed-effects model fits that show (a) a non-significant relationship of population stability with the CWM of stomatal control (CWM stomatal control; P=0.52) and (b) a non-significant relationship of population stability with the CWM of resistance-acquisition (CWM resistance-acquisition; P=0.72) in mixtures. Grey bands represent a 95% confidence interval. See Table S2 for details on the fitted models.  S8. Hypotheses driven framework of the direct and indirect drivers of community stability in mixed-species tree communities. Directional arrows represent hypothesized causal relationships. We explored whether the data supported our first hypothesis through including indirect pathways that tested for effects of species richness on community stability that are mediated via asynchrony and population stability (7,28). We tested our second hypothesis through including indirect pathways that tested for effects of functional diversity of stomatalcontrol traits and functional diversity of resistance-acquisition traits on community stability through effects mediated via asynchrony (7,36,41). Similarly, we tested our third hypothesis through including indirect pathways that tested for the effects of the CWM of stomatal-control traits and the CWM of resistance-acquisition traits on community stability through effects mediated via population stability (28,36,38,41). As the experimental manipulation of species richness may directly affect the functional diversity of a community (39), we included pathways from species richness to functional diversity of stomatal control and functional diversity of resistance-acquisition. We further included direct pathways from the diversity facets and the CWMs of both trait gradients to community stability, to test for remaining effects not mediated by asynchrony or population stability. This further allowed us to test our second and third hypothesis separately through either including asynchrony or population stability (Figs. S9-S10). In the absence of population stability, these direct pathways could for example account for performance enhancing effects that increase temporal mean productivity in mixtures (7,13,16), an effect that should otherwise operate via population stability (28). The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few versus abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). This framework was tested with data from experimental tree communities from the BEF-China experiment (16,39), that span a long gradient of planted tree species richness with mixtures of up to 24 different tree species using piecewise structural equation models (SEMs) (40).

Fig. S9 Direct and indirect effects of tree species richness, drought-tolerance diversity and
CWMs of drought-tolerance traits on community stability. The structural equation model (SEM) tests the direct effects of tree species richness, functional diversity of stomatal control (FD stomatal control) and functional diversity of resistance-acquisition (FD resistanceacquisition) as well as their indirect effects mediated via asynchrony on community stability in the absence of population stability. Effects of CWM traits are explored through testing the effect of the CWM of stomatal control (CWM stomatal control) and the CWM of resistance-acquisition (CWM resistance-acquisition) on community stability. The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few vs abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). Functional diversity was calculated as abundance-weighted functional dispersion. The SEM fit the data well (Fisher's C=9.7, P =0.28, d.f.=8, n=218 plots). Data is based on a long, experimental species richness gradient with mixtures of 2, 4, 8, 16 and 24 tree species. Examined variables are shown as boxes and relationships as directional arrows with significant positive effects in blue, significant negative effects in red and non-significant paths in dotted grey based on a hypothesis driven SEM framework (Fig. S8). Standardized (significant) path-coefficients are shown next to each path with asterisks indicating significance (* P<0.05, ** P<0.01, *** P<0.001), path-width is scaled by coefficient size. Significant partial correlations are shown through grey, bi-directional arrows. Species richness was log2 transformed while asynchrony and community stability were square-root transformed. The variation in asynchrony and community stability explained by fixed (left, marginal R 2 ) and fixed together with random model effects (right, conditional R 2 ) is shown in the corresponding boxes. Note that the unexpected direct negative effect of functional diversity of drought tolerance on stability can be attributed to its positive effect on the temporal standard deviation (Fig. 5, marginally significant, P=0.06). This could be because the likelihood of including highly drought sensitive species increases with increasing tree species richness. Increased mortality as a result of including these sensitive species may cause increased temporal variation in community productivity at high species richness, as observed here (see Fig. 5).

Fig. S10 Direct and indirect effects of tree species richness, CWMs of drought-tolerance traits and drought-tolerance diversity on community stability. The structural equation model
(SEM) tests the direct effects of tree species richness, the CWM of stomatal control (CWM stomatal control) and the CWM of resistance-acquisition (CWM resistance-acquisition) as well as their indirect effects mediated via population stability on community stability in the absence of asynchrony. Effects of drought-tolerance diversity are explored through testing the effect of stomatal control (FD stomatal control) and functional diversity of resistance-acquisition (FD resistance-acquisition) on community stability. The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few vs abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). Functional diversity was calculated as abundance-weighted functional dispersion. The SEM fit the data well (Fisher's C=5.2, P =0.73, d.f.=8, n=218 plots). Data is based on a long, experimental species richness gradient with mixtures of 2, 4, 8, 16 and 24 tree species. Examined variables are shown as boxes and relationships as directional arrows with significant positive effects in blue, significant negative effects in red and non-significant paths in dotted grey based on a hypothesis driven SEM framework (Fig. S8). Standardized (significant) path-coefficients are shown next to each path with asterisks indicating significance (* P<0.05, ** P<0.01, *** P<0.001), path-width is scaled by coefficient size. Significant partial correlations are shown through grey, bi-directional arrows. Species richness was log2 transformed while population stability and community stability were square-root transformed. The variation in population stability and community stability explained by fixed (left, marginal R 2 ) and fixed together with random model effects (right, conditional R 2 ) is shown in the corresponding boxes.

Fig. S11. Hypotheses driven framework for the partitioning of the direct and indirect effects of species richness, drought-tolerance diversity and CWMs of drought-tolerance traits into overyielding and variance buffering effects of diversity.
Directional arrows represent hypothesized pathways. The framework separates the hypothesized effects of diversity on community stability (Fig. S8) into effects on the temporal mean (µAWP) and the temporal standard deviation of productivity (σAWP). Increases in µAWP would enhance community stability through overyieldinga higher productivity in mixtures vs monoculturesand decreases in σAWP would enhance community stability through buffered variations in productivity. All drivers hypothesized to influence community stability (Fig. S8)species richness, functional diversity of stomatal control (FD stomatal control), functional diversity of resistance-acquisition (FD resistance-acquisition), the CWM of stomatal control (CWM stomatal control), the CWM of resistance-acquisition (CWM resistance-acquisition) and asynchronyare partitioned here into their effects on µAWP and σAWP. Population stability was not included in this framework, as it did not respond to diversity nor CWM traits (Fig. 4). The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few versus abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). This framework was tested with data from experimental tree communities from the BEF-China experiment (16,39), that span a long gradient of planted tree species richness with mixtures of up to 24 different tree species using piecewise structural equation models (SEMs) (40). (72) captures the annual standardized water balance (precipitation minus potential evapotranspiration) and its variation during our study period. Negative values indicate climatic water deficits (brown coloured points) and positive values a water surplus (blue coloured points). Years with values above 1 or below -1 can be considered exceptionally wet and dry, respectively. The SPEI index was calculated based on monthly resolved precipitation and potential evapotranspiration data (CRU TS 4.04) (73) and is expressed here as annual water balance between two tree census intervals (SPEI12, water balance from September to October of the preceding year). The SPEI index was calculated with the SPEI package (74) in R using climate data from 1901-2019 as a reference period. (B) The residuals of this regression represent the annual variation in productivity without a directional stand development trend. The standard deviation of these residuals gives the detrended standard deviation that was subsequently used to calculate detrended community stability as temporal mean productivity (µAWP) divided by its detrended standard deviation (σAWP).

Fig. S14. Direct and indirect effects of tree species richness, drought-tolerance diversity
and CWMs of drought-tolerance traits on community stability. The structural equation model (SEM) is fit to data including monocultures and tests the direct effects of tree species richness as well as its indirect effects mediated via asynchrony and population stability on community stability (H1). Effects of functional diversity are explored through testing the effect of functional diversity of stomatal control (FD stomatal control) and functional diversity of resistance-acquisition (FD resistance-acquisition) as well as their indirect effects mediated via asynchrony on community stability (H2). Effects of CWM traits are explored through testing the effect of the CWM of stomatal control (CWM stomatal control) and the CWM of resistanceacquisition (CWM resistance-acquisition) as well as their indirect effects mediated via population stability on community stability (H3). The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few versus abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). Functional diversity was calculated as abundance-weighted functional dispersion. The SEM fit the data well (Fisher's C=11.6, P=0.48, d.f.=12, n=375 plots). Data is based on a long, experimental species richness gradient ranging from monocultures to mixtures of 24 tree species. Examined variables are shown as boxes and relationships as directional arrows with significant positive effects in blue, significant negative effects in red and non-significant paths in dotted grey based on a hypothesis driven SEM framework (Fig. S8). Standardized (significant) path-coefficients are shown next to each path with asterisks indicating significance (* P<0.05, ** P<0.01, *** P<0.001), path-width is scaled by coefficient size. Significant partial correlations are shown through grey, bidirectional arrows. Species richness was log2 transformed while asynchrony, population stability and community stability were square-root transformed. The variation in asynchrony and community stability explained by fixed (left, marginal R 2 ) and fixed together with random model effects (right, conditional R 2 ) is shown in the corresponding boxes. (Fig.  S8) into effects on the temporal mean (µAWP) and the temporal standard deviation of productivity (σAWP). Increases in µAWP enhance community stability through overyieldinga higher productivity in mixtures vs monocultures -and decreases in σAWP enhance community stability through buffered variations in productivity. All drivers hypothesized to influence community stabilityspecies richness, functional diversity of stomatal control (FD stomatal control), functional diversity of resistance-acquisition (FD resistance-acquisition), the CWM of stomatal control (CWM stomatal control), the CWM of resistance-acquisition (CWM resistanceacquisition) and asynchrony -were tested for their effects on µAWP and σAWP. Only significant pathways (P<0.05) are shown here to avoid overplotting (see Fig. S11 for the full model). Population stability was not included in this analysis, as it did not respond to diversity nor CWM traits (Fig. 4). The sketches schematically illustrate the trait gradients: water-spending vs watersaving stomatal control (few versus abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). The SEM fit the data well (Fisher's C=10.7, P=0.22, d.f.=8, n=375 plots). Data is based on a long, experimental species richness gradient ranging from monocultures to mixtures of 24 tree species. Examined variables are shown as boxes and relationships as directional arrows with significant positive effects in blue, significant negative effects in red and non-significant paths in dotted grey. Standardized (significant) pathcoefficients are shown next to each path with asterisks indicating significance (* P<0.05, ** P<0.01, *** P<0.001), path-width is scaled by coefficient size. Significant partial correlations are shown through grey, bi-directional arrows. Species richness was log2 transformed while asynchrony, µAWP and σAWP were square-root transformed. The variation in asynchrony, µAWP and σAWP explained by fixed (left, marginal R 2 ) and fixed together with random model effects (right, conditional R 2 ) is shown in the corresponding boxes.

Fig. S16. Direct and indirect effects of tree species richness, drought-tolerance diversity
and CWMs of drought-tolerance traits on community stability. The structural equation model (SEM) tests the direct effects of tree species richness as well as its indirect effects mediated via asynchrony and population stability on community stability (H1) and is fit without square-root transformed asynchrony, population stability and community stability but with an exponential variance structure (varExp) (67) for log2 species richness. Effects of functional diversity are explored through testing the effect of functional diversity of stomatal control (FD stomatal control) and functional diversity of resistance-acquisition (FD resistance-acquisition) as well as their indirect effects mediated via asynchrony on community stability (H2). Effects of CWM traits are explored through testing the effect of the CWM of stomatal control (CWM stomatal control) and the CWM of resistance-acquisition (CWM resistance-acquisition) as well as their indirect effects mediated via population stability on community stability (H3). The sketches schematically illustrate the trait gradients: water-spending vs water-saving stomatal control (few vs abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). Functional diversity was calculated as abundance-weighted functional dispersion. The SEM fit to the data was only marginally significant (Fisher's C=9.2, P=0.51, d.f.=10, n=218 plots). Data is based on a long, experimental species richness gradient with mixtures of 2, 4, 8, 16 and 24 tree species. Examined variables are shown as boxes and relationships as directional arrows with significant positive effects in blue, significant negative effects in red and nonsignificant paths in dotted grey based on a hypothesis driven SEM framework (Fig. S8). Standardized (significant) path-coefficients are shown next to each path with asterisks indicating significance (* P<0.05, ** P<0.01, *** P<0.001), path-width is scaled by coefficient size. Significant partial correlations are shown through grey, bi-directional arrows. Species richness was log2 transformed. The variation in asynchrony and community stability explained by fixed (left, marginal R 2 ) and fixed together with random model effects (right, conditional R 2 ) is shown in the corresponding boxes.  (Fig. S8) into effects on the temporal mean (µAWP) and the temporal standard deviation of productivity (σAWP) and is fit without square-root transformed asynchrony, µAWP and σAWP but with an exponential variance structure (varExp) (67) for log2 species richness. Increases in µAWP enhance community stability through overyieldinga higher productivity in mixtures vs monocultures -and decreases in σAWP enhance community stability through buffered variations in productivity. All drivers hypothesized to influence community stabilityspecies richness, functional diversity of stomatal control (FD stomatal control), functional diversity of resistance-acquisition (FD resistance-acquisition), the CWM of stomatal control (CWM stomatal control), the CWM of resistance-acquisition (CWM resistanceacquisition) and asynchronywere tested for their effects on µAWP and σAWP. Only significant pathways (P<0.05) are shown here to avoid overplotting (see Fig. S11 for the full model). Population stability was not included in this analysis, as it did not respond to diversity nor CWM traits (Fig. 4). The sketches schematically illustrate the trait gradients: water-spending vs watersaving stomatal control (few versus abundant stomata) and resistant vs acquisitive (high versus low cavitation resistance). The SEM fit to the data was only marginally significant (Fisher's C=14.1, P=0.079, d.f.=8, n=218 plots). Data is based on a long, experimental species richness gradient with mixtures of 2, 4, 8, 16 and 24 tree species. Examined variables are shown as boxes and relationships as directional arrows with significant positive effects in blue, significant negative effects in red and non-significant paths in dotted grey. Standardized (significant) pathcoefficients are shown next to each path with asterisks indicating significance (* P<0.05, ** P<0.01, *** P<0.001), path-width is scaled by coefficient size. Significant partial correlations are shown through grey, bi-directional arrows. Species richness was log2 transformed. The variation in asynchrony, µAWP and σAWP explained by fixed (left, marginal R 2 ) and fixed together with random model effects (right, conditional R 2 ) is shown in the corresponding boxes.  (33,34). See these studies for a detailed explication of the selected traits and Fig. S1 for additional information. In brief, we quantified stomatal control as species-specific stomatal sensitivity to water shortage via modelled curves of stomatal conductance (gs) under increasing water pressure deficits (VPD) (34). In these curves, VPDMAXFIT is the VPD at maximum stomatal conductance of a species and represents the threshold at which a species starts to limit stomatal conductance due to increasing VPD (34). VPDPOI is the VPD at the second point of inflection of this modelled curve of stomatal conductance, that is the point when the slope of the gs and VPD curve turns from negative to positive, and thus, can be taken as a measure of how quickly stomata close at high VPD, or in other words, of stomatal sensitivity. We consider both as physiological traits indicative of species-specific stomatal-control strategies: water spenders downregulate stomatal conductance only at high VPD while water savers close their stomata and down regulate stomatal conductance fast during increasing water shortage (Fig. S1). We have also included gsmax in the analysis (called CONMAXFIT, which is the modelled maximum conductance, see Kröber & Bruelheide (34)). However, in the BEF-China experiment, tree species' CONMAXFIT was uncorrelated to stomata size and stomata density but positively related to leaf area (i.e. to the leaf economics spectrum, see Fig. 4a in Kröber & Bruelheide (34)). For that reason, CONMAXFIT (as well as observed gsmax values) does not align well with our two PCA axes in Fig. S1. In contrast, VPDMAXFIT, which is the corresponding x value of the modelled gs curve, perfectly aligns with PC1, which is the water spender -water saver axis. Morphological stomata traits are stomata density and stomata index, the product of stomata size and stomata density. Note: Significant fixed effects (P<0.05) printed in bold. Drought-tolerance diversity was quantified as functional diversity of stomatal control (FD stomatal control) and functional diversity of resistance-acquisition strategies (FD resistance-acquisition) and CWMs of droughttolerance traits as CWM of stomatal control (CWM stomatal control) and of resistanceacquisition traits (CWM resistance-acquisition). Data is based on a planted diversity gradient ranging from monocultures up to mixtures of 24 tree species. Tree species richness was log2 transformed in all models. ddf are the denominator degrees of freedom; t the ratio between the estimate and its standard error; P-value from a t-distribution; n the number of plots; marginal R 2 values (R 2 m) show the variance explained by fixed effects and conditional (R 2 c) values the variance explained by fixed and random effects (71). Note: Shown are species and family names, species identity codes, leaf habit (E, evergreen; D, Deciduous) and the site at which the species were planted. For more details on the characteristics and taxonomy of the tree species see ref. (16) and for the experimental design refs. (16,39).