Drivers of productivity and its temporal stability in a tropical tree diversity experiment

There is increasing evidence that mixed‐species forests can provide multiple ecosystem services at a higher level than their monospecific counterparts. However, most studies concerning tree diversity and ecosystem functioning relationships use data from forest inventories (under noncontrolled conditions) or from very young plantation experiments. Here, we investigated temporal dynamics of diversity–productivity relationships and diversity–stability relationships in the oldest tropical tree diversity experiment. Sardinilla was established in Panama in 2001, with 22 plots that form a gradient in native tree species richness of one‐, two‐, three‐ and five‐species communities. Using annual data describing tree diameters and heights, we calculated basal area increment as the proxy of tree productivity. We combined tree neighbourhood‐ and community‐level analyses and tested the effects of both species diversity and structural diversity on productivity and its temporal stability. General patterns were consistent across both scales indicating that tree–tree interactions in neighbourhoods drive observed diversity effects. From 2006 to 2016, mean overyielding (higher productivity in mixtures than in monocultures) was 25%–30% in two‐ and three‐species mixtures and 50% in five‐species stands. Tree neighbourhood diversity enhanced community productivity but the effect of species diversity was stronger and increased over time, whereas the effect of structural diversity declined. Temporal stability of community productivity increased with species diversity via two principle mechanisms: asynchronous responses of species to environmental variability and overyielding. Overyielding in mixtures was highest during a strong El Niño‐related drought. Overall, positive diversity–productivity and diversity–stability relationships predominated, with the highest productivity and stability at the highest levels of diversity. These results provide new insights into mixing effects in diverse, tropical plantations and highlight the importance of analyses of temporal dynamics for our understanding of the complex relationships between diversity, productivity and stability. Under climate change, mixed‐species forests may provide both high levels and high stability of production.


| INTRODUC TI ON
Forest restoration has been identified as the most important natural solution for climate change mitigation (Griscom et al., 2017). For example, the Bonn Challenge advocates for +150 Mil ha of restored forests by 2020. Tropical and subtropical regions, where ongoing deforestation and forest degradation have left extensive areas of degraded land, provide a unique opportunity for restoring productive forests (Bauhus, van der Meer, & Kanninen, 2010). Planted forests can provide many of the ecosystem functions and services of natural tropical forests, albeit some at a lower level (Bauhus et al., 2010;Pawson et al., 2013), while making a substantial contribution to satisfy the increasing demand for global roundwood (Kanninen, 2010).
Planted forests are, however, in most cases still established as monocultures, often with non-native tree species , despite the fact that mixed-species forests (either planted or from natural regeneration) are considered important for adaptation of forests in the face of global change (Messier, Puettmann, & Coates, 2013;Pawson et al., 2013). Tree species mixtures can provide multiple ecosystem services at higher levels than their monospecific counterparts, although this may not be the case for all ecosystem services (Gamfeldt et al., 2013;van der Plas et al., 2016). The strongest evidence for such positive mixing effects exists for productivity and C sequestration (e.g. Forrester & Bauhus, 2016;Paquette & Messier, 2011;Piotto, 2008). While there is some evidence regarding higher ecological stability and resilience of mixed-species stands in relation to specific stress and disturbance factors (Bauhus, Forrester, Gardiner, et al., 2017;Hutchison, Gravel, Guichard, & Potvin, 2018;Jactel et al., 2017), there are very few long-term analyses of stability of productivity in relation to tree diversity.
Understanding the mechanisms behind the relationship of biodiversity with ecosystem functioning (BEF) is crucial for designing and implementing diverse, resilient and productive planted forests.
Studies of BEF relationships in forests have employed various approaches over the last two decades, ranging from analysis of forest inventories to experimental plantations specifically designed to test BEF relationships Nock et al., 2017). While each approach has its specific strengths and drawbacks, experiments provide the strongest test of BEF effects by controlling for underlying environmental effects and directly comparing tree performance in monocultures and in mixtures (see . Given experiments are generally inventoried more frequently than forest plots, they offer unique opportunities to study temporal developments of diversity effects (Huang et al., 2018), as well as the buffering effects of diversity from disturbance and environmental variation (Isbell et al., 2015).
To date, few such analyses have been conducted, as most forest BEF experiments are still young. One exception is 'Sardinilla' in Panama, the oldest BEF experiment in the tropics (Scherer-Lorenzen et al., 2005), which was used here to analyse diversity-stability relationships (DSRs) and the temporal development of diversity-productivity relationships (DPRs).
Net overyielding occurs when productivity in mixtures is higher than in monocultures. This can be attributed to the combined effect of competitive reduction (+) and facilitation (+) versus competition (−) and is also referred to as 'complementarity effect' (Forrester & Pretzsch, 2015). However, enhanced mixture productivity might also result from the dominance of one or few species caused by selection or mass ratio effects (Fotis et al., 2018;Grime, 1998). Indeed treetree interactions that scale up to community-level responses can be positive or negative, depending on species' assemblage and environmental influences (Forrester & Bauhus, 2016;Forrester & Pretzsch, 2015) with tree size and competition by neighbouring trees strongly influencing diversity effects on single-tree productivity (Dănescu, Albrecht, & Bauhus, 2016;Fichtner et al., 2018). Hence, to develop resilient plantations, it is crucial to clarify the context dependency of DPRs in forest ecosystems: under which climatic conditions, during which stage(s) of stand development and at what levels of diversity forest managers can expect beneficial effects of mixtures on productivity (Forrester & Bauhus, 2016).
In the face of future climatic stress and disturbances, it will become increasingly important to design plantations not only to increase productivity but also to stabilize it. The effects of tree species diversity on the resistance to drought, wind, fire, pests and pathogens appear equivocal and, in most cases, except for herbivorous insects, the evidence base is weak (Bauhus, Forrester, Gardiner, et al., 2017). Even less is known about the effects of diversity on the temporal stability (Lehman & Tilman, 2000;Tilman, 1999) of community productivity in forests, that is, the fluctuation of productivity around its long-term mean. In grassland ecosystems, there is abundant evidence that interannual fluctuations of community-level productivity are smaller in more diverse compared to less diverse communities, resulting in a net positive DSR (Hautier et al., 2014;Isbell et al., 2015). The few studies that have analysed temporal stability in temperate and boreal forest ecosystems support the hypothesis that diversity can stabilize community-level productivity (Aussenac, Bergeron, Gravel, & Drobyshev, 2019;del Río et al., 2017;Jucker, Bouriaud, Avacaritei, & Coomes, 2014;Morin, Fahse, de Mazancourt, Scherer-Lorenzen, & Bugmann, 2014 2018) but no detailed analysis of the underlying drivers of this phenomenon exist.
The overall aim of this study was to test not only whether diversity increases productivity and its temporal stability in mixed stands but also to identify whether this might be driven by species diversity or structural diversity, whether stability was more influenced by overyielding or asynchronous growth of tree species (Jucker et al., 2014), and how these influences change with stand development. We expect that the strength of positive diversity effects and consequently overyielding in mixed species stands, which were reported for the first half of the experiment's lifespan Sapijanskas, Potvin, & Loreau, 2013), increases with stand development. Moreover, the period of development of the 'Sardinilla' plantation has been characterized by contrasting climatic conditions, including an exceptionally wet and an exceptionally dry La Niña and El Niño period, respectively (Detto, Wright, Calderón, & Muller-Landau, 2018;Hutchison et al., 2018). Here, we used this climatic variation to examine whether DPRs, as previously hypothesized (Forrester & Bauhus, 2016), change along a gradient of climate-induced water variability. Hutchison et al. (2018) showed that tree mortality in the monocultures of Sardinilla was modulated by extreme climatic events while species mixing buffered against this effect. We expect that lower climatic sensitivity of mixtures is driven by species asynchrony, that is, the fluctuating responses of species to contrasting climatic conditions (Jucker et al., 2014) and overyielding. We test whether these mechanisms translate into an overall positive DSR, expressed here as one aspect of stability, the single and intuitive metric 'temporal stability'.
While tree community-level analyses are common in forestry studies aiming to produce results for a management-relevant scale, it is increasingly recognized that community responses in mixed stands are driven by tree-tree interactions at the neighbourhood level (Dănescu et al., 2016;Fichtner et al., 2018;Potvin & Dutilleul, 2009). Importantly, neighbourhood analyses allow to accurately describe variability in stand density, mortality and stand structure (Forrester & Bauhus, 2016;Forrester & Pretzsch, 2015). Here, we employ a combined analysis of DPRs at both the community and tree neighbourhood scale to provide insight into the complex interplay of complementarity effects during stand development.
Finally, most studies on BEF relationships in forests simply use species richness or diversity as the measure of tree diversity (Forrester & Bauhus, 2016). However, structural diversity is increasingly recognized as another key attribute influencing productivity and stress tolerance of trees (Dănescu et al., 2016;Pretzsch, Schütze, & Biber, 2018). Applying this perspective, trees of different sizes occupy distinct niches and could behave, at least to a certain degree, like functionally different species (Dănescu et al., 2016). Recently, an indirect positive effect of species diversity on productivity via changes in size structure (Zhang, Chen, & Coomes, 2015) and also a direct positive, nonmediated, effect of structural diversity (Dănescu et al., 2016) have been described.
In contrast, structural diversity had negative effects on productivity in monospecific, clonal eucalypt stands established through staggered planting to create structural diversity (Binkley, Stape, Bauerle, & Ryan, 2010;Ryan et al., 2010). In tree mixtures with highly divergent tree growth rates, effects of both components of diversity should be tested. This has not been done so far under controlled conditions. We are also not aware of any study that examined effects of structural diversity on temporal stability of tree growth.
Thus, this paper addresses the following hypotheses: 1. Overyielding increases with stand development and is highest in the most diverse tree neighbourhoods and stands (plots); 2. Both tree species diversity and structural diversity increase productivity and its temporal stability during a period with contrasting climatic features; 3. Temporal stability of productivity in mixtures is driven by species asynchrony and overyielding.

| Study design
The study was conducted in an experimental planted forest, and mixtures of all six-species (N = 6) with an average of 233 individuals per plot (equalling 1,150 trees/ha; Potvin & Dutilleul, 2009;Scherer-Lorenzen et al., 2005). To ensure trait divergence in each mixture, species were allocated based on their relative growth rates. In each three-species mixture, one fast growing pioneer species, either Luehea seemannii (LS) or Cordia alliodora (CA), one light-intermediate species, either Anacardium excelsum (AE) or Hura crepitans (HC) and one slow-growing and shade-tolerant species, either Tabebuia rosea (TR) or Cedrela odorata (CO) were chosen randomly based on their relative growth rates in nearby natural forests (Scherer-Lorenzen et al., 2005). One species, CA, suffered mortality rates >90% in the first years after planting likely a result of site properties, possibly due to the compacted and undrained soil Sapijanskas et al., 2013) or root herbivory by beetle larvae (Healy, Gotelli, & Potvin, 2008), but not a diversity effect. We, therefore, excluded it from our analysis (for details see Appendix S2). In this study, we consequently refer to the 'realized' species richness levels in Sardinilla, which comprise monocultures (N = 10), two-(N = 3), three-(N = 3) and five-species mixtures (N = 6).

| Performance proxies of tree growth
At the time of analysis, trees in the plantation were 16 years old, with the tallest trees over 25 m, old enough, therefore, to examine diversity-stability relationships (DSRs) and the temporal development of diversity-productivity relationships (DPRs). We used two key response variables, diameter and height growth, as proxies of tree performance. Diameter and height were measured annually from 2002 to 2017 for all trees in the plantation at the end of each growing season (December-January). Diameter was measured at breast height (1.3 m) for trees with a total height >2 m for each stem (i.e. in case of multistemmed trees). We chose 2006, when 85% of all trees had reached a height of >2 m, as start year of our analysis to ensure a complete and consistent data set. Our inventory data set thus comprises a complete, spatially explicit, inventory of all trees measured annually from January 2006 to January 2017. Basal area increment (BAI) between successive years was used as a proxy for productivity and was calculated as: where dbh is diameter at breast height, j is an index for the n stems of each tree (i.e. for multistemmed trees) and t is an index for the year of survey.
We corrected negative BAI values (e.g. caused by measurement errors or wind damage) in order to avoid biased model estimates.
As described in the previous work, negative increments represent a common challenge in experiments utilizing inventory data sets with a high temporal resolution (Fichtner et al., 2018). The applied procedure included predicting basal area at dbh with additionally available basal diameter data via species-specific allometric relationships and the exclusion of wind-damaged trees, for example, as after the tropical storm 'Otto' in 2016 (see Appendix S2 for details).
To avoid edge effects, we excluded the outer border row of trees of each plot for calculation of response variables. However, to best reflect the actual growing conditions of trees, we calculated all predictor variables (e.g. diversity indices) using data of all trees in the plantation including border trees.

| Community-and tree neighbourhoodlevel analyses
First, we analysed DPRs and DSRs at the community level in relation to community tree species richness, the most common diversity predictor in BEF research (Forrester & Bauhus, 2016;Jucker et al., 2014). Secondly, we used a wider set of diversity measures to assess the impact of species and structural diversity on productivity and its stability at the community level. We hypothesize that DPRs at the community level (here the plot) may be driven by tree-tree interactions in neighbourhoods. In a third step, we explored the underlying mechanisms of community-level DPRs through modelling the influence of the same set of candidate indices on growth of individual trees at the neighbourhood level (Figure 1).
At the community level we calculated annual productivity (BAI plot , m 2 ha −1 year −1 ) as the sum of BAI of all trees per plot or species that were alive in a particular year (N = 2,596). In contrast, at the neighbourhood level, we analysed single-tree growth (BAI tree , cm 2 /year) of all trees that were alive at the end of the observation , F I G U R E 1 Design of the community and tree neighbourhood level analyses. On the right, the plots in the Sardinilla plantation are shown, coloured according to their species richness level (N = 22). The black points represent the position of individual, living focal trees (2006-2016) whose productivity was modelled in response to their immediate neighbours. Grey points show all trees planted in 2001 (including dead individuals and border trees). On the left, the design of the tree neighbourhood analysis is illustrated. The central black tree represents the focal tree with its immediate neighbours, up to a maximum of eight living trees. Community-and neighbourhood-level productivity and predictor variables (e.g. diversity indices) were calculated for each year (2006-2016) based on annually resolved values of tree basal areas and heights period, hereafter called 'focal trees' (N = 2,159, Figure 1). Annually resolved values of species diversity and structural diversity were calculated for the whole community (each plot) and for each tree's neighbourhood.
The neighbourhood of focal trees comprised of its immediate neighbours, that is, all living trees within a radius of 5 m. This resulted in a maximum of eight neighbours, considering the fixed planting design of the plantation (3 × 3 m, Figure 1). Compared to Potvin and Dutilleul (2009)

| Measures of species and structural diversity
To improve our understanding of the processes that drive DPRs and DSRs in forests, we used a wide set of species diversity and structural diversity indices, because previous studies have shown that the choice of indices can strongly influence the outcome of analyses (Dănescu et al., 2016;Schnabel, Donoso, & Winter, 2017).
In other studies, species diversity has been considered a component of forest structure (Dănescu et al., 2016). Here, we quantified tree species diversity using three conventional indices: (a) species richness (i.e. the number of tree species), (b) the Shannon diversity index (Shannon, 1948) using relative basal area to quantify species proportions and (c) evenness, calculated as Shannon index divided by its theoretical maximum (Table 1).
To quantify diameter and height diversity we calculated widely used metrics of forest structural diversity: (a) standard deviation (sd), (b) coefficient of variation (CV) and (c) Gini coefficient (GC; Gini, 1912; Table 1). Higher index values reflect higher structural diversity for all indices (see Lexerød & Eid, 2006 for a detailed index comparison). Hereafter, we refer to these indices as measures of 'species diversity' and 'structural diversity'. We acknowledge that these indices reflect only a small subset of all aspects of forest structure, namely species diversity and the variation in tree diameters and heights and measure different aspects of diversity (e.g. variation, diversity and inequality).

| Community-level overyielding
To quantify overyielding or underyielding of mixtures, we where BAI mix is the BAI plot of all species in the mixture, BAI mono is the BAI plot of the respective monoculture of species i and m i is the proportion of species i in mixture corresponding to its initial planting density (N trees/ha). (2) TA B L E 1 Summary of the species diversity and structural diversity indices used in this study. Data are for the community level (N = 22 plots × 11 years) and the tree neighbourhood level (N = 2,159 trees × 11 years) Component Index Acronym and equation

Range Mean Range
Species diversity  Gini coefficient

Species richness
18 0.08-0.28 0.16 0.00-0.65 Note: All variables and indices were calculated for annually resolved values per plot (all living trees including border and snapped trees) and per tree neighbourhood (alive, immediate neighbours). ba is the basal area measured at 1.3 m (cm 2 ) and h is the height (cm) of tree i, P the proportion of ba for species i, x ba and x height are the mean tree ba and height, j is a tree's rank is ascending order from 1 to n, ba j is the basal area and h j is the height of the tree with rank j. Mean values show the temporal mean and range covers the respective minimum and maximum values for the observation period 2006-2016.
As already mentioned, the failure of CA led to a decrease in tree density in some mixtures, which resulted in more growing space for the remaining individuals as compared to the denser monocultures. To avoid an overestimation of diversity effects, we considered the realized richness by excluding CA individuals when calculating mixing proportions (i.e. m i equal to 0.50 and not 0.33 in the former three-and now two-species mixtures). As m i determines the monoculture productivity, here presented RP values should be considered as conservative estimates of overyielding.
RP of a given species i was calculated for each year as: where BAI i mix is the total BAI per plot for species i in the mixture and BAI i mono is the total BAI per plot for species i in the monoculture.
Comparing the RP among species allowed us to disentangle the possibly contrasting mixing-effects of single species that translate into a net complementarity effect of the entire community. To calculate comparable productivity estimates per species, we used the mixing proportion m i as BAI plot species /m i . A one-sample t test was used to compare the temporal mean RP community and RP species against 0 (the expected monoculture yield).

| Growth stability
We tested for DSRs at the community level by analysing the temporal stability of tree productivity, hereafter 'stability', following Jucker We calculated species asynchrony at the community level using the species synchrony measure φ (Loreau & de Mazancourt, 2008) as 1 − φ: where BAI species i is the standard deviation of productivity of species i in a community of n species (Hautier et al., 2014;Jucker et al., 2014).

| Modelling the drivers of diversity-productivity relationships
To understand the underlying drivers of DPRs, we modelled tree productivity at the community and tree neighbourhood level in relation to species and structural diversity (Table 1), while also accounting for other factors relevant for tree growth ( Table 2). The latter included tree mortality at the community level and tree size and competition at the neighbourhood level (Dănescu et al., 2016;Fichtner et al., 2018;Hutchison et al., 2018;).

| Community-level growth models
At the community level, we accounted for the effect of tree mortality on productivity (BAI plot ) by considering relative mortality in plots (relM; Table 2). To model temporal trends in DPRs we incorporated year, its squared form and tested interactions between year and each candidate diversity index (Table 2). We did not include a stand density proxy like plot basal area due to the high correlation with year to avoid collinearity. Alternative model runs with this proxy produced similar results. Annually resolved values (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) were used for all response and predictor variables (Tables 1 and 2).
TA B L E 2 Overview of response variables and the final set of nondiversity growth predictors used in the community-and tree neighbourhood-level models

| Neighbourhood-level growth models
At the tree neighbourhood level, we built species-specific models to address differences among species. In addition to diversity indices, we included the following growth relevant factors (see Appendix S4 for details): (a) relative mortality (relM) of immediate neighbours to account for neighbourhood mortality, (b) a focal tree's logtransformed ba to account for the effect of tree size and (c) ba of the immediate neighbours j larger than the focal tree i (BAL) to account for competition calculated as ∑ j≠i ba j (Table 2). We modelled the mean annual BAI tree of individual trees during

| Diversity-productivity relationships
We used linear mixed-effects models (LMMs) to gain insight into the temporal development of diversity-productivity relationships while accounting for the inherently correlated errors in our data. Factors identifying replication and spatial structure (tree species composition at the community and additional plots, subplots and focal tree identity at the tree neighbourhood level) were modelled as nested random effects. We followed the same model selection procedure at the community and neighbourhood level to ensure comparability of results. To test whether species diversity or structural diversity or a combination of both affect tree growth, we developed a null model of tree productivity without any diversity index and a series of models, each incorporating one diversity component. Model development involved the stages proposed by Zuur, Ieno, Walker, Saveliev, and Smith (2009): 1. Specifying a null model (i.e. excluding diversity variables) with a beyond-optimal selection of fixed effects.
2. Optimizing the random structure (random effects, temporal autocorrelation, variance structure) in the presence of the beyondoptimal model specification.
3. The optimal null model structure was chosen via removing all nonsignificant fixed effects (see Table 2 for the final set) and was kept fixed in the subsequent analysis.
4. Testing diversity indices: we included the species and structural diversity indices (Table 1)  Due to the non-normal distribution of the response we applied a fourth root transformation at the community level and a Box-Cox transformation (Box & Cox, 1964) at the tree neighbourhood level. We only included predictors into the models that were not collinear (Spearman's rho < 0.6 and a variance inflation factor for mixed models <5, which is a conservative choice (Dormann et al., 2013).

| Diversity-stability relationships
We used linear regression to test our hypothesis that species and structural diversity stabilize productivity in mixtures via regressing community stability against diversity indices (Table 1). As described in the LMM framework above, Akaike weights were used to determine the best candidate diversity indices. To determine whether diversity effects on stability resulted from increased µ BAI or decreased σ BAI , we regressed both against the best-performing index.
Finally, we tested for the effect of species asynchrony by regressing it against stability. If residuals were not normally distributed, we applied a log-transformation to the response variable. Simple linear regression was used to analyse stability, quantified as temporal means without repeated measurements.
All analysis were performed in R (version 3.5.0) using the packages nlme (Pinheiro, Bates, DebRoy, Sarkar, & R Core Team, 2018), piecewiseSEM (Lefcheck, 2016), AICcmodavg (Mazerolle, 2019) and ggplot2 (Wickham, 2016) for graphics. Validity of model assumptions was tested via graphical tools (quantile-quantile, residual, autocorrelation and partial autocorrelation plots). LMMs were fit with the package nlme to allow for the specification of variance functions, to address heteroscedasticity and to model temporal autocorrelation (see Appendix S4 for technical details on LMM development, structure and model evaluation). Finally, we explored whether the overlap between the considered tree neighbourhoods ( Figure 1) might influence our results. We tested two alternatives: models based on a subset of data with strictly nonoverlapping neighbourhoods or with spatial instead of temporal autocorrelation, which yielded consistent results (see Appendix S7 for details).

| Community-level overyielding and stability
All mixtures were more productive than monocultures and this overyielding increased with species richness (Figure 2a). Across all years, mean overyielding was 25%-30% in the two-and three-species mixtures and nearly 50% in the five-species mixtures. Overyielding was only significant for the two-and five-species mixtures. Calculating overyielding based on the 'initially planted species richness', in contrast to the here used 'realized species richness', increased overyielding estimates to ~80%-90% in the two-and five-species mixtures, while overyielding in the three-species mixtures remained unchanged (no CA individuals were planted in these plots; Figure S6).
Overyielding differed among the five tree species: the fastest growing species, LS (pioneer) and AE (light-intermediate) showed highly significant overyielding, while the two shade-tolerant species (CO, TR) did not (Figure 2b). In CO, the response was highly variable but tended to be positive (p = .08, Figure 2b). HC (intermediate) was the only species with significant underyielding in mixtures (Figure 2b).
Overyielding significantly increased (p = .0082) over time in the five-species mixtures, while no significant trend was detected in two-and three-species mixtures (Figure 3). An apparent decline in overyielding in three-species mixtures from high values in the first years (2006)(2007)(2008) (Figure 3; Figure S7) was not significant. While we found overyielding for all richness levels and in most years between 2006 and 2016, it was lowest and even negative around 2010 ( Figure S7). This drop, during a particularly wet period, coincided with a peak in monoculture productivity but no consistent change in mixtures (Figure 4; Figure S3). Differences in species productivity F I G U R E 2 Temporal mean (2006-2016) relative productivity (RP) of mixtures compared to monocultures. RP was calculated based on community-level productivity (m 2 ha −1 year −1 ) with Equations (2) and (3). RP is presented per realized species richness levels (a) and individual species (b). Species are ordered according to their relative growth rates in natural forests, from fast (left) to slow (right). The five species are Luehea seemannii (LS), Anacardium excelsum (AE), Hura crepitans (HC), Tabebuia rosea (TR) and Cedrela odorata (CO). The zero line represents the monoculture yield. Significant differences between mixtures and monocultures are shown as stars with *p ≤ .05, **p ≤ .01, ***p ≤ .001, n.s., not significant  (Figure 5a). The stabilizing effect of diversity resulted mainly from a significant positive effect of richness on the temporal mean of productivity (µ BAI ), while there was no significant effect on the temporal variation in productivity (σ BAI ; Figure S8). Species asynchrony had a strong and highly significant positive effect on stability, consistent with a trend towards higher asynchrony at higher richness levels ( Figure 5b). µ BAI had a positive and σ BAI a negative relationship with asynchrony, albeit both relationships were not significant (results not shown). All candidate species diversity indices had a similar impact on stability, with the highest Akaike weight for species richness (w = 0.31), followed by the Shannon diversity (w = 0.24) and evenness index (w = 0.15; results not shown). We found no significant effect of structural diversity on stability and low Akaike weights for all structural diversity indices (w = 0.03 or below; results not shown).

| Drivers of community productivity
Species and structural diversity both had overall positive influences on productivity but displayed contrasting patterns over time.
While the positive effect of species diversity significantly increased (p = .0095), the opposite was true for structural diversity. For the latter, the strong positive effect in the first years (tree age >5) F I G U R E 4 Temporal development of community-level productivity (BAI plot , m 2 /ha) per year and species, in connection with four levels of species richness (1, 2, 3 and 5 species; see boxes above panels). Temporal patterns are described by local polynomial regressions. Points represent observed values of productivity per plot for 11 years (N = 10 × 11 for monocultures, N = 3 × 11 for two-and three-species mixtures and N = 6 × 11 for five-species mixtures). The productivity of species in mixtures was scaled to the monoculture yield with the mixing proportion m i as BAI plot species /m i . In the legend, species are ordered according to their relative growth rates in natural forests, from fastest (LS) to slowest (CO). AE, Anacardium excelsum; CO, Cedrela odorata; HC, Hura crepitans; LS, Luehea seemannii; TR, Tabebuia rosea F I G U R E 5 Temporal stability of community productivity as a function of species richness (a) and species asynchrony (b). Stability and species asynchrony were calculated with Equations (4) and (5) respectively. For (a) fitted values were back transformed from a log-scale to match the original values per plot (grey points). For (b) points represent values per plot and colours the respective species richness. The grey-shaded areas show a 95% confidence interval for the fitted models significantly declined over time (p = .0049; Figure 6; Table S3). Standlevel productivity significantly decreased over time and with higher mortality of trees (p ≤ .0001; Figure 6, Table S3). All species diversity indices showed similar and significantly positive effects, but species richness led to the most parsimonious model (Table 3 and Table S2).
The standard deviation of tree height led to the most parsimonious structural diversity index model (Table 3 and Table S2), but also the other diameter and height diversity indices that increased model performance had overall positive effects. Akaike weights and direct standardized effect sizes clearly supported the relative superiority of species diversity over structural diversity, but changes in productivity over time were stronger for structural diversity (Table 3; Table S3).

| Drivers of individual tree productivity at the neighbourhood level
All species-specific null models had a similar set of fixed effects.
Significant nondiversity growth predictors unrelated to diversity were (in decreasing order of their effect size): (a) size of the focal tree with a positive effect on productivity; (b) competition (expressed as BAL), which had a negative effect on productivity; and (c) neighbour mortality that increased productivity of focal trees (Table 3;   Tables S4-S8;  Indices of species and structural diversity were significantly related to individual tree productivity for three out of five species. Growth in the species HC and TR responded only to species diversity but not to structural diversity. To capture the partly contrasting effects of species and structural diversity on growth performance and to simultaneously focus on the most relevant effects, we report here only the highest ranking species and structural diversity index (see Table S2 for a detailed index ranking). All presented models explained high shares of the variation in tree productivity (marginal R 2 ≥ 0.74; Table 3). Finally, as interactions of the diversity component F I G U R E 6 Changes in species diversity (a) and structural diversity (b) effects over time. Diversity effects show the change in community-level productivity (BAI plot , m 2 /ha) at the start (2006)  Note: Only the most parsimonious (highest Akaike weight) models with significant diversity components are presented (see Appendices S4 and S6 for detailed selection criteria and index ranking). Akaike weights (w) show the relative support for a candidate diversity model (with 0 for a low and 1 for the highest relative likelihood). Spearman's rho shows the correlation between the species and structural diversity indices. Marginal R 2 values (R 2 m ) represent the variance explained by the fixed and conditional R 2 values (R 2 c ) the variance explained by fixed and random effects (Nakagawa, Schielzeth, & O'Hara, 2013) of the final model. For a detailed description of the diversity indices, see Table 1 and for other variable names, Table 2. with period and/or a tree's dominance class were significant (p ≤ .05) for all species, we only report on these higher-order model terms (Table 4).
Neighbourhood diversity effects on individual tree productivity strongly varied between the examined periods (Table 4). We found the overall strongest species diversity and structural diversity effects in the third period. In contrast, in the second period, measures of diversity displayed an overall low, nonsignificant influence on tree productivity. In addition, relative size of focal trees modulated species diversity but not structural diversity effects on single-tree productivity (Table 4). Productivity in focal trees of different species responded differently to species diversity, but overall positive effects prevailed (Table 4). Effects were consistently positive for the fast-growing species LS, the intermediate AE and the shade-tolerant CO that experienced the strongest effects in the first (LS) and third period (AE, CO) respectively. For these three species, positive species diversity effects on productivity increased with a focal tree's dominance (Table 4). Only the slow-growing species HC showed a consistently and increasingly negative relationship with species diversity (Table 4;   Table S6). Interestingly, species diversity effects on productivity of the shade-tolerant TR were comparably variable over time (Table 4).
Effects of structural diversity on individual species changed over time and became more contrasting (Table 4). In the first period, effects were consistently positive for all species except HC (Tables S4-S8), but this effect was significant only for the shadetolerant species CO (Table 4). In the third period, structural diversity had especially a positive effect on the intermediate species (AE), while the fast-growing species LS was negatively influenced. Finally, we found a moderate positive correlation among species and structural diversity for most species (LS, AE, TR, CO, Table 3).

| Temporal dynamics of species diversity-productivity relationships
Results of this study clearly support our hypothesis of a positive diversity-productivity relationship (DPR) that increases in strength over time. While earlier studies in Sardinilla concluded that twoand three-species mixtures were the most productive (e.g. Healy et al., 2008;, our analysis of a longer temporal record shows that overyielding increased with stand development in the five-species mixtures, which overall outperformed the less diverse mixtures. Thus, if only early years of the experiment were analysed, very different conclusions would be drawn regarding optimal richness levels for the productivity of mixed-species plantations, that is, few versus many species. Several processes could explain this finding: differences between species likely required time to develop into complementary TA B L E 4 Significant species diversity and structural diversity effects (p ≤ .05) on individual tree productivity at the tree neighbourhood level  ). The diversity component (species diversity or structural diversity) with the higher Akaike weight (w) for each species is printed in bold (Table 3). For linear mixed-effects model parameter estimates, see Tables S4-S8. Species are ordered according to their relative growth rates in natural forests, from fastest (LS) to slowest (CO). Abbreviations: AE, Anacardium excelsum; CO, Cedrela odorata; HC, Hura crepitans; LS, Luehea seemannii; n.s., not significant; TR, Tabebuia rosea.
interactions, for example to realize crown packing induced overyielding (Jucker, Bouriaud, & Coomes, 2015). Furthermore, species with similar growth rates were planted in direct adjacency in the five-species mixtures of Sardinilla, which likely caused an initial spatial dilution of complementary species characteristics (Sapijanskas et al., 2013;Scherer-Lorenzen et al., 2005). This would lead to the observed pattern of overyielding increasing with stand development as contrasting species started to interact. Finally, the fast-growing species LS likely benefitted from neighbourhood species diversity especially in the first years when it experienced no or little competition for light from slower growing competitors (Table 4; Figure 3). Community overyielding was consequently high in the three-species mixtures (one third of individuals are LS trees) in the first years ( Figure S7).
Current knowledge on the mechanisms driving the temporal dynamics of DPRs in forests is so far based nearly exclusively on results from two-species mixtures (Forrester & Bauhus, 2016). Only very recently, a similarly increasing DPR over time has been reported for the BEF-China and FORBIO tree diversity experiments (Huang et al., 2018;van de Peer, Verheyen, Ponette, Setiawan, & Muys, 2018) and for a replacement-series mixture experiment in Costa Rica (Ewel, Celis, & Schreeg, 2015), while most other forest BEF experiments are still too young for such temporal analysis. Our results add to the increasing evidence that species diversity is a key tool to increase forest productivity (Huang et al., 2018;Jactel et al., 2018) and caution against conclusions drawn from early stages of experiments and afforestation trials, which cannot take into account the temporal dynamics of DPRs.
The importance of long-term observations in diversity experiments is also supported by studies of commercial mixed-species plantations with native tree species in Costa Rica (Piotto, Craven, Montagnini, & Alice, 2010;Redondo-Brenes & Montagnini, 2006). Here, several species that grew well at a juvenile stage suffered high mortality with progressing stand development (Piotto et al., 2010).
To mechanistically explain our findings regarding overyielding, we modelled tree productivity as a function of species and structural diversity both at the community and tree neighbourhood level. In tree neighbourhoods, species and structural diversity effects on single-tree productivity varied (a) among species, (b) with progressing stand development, (c) according to contrasting climatic conditions and (d) with relative dominance of focal trees. The general patterns found for neighbourhoods, however, were consistent with those found at the community level. We conclude that tree-tree interactions in neighbourhoods are the principal drivers of diversity effects observed at the community level, confirming earlier results from Sardinilla (Potvin & Dutilleul, 2009) and other tree diversity experiments (Fichtner et al., 2018;van de Peer et al., 2018). In line with former results from Sardinilla (Potvin & Dutilleul, 2009;Sapijanskas et al., 2013), tree size was the strongest determinant of tree-level productivity, followed in our study by competition and mortality. Importantly, species diversity effects on productivity increased with dominance of individual trees for species with different shade-tolerances. This indicates that one important diversity effect is competitive reduction for fast growing tree species, while slower growing and shade-tolerant species may also benefit from species diversity if they are less suppressed by their neighbouring companion species. This is consistent with results from the BEF-China experiment, where small trees that experienced comparatively little competition benefitted most from neighbourhood species richness (Fichtner et al., 2018).
The combined contribution of tree size and competition by neighbours to diversity effects highlights the increasingly recognized role of local tree neighbourhoods for understanding DPRs in forest ecosystems (Forrester & Pretzsch, 2015;Stoll & Newbery, 2005).
Here, we calculated the RP of mixed-species systems (compared to monocultures) for the whole community and for individual species to take into account selection versus complementarity effects (Forrester & Bauhus, 2016;Forrester & Pretzsch, 2015). Sapijanskas, Paquette, Potvin, Kunert, and Loreau (2014) showed that tree species diversity enhanced community-level light capture and hence productivity. Importantly, complementary water use was found in mixtures in Sardinilla, caused by three distinct water uptake depths of participating species (Schwendenmann, Pendall, Sanchez-Bragado, Kunert, & Hölscher, 2015), which could cause competitive reduction through niche differentiation (spatial segregation) and possibly also facilitation due to hydraulic redistribution of water from deeper to shallower soil layers (Forrester, 2017

| Structural diversity effects on productivity
We showed for the first time that structural diversity is an important driver of ecosystem functioning in forests under the controlled conditions of a forest BEF experiment. The structural diversity indices used here may be regarded as measures of canopy complexity (Dănescu et al., 2016;McElhinny, Gibbons, Brack, & Bauhus, 2005), which is increasingly recognized as an important determinant of forest productivity, especially in mixed stands via light-mediated tree-tree interactions (Forrester & Bauhus, 2016). In Sardinilla, three processes were shown to increase light capture and hence productivity of trees growing in mixtures: (a) architectural niche separation, (b) plastic changes in crown shape and (c) temporal niche differentiation among species driven by different leaf phenologies (Sapijanskas et al., 2014).
Our diameter and height diversity indices likely reflect architectural differences and to a certain degree plastic changes in crown shape, the two processes most often evoked as drivers of canopy complexity (Dănescu et al., 2016;Sapijanskas et al., 2014). Phenological differences among tree species played an important role in Sardinilla (Sapijanskas et al., 2014) (Pretzsch et al., 2018), a process that has been found to strongly determine DPRs especially of conservative species (Fichtner et al., 2017). In contrast, the negative neighbourhood structural diversity effects on productivity of fast-growing species (LS) might be caused by greater crown exposure to wind and radiation under dry conditions, which has been described for taller species in other mixtures (Vitali, Forrester, & Bauhus, 2018).
Increasingly divergent effects of neighbourhood structural diversity on individual species may have led to the declining effect of structural diversity on community-level productivity.
Structural diversity had less influence on community productivity than species diversity. While species diversity leads to complementary above-and below-ground interactions, structural diversity effects on forest productivity likely result from above-ground niche partitioning, that is, a subset of the effects of species diversity. As suggested by work in natural forests (Jucker et al., 2015;Zhang et al., 2015), structural diversity effects in our species-rich but rather young tropical plantation are thus likely (partially) mediated effects of species diversity. Crown complementarity as result of intrinsic differences between species is such an effect, which was recently shown to strongly drive DPRs in young tree mixtures (Williams, Paquette, Cavender-Bares, Messier, & Reich, 2017). It is unlikely that at a certain point in time DPRs are driven either by species diversity or by structural diversity. Instead, our results support the idea that the relative contribution of these complementarity effects changes over time. This conclusion is consistent with recent theoretical work and results from many studies, even though the latter were almost exclusively based on two-species mixtures (Forrester, 2014(Forrester, , 2017Forrester & Bauhus, 2016 (Williams et al., 2017) and thus enhance light capture and light-use efficiency (Forrester, 2017;Potvin et al., 2011). Similarly, structural diversity effects on productivity were positive in stands characterized by distinct vertical layering, species with different shade tolerances (Hardiman, Bohrer, Gough, Vogel, & Curtis, 2011;Zhang et al., 2015) or only shade-tolerant species (Dănescu et al., 2016) but negative for shadeintolerant Eucalypt species Ryan et al., 2010).
Structural diversity may hence act as an important driver of positive DPRs if the above prerequisites are met.

| Tree diversity increases production stability
Under climate change, stability of production and other functions of forests are likely to become a key issue in the 21st century. Our results for the oldest tropical tree diversity experiment, like those of Hutchison et al. (2018), clearly support the idea that diversity exerts a positive influence on stability. Tropical mixed-species plantations showed a more stable productivity across periods of contrasting climatic conditions. Hutchison et al. (2018) reported that monoculture mortality in Sardinilla was strongly driven by climatic conditions and that mixing species buffered this effect. Whereas Hutchison et al.
(2018) separated the influence of extreme climatic events on growth and mortality, we focussed on living trees to express stability as an integrated metric, 'temporal stability' (Tilman, 1999), over 10 years of growth. In agreement with our expectation, based on both ecological theory (Loreau & de Mazancourt, 2013) as well as results from the few studies that examined temporal stability (del Río et al., 2017;Jucker et al., 2014), we found species asynchrony to be the strongest predictor of stability. The complementary species interactions that have been described for the Sardinilla experiment (see above) provide the basis for niche differentiation among species and facilitate asynchronous species responses to changing environmental conditions (Loreau & de Mazancourt, 2008). Since asynchrony depends mostly on species characteristics, it may not surprise that we did not find a significant contribution of structural diversity to stability of productivity.
In addition to species asynchrony, stability increased through overyielding (higher mean productivity µ in mixtures) but not through reduced variation (σ) of productivity ( Figure S8). Similarly, Jucker et al. (2014) found a consistent stabilizing effect of overyielding on growth stability across European forest biomes. The lack of a mixing effect on variability of productivity in our study might be partly attributable to the delayed complementarity effects in fivespecies mixtures, combined with the random positioning of the fivespecies plots on especially variable soils (Healy et al., 2008). Finally, species diversity can increase stability via enhanced growth in mixtures but also by reducing drought-induced mortality compared to monocultures (Hutchison et al., 2018).
Overyielding of mixed stands is expected to increase with harsher climatic conditions and more limited water resources, but only if species interactions increase water availability and/or wateruse efficiency (Forrester, 2014;Forrester & Bauhus, 2016). In line with theory, we found the strongest effects of neighbourhood-level species diversity on productivity during a dry period, characterized by a strong and prolonged El Niño drought, while diversity effects were negligible during a wet climatic period. This pattern was apparent at the tree neighbourhood and the community level. The nonlinear drop in overyielding during the extremely wet climate period (P2) likely reflects the absence of a limiting resource (Forrester & Bauhus, 2016) and underlines the climate-induced water availability dependence of overyielding in Sardinilla. The influence of climate on overyielding in the most recent years, however, cannot be disentangled from the effect of increasing species interactions as trees grow. We, therefore, assume that the highest overyielding in the last years was attributable to both increasing strength of interactions with progressing stand development and amplified complementarity during drought. These climate-driven changes in complementary tree-tree interactions at the tree neighbourhood level are likely the principle mechanisms behind the community-level growth responses to contrasting climatic conditions that were reported by Hutchison et al. (2018). Complementary neighbourhood interactions for water are also the likely underlying reason for lower mortality in mixtures when compared to monocultures during the dry period (Hutchison et al., 2018). Complementary water uptake strategies (Schwendenmann et al., 2015) in addition to distinctly different leaf phenologies of the assembled species (Kunert, Schwendenmann, & Hölscher, 2010) may have allowed mixtures to outperform monocultures during drought.
In summary, we found that species and structural diversity enhanced both productivity and its stability in mixed compared to monospecific stands. We show that beneficial effects of diversity in this tropical tree plantation increased with stand development, were highest at the highest levels of diversity and strongest under drought conditions. Results of this study regarding increased productivity in mixtures are consistent with findings from tropical and boreal forests but may not similarly hold in temperate forests or at larger spatial scales, as competitive species interactions and environmental gradients can outweigh beneficial complementarity effects (Chisholm et al., 2013;Fotis et al., 2018;Paquette & Messier, 2011).
Tree-tree interactions in local neighbourhoods were the principle drivers of these diversity effects. For forest restoration initiatives tree-by-tree mixing, compared to commonly used group planting, might, therefore, facilitate positive effects of mixed-species systems on productivity during early stages of stand development (van de Peer et al., 2018). These results support the idea that mixed plantations with species of complementary resource use are a promising strategy for combining high productivity and production stability in the face of unprecedented climate changes.