Multispecies for multifunctions: combining four complementary species enhances multifunctionality of sown grassland

We used a dataset from a grassland diversity experiment at Zürich-Reckenholz, Switzerland, in the Atlantic central climatic zone of Europe. The data contain measurements on many functions from 78 plots that comprised monocultures and mixtures sown at a wide range of species relative abundances, set up at three levels of N fertiliser application and maintained for 3 years following establishment, which is a typical time in grassland-crop rotations.

Monocultures and mixtures were sown following a simplex design⁶⁴. Four perennial species, known to be key forage species in ruminant production, were selected based on the factorial combination of their functional traits related to temporal establishment (fast-establishing vs. temporally persistent), and N acquisition (non-fixing for grasses, N₂-fixing for legumes). The species were Lolium perenne L. cultivar (cv.) Lacerta (fast-establishing grass), Dactylis glomerata L. cv. Accord (temporally persistent grass), Trifolium pratense L. cv. Merviot (fast-establishing legume), and Trifolium repens L. cv. Milo (temporally persistent legume). The type of stands were: monocultures (100% of one species), binary mixtures (50% of each of two species), an equi-proportional mixture (25% of each of the four species), dominant mixtures (70% of the dominant species, 10% of each of the other three), and co-dominant mixtures (40% of each of two species, 10% of each of the other two; see Supplementary Table S1). All types of stands were sown at two levels of overall sown density, with the high level being the recommended seed weight (100%) under conditions typical of Switzerland, and the low level being 60%.

The experiment was sown in August 2002 on plots of 3 m × 6 m and was maintained from 2003 (year 1) to 2005 (year 3). The plots were fertilised with N fertiliser (as NH₄NO₃) at rates following a geometric series: 50, 150, or 450 kg N ha⁻¹ yr⁻¹ (N50, N150, and N450, respectively), split into five equal applications. In early spring, all plots received phosphorus and potassium in amounts expected to be non-limiting for intensively managed grasslands on fertile soils in Switzerland. At the N150 treatment, all types of monocultures and mixtures were established, whereas the N50 and N450 treatments only included the monocultures, the equi-proportional mixture, and the dominant mixtures. The 78 plots were arranged in a fully randomised design. Consult Nyfeler et al.³⁷ for full details of the experimental design, establishment, and maintenance.

Ten functions were measured representing (i) forage production: aboveground biomass yield, standard deviation of yield, temporal stability, weed biomass; (ii) N cycling: symbiotic N₂ fixation, N efficiency, NO₃ in soil solution; and (iii) forage quality: crude protein content, organic matter digestibility, metabolisable energy content (Table 1). To date, detailed analyses from the experiment have been published on two functions, namely biomass yield³⁷ and symbiotic N₂ fixation³³.

Measurement of functions

Aboveground biomass yield and weed biomass

All plots were harvested five times annually at 5 cm above ground surface. Aboveground biomass yield at each harvest was determined by drying a representative subsample to constant weight (65° C for 48 h), and this data was summed to give total annual biomass yield. Biomass proportions of the four sown and pooled unsown species (weeds) were measured by manually separating samples from permanent sub-plots (0.8 m × 0.3 m), which was done at the first, third, and fifth harvest of each year. These data allowed for calculation of weed biomass per ha and year.

Standard deviation and stability of yield

Year-to-year standard deviation of yield (SD_yield) was calculated from the annual yields of the three experimental years, and stability was defined as the ratio of averaged annual yields to year-to-year SD_yield (following Lehman and Tilman⁶⁵). To measure yield variation within each year, seasonal SD_yield was calculated from the five annual harvests, and seasonal stability was defined as the ratio of total annual yield to seasonal SD_yield. We purposely use both SD_yield and stability as both measures are essential to evaluate yield variation⁶⁶.

Symbiotic N ₂fixation

Symbiotic N₂ fixation (N_sym) was determined by the isotope dilution method⁶⁷. Double-labelled ¹⁵N-enriched ¹⁵NH₄¹⁵NO₃ was applied on a permanently defined, central part of each plot (1.4 m × 1.5 m). Plant samples were analysed for ¹⁵N and ¹⁴N abundance by gas isotope ratio mass spectrometry and by thermal conductometry. N_sym in the sward, as calculated here, comprises legume N derived from the atmosphere (N_dfa) plus N derived from apparent N_dfa transfer to the grass (N_trans). See Supplementary Appendix S1 and Nyfeler et al.³³ for full details of measurements and calculations.

N efficiency

N efficiency was defined as the ratio of total N yield to the amount of applied fertiliser N and therefore measures the total N output of the system in relation to the fertiliser N input. Total N yield was calculated by first multiplying N content from biomass samples with their total dry mass to give the N yield per harvest. Annual total N yield was then computed as the sum of all harvests.

NO ₃in soil solution (NO₃)

Porous cup tension lysimeters were installed to extract soil water from a depth of 60 cm below ground surface. In 2-week intervals from October 2004 to April 2006, a suction of 80 kPa was applied 1 day prior to sampling, and concentrations of nitrate–N (NO₃-N) were determined by spectrophotometry. We note that NO₃ data were only available for years 2 and 3. See Supplementary Appendix S1 for details of the measurements.

Crude protein content 

Crude protein content (CP) in stand biomass was calculated from the N content in biomass samples, multiplied by 6.25. The justification for the multiplicative factor is given by the fact that all biological proteins contain on average 16% N⁶⁸.

Organic matter digestibility

Organic matter digestibility (OM digestibility) was determined from biomass samples of the second and fourth harvest following the two-stage in vitro fermentation process with rumen liquor and acidic pepsin solution according to Tilley and Terry⁶⁹; see Supplementary Appendix S1 for details. Information on OM digestibility was only available for years 2 and 3 of the study.

Metabolisable energy content

Metabolisable energy content (ME) of stand biomass was calculated based on OM digestibility and CP following a reference manual of Agroscope⁷⁰; see Supplementary Appendix S1 for calculation. Due to the connection with the measurement OM digestibility, ME data were only available for years 2 and 3.

Data for each function were computed at the plot level for each of three experimental years (the three exceptions as noted). For analyses across years, data was averaged across available years, except SD_yield and stability (see above).

Data analyses

We applied the multivariate modelling framework¹³ to estimate simultaneously species identity and diversity effects of the ten functions along with effects of N fertilisation. To allow direct comparisons of the model terms, all functions’ data were standardised to a common scale by dividing them by their maximum value (at a single year) over the 3-year experiment and N fertilisation treatments. This scaling allowed for a direct comparison of results among years. Note that the multivariate approach is a generalisation of the univariate diversity interaction model⁶¹, and we refer to Supplementary Appendix S1 for a summary to the univariate regression.

In the following, we generally refer to the analysis of data averaged across experimental years, and all equations model the response at a single plot (plot subscripts are omitted). A preliminary regression equation was specified for the kth function (k = 1–10) with:

The α coefficient denotes the effect of changing sowing density on the response variable y_k, for example biomass yield, with DENS being coded as − 0.5 and 0.5 for the low and high sowing density, respectively, so that all other parameters give the response y_k at average density. Variables P_i denote the species’ sown proportions in a stand. Coefficients β_1fk to β_4fk estimate the effects of the four species’ proportional contributions on y_k (identity effects) and, if P_i = 1, β coefficients estimate the response y_k of species’ monocultures. Identity effects β_ifk are estimated at each N fertilisation treatment f (factor N_Treat with three levels: N50, N150, and N450), which is equivalent to specifying an interaction between N_Treat and P_i. Coefficients δ_ijk estimate the six possible pairwise interactions among the four species to evaluate diversity effects. The residual term ε_k is assumed to be normally distributed with constant variance σ²_k.

Dooley et al.¹³ extended Eq. (1) to a multivariate model using matrix notation, where response variables y_k constitute a matrix of k columns. Yet for parameter estimation, the multivariate matrix notation can be re-written applying principles of linear mixed-effects regression⁷¹, leading to:

Here, the response variable y_k denotes a column vector, in which performances of all functions k are listed. The variable FUNC is a factor with ten levels, one for each function k. Consequently, predictor variables DENS, P_i, P_iP_j, and N_Treat (with meanings as explained) are repeated k times within columns of the design matrix, and corresponding coefficients are estimated as fixed parameters (see Dooley et al.¹³ for an example). The variable ‘Plot’ (also repeated k times) estimates the plot-specific, common variance of the functions per plot (random intercept). The term e_k ~ MVN(0, Σ), where MVN denotes multivariate normal, with mean 0 and co-variance matrix Σ among functions. For parameter estimation, the residual parameter was defined as Var(e_k) = σ²δ_k², with δ being a ratio to represent k variances (see Pinheiro and Bates⁷¹ p. 209 for details), and an unstructured co-variance matrix was imposed on the residuals to estimate Σ.

Applying Eq. (2), it turned out that sowing density had no significant effect (t < 1.65, P > 0.10 for all functions but one), and it was omitted from all further models. Moreover, to achieve a multivariate normal residual distribution, the functions stability, weed biomass and NO₃ were first natural log transformed and then divided by their maximum value to range between 0 and 100%. Given this amendment, residuals showed no evidence of a deviation from multivariate normality (approved by Mardia’s multivariate normality test⁷²).

Equation 2 estimates six δ_ijk coefficients per function (diversity effects). To increase parsimony, a series of hierarchical models were constructed as described in detail in Nyfeler et al.³⁷ and Helgadóttir et al.⁵⁵. Applying likelihood ratio tests for the comparison of nested models (see Pinheiro and Bates⁷¹ p. 83), it appeared that the six diversity effects could be grouped together to represent specific interactions between grass and legume species (D_BGL), and interactions between the two grass and between the two legume species. Moreover, given known effects of N fertilisation on species diversity effects^37,50, it was tested whether the N fertilisation treatment interacts with the (pooled) diversity effects, in which case it turned out that interactions of N_Treat with the D_BGL term were highly significant, but not with the grass-grass and legume-legume terms, which led to:

with Lp: L. perenne, Dg: D. glomerata, Tp: T. pratense, Tr: T. repens, and D_BGL = P_LpP_Tp + P_LpP_Tr + P_DgP_Tp + P_DgP_Tr, representing the four pooled pairwise interactions between grass and legume species. All other variables and their related regression coefficients have meanings as explained. The marginal and conditional R² (following Nakagawa and Schielzeth⁷³) of this regression was 0.876 and 0.881, respectively, which led us to conclude that predictions based on Eq. (3) were highly reliable (see also Supplementary Fig. S6 for observed versus predicted values of the ten functions based on Eq. 3, and Supplementary Table S2 for goodness-of-fit measures for selected models).

We choose the four-species equi-proportional mixture as a reference mixture to evaluate beneficial effects on functional performance in mixtures as compared to the average of monocultures (overyielding: OY) using the estimated coefficients of the final model (Eq. 3):

with ŷ_{equi_k} being the predicted functional performance of function k at the four-species equi-proportional mixture and ŷ_{avemono_k} the predicted performance of the average of monocultures. The 95% confidence interval (CI) to overyielding was calculated by parametric bootstrapping⁷⁴. Because we intended to achieve an approximate multivariate normal distribution of the bootstrap sample, the procedure was performed with the closely related log response ratio (LRR):

with meanings of components as explained, and lastly the LRR was rescaled to OY to give the confidence intervals (see Supplementary Appendix S1 for details to the bootstrap sampling). We note that eqs. (1–5) were also applied to data of each individual year (year 1–3), and we refer to Supplementary Appendix S1 for details of the single years’ analyses.

Finally, the diversity-multifunctionality relationship was evaluated using the mean LRR (MLRR) across all functions. This follows the reasoning that a greater number of functions with higher LRRs indicates enhanced multifunctionality in mixtures as compared to monocultures. To this aim, the LRR as defined in Eq. (5) was generalised to:

Here, ŷ_{mix_k} is the predicted functional performance of function k of any mixture and ŷ_{ave_w_mono_k} the predicted performance of the weighted average of monocultures (all based on Eq. 3), the weights being the species proportions in the mixture. For functions where minimal values were regarded as of positive benefit (SD_yield, weed biomass, and NO₃), their LRR was multiplied by − 1. Calculation of LRR_k was then followed by computation of the MLRR across functions:

Note that the MLRR_D reflects a change in multifunctionality in dependence on the components of the ratio (here: mixtures versus averaged monocultures, reflecting a change in species diversity). We prefer such a metric over absolute measures of multifunctionality, as absolute measures are highly context specific and have little value for comparison among systems. To provide a statistical inference to the MLRR, t tests against zero were not applicable because the single LRRs must be assumed to be correlated among functions for any given community. Instead, we used generalised least squares regression and implemented the correlation matrix among single LRRs that could be derived from the bootstrap sampling.

Justified by the outcome of Eq. (3), the MLRR_D was calculated for a range of overall legume proportions with equal proportions of the two grass and the two legume species to display the effect of species diversity and grass–legume interactions on multifunctionality. The range of legume proportions for which the MLRR_D was significantly different from zero was calculated using the Johnson Neyman technique³⁸. See Supplementary Appendix S1 for details of the procedure, including the range test.

The MLRR approach is flexible in that it allows the consideration of different comparisons, reflected by the components of the corresponding LRRs. We directly tested the effect of N fertilisation on multifunctionality using comparisons of the four-species reference mixture at N50 with three selected types of communities at N450:

where ŷ_{equi_k_N50} is the predicted functional performance of function k at the four-species equi-proportional mixture at N50, and ŷ_{community_k_N450} is the predicted functional performance of a community at N450 (both based on Eq. 3), namely either (i) the average of the two grass monocultures, (ii) the average of all monocultures, or (iii) the four-species equi-proportional mixture. For functions in which minimal values were targeted, their LRR was multiplied by − 1. See Supplementary Appendix S1 for details of the inference to the MLRR_N. All analyses were performed using the statistical software R, version 3.6.1⁷⁵ and the package nlme for linear mixed-effects models⁷⁶.

Source: Ecology - nature.com