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

Assessing the overall performance of ecosystems requires a quantitative evaluation of multifunctionality. We investigated plant species diversity effects on individual functions and overall multifunctionality in a grassland experiment with sown monocultures and mixtures comprising four key grass and legume species. Nitrogen fertilisation rates were 50, 150, and 450 kg N ha−1 yr−1 (N50, N150, N450). Ten functions were measured representing forage production, N cycling, and forage quality, all being related to either productivity or environmental footprint. Multifunctionality was analysed by a novel approach using the mean log response ratio across functions. Over three experimental years, mixture effects benefited all forage production and N cycling functions, while sustaining high forage quality. Thus, mixture effects did not provoke any trade-off among the analysed functions. High N fertilisation rates generally diminished mixture benefits. Multifunctionality of four-species mixtures was considerably enhanced, and mixture overall performance was up to 1.9 (N50), 1.8 (N150), and 1.6 times (N450) higher than in averaged monocultures. Multifunctionality of four-species mixtures at N50 was at least as high as in grass monocultures at N450. Sown grass–legume mixtures combining few complementary species at low to moderate N fertilisation sustain high multifunctionality and are a ‘ready-to-use’ option for the sustainable intensification of agriculture.

with the legume being either T. pratense or T. repens, and L. perenne grown in the same plot acting as the reference plant. Apparent N transfer from legume to grass species was calculated following Vallis et al. 3 : where 'grassmix' are the grass species grown in mixtures, and 'grassmono' are the grass species grown in monocultures adjacent to the mixture plot and at the same N level. The amount of N (kg N ha -1 yr -1 ) of Ndfa and Ntrans was calculated using the N content in the species and their respective dry mass. The two fractions Ndfa and Ntrans were summed to give Nsym. See Nyfeler et al. 4 for full details of calculations and a critical appraisal of the methods.

NO3 in soil solution
Sampling plots for NO3 in soil solution were monocultures, the four-species equi-proportional mixture, and the dominant mixtures (Table S1,

Metabolisable energy
Metabolisable energy content (ME) is defined following a reference manual of Agroscope 5  where CP denotes crude protein content and CP OS denotes crude protein content per organic matter (OM) (i.e. CP OS = g CP/kg OM). CP and OM are determined as described in the main text.
ME is then computed as:

Data analyses
The univariate diversity interaction model Kirwan et al. 6 introduced the univariate diversity interaction model using multiple linear regression. The approach allows modelling of the response variable as a function of species proportions, the latter summing to unity in each stand (see also Cornell 7 for extensive information). The community-level response is modelled as a linear combination of (1) identity effects of species as given by their monoculture performance, (2) species net interactions, termed diversity effects (D: being positive, negative or neutral), which are defined as the difference between the actual mixture performance and that expected from the relative contribution of the constituent monocultures, and (3) any further variables, such as overall sowing density or N fertilisation. Given this framework, a basic univariate diversity interaction model to one function of our data is: The α coefficient denotes the effect of changing the overall sowing density on the response variable y, e.g. yield. With P i denoting the species' proportions in a stand, coefficients β 1 to β 4 estimate the effects of species' proportional contributions on y and, if P = 1, β coefficients estimate the response y of species' monocultures. The γ coefficient estimates the effect of N fertilisation, with N_Treat being a factor with three levels: N50, N150, and N450. Coefficients δ 1 to δ 6 estimate the six possible pairwise interactions among the four species to model diversity effects. The residual term ε is assumed to be normally distributed with constant variance σ 2 . Interactions between e.g. N_Treat and P i can be added to eqn. S1.

Bootstrap sampling
Parametric bootstrapping was based on the estimated coefficients of the final model (eqn . 3 main text) and related residual variances of functions with corresponding co-variance matrix Σ.
Because residuals were correlated, we randomly sampled from a multivariate normal distribution with a zero-vector of means and given the residual variances and Σ, as estimated by eqn. (3). This allowed for calculation of new data (fitted values from model + randomly sampled 'residuals') to which eqn. (3) was fitted, followed by calculation of the log response ratio (LRR) and overyielding as defined in the main text. The procedure was repeated 1000 times. Functions that were natural log transformed in eqn. (3), namely stability, weed biomass and NO3 in soil solution, were back-transformed to linear scale before calculating the LRR and overyielding.

Analyses of single year's data
Analyses of the data of each experimental year followed the same principles as the data averaged across years (eqs. 1-5 main text). The following details need to be mentioned: Regarding the variable P i , which denotes species proportions in a stand (eqs. 1-3), the sowing proportions were used as a predictor for analyses in year 1, whereas observed species proportions in annual biomass yield of the preceding year were used in years 2 and 3. Using the previous year's species proportions as predictors was done to overcome effects of year-to-year changes in community composition.
In all analyses of single years, the functions weed biomass and NO3 were first natural log transformed and then scaled to range between 0 and 100% to achieve a multivariate normal residual distribution.
Regarding sowing density, increased density only had an effect in year 1 and only on weed biomass. Consequently, the variable sowing density was omitted from the models for years 2 and 3. At year 1, increased sowing density resulted in 24% less weed biomass in the fourspecies equi-proportional reference mixture at all N fertilisation treatments (t361 = 2.45, P = 0.015). Notably, the weed reduction occurred on low absolute levels (

Effect of log transformation of functions on multifunctionality
We further evaluated the effect of the natural log transformation of the functions stability, weed biomass, and NO3 in soil solution on the outcome of the MLLRD. Log transformation of data in regression analyses can introduce bias in that the geometric mean is predicted on the back-transformed linear scale, the geometric mean being always equal or smaller than the linear mean of untransformed data. Because the MLLRD as defined in eqn. (7), main text, is a mean of several log response ratios, the effect of log transformation of single functions on the MLLRD is not directly obvious. To assess this aspect, we computed the functional responses and the MLLRD as described in the main text but left the functions stability, weed biomass, and NO3 untransformed. We calculated the MLLRD for overall legume proportions of 0, 0.25, 0.5, 0.75 and 1 in mixtures, at the three N fertilisation treatments. It turned out that the regression with untransformed data resulted in (unrealistic) negative mixture predictions for weed biomass and NO3 in one fifth of cases, which made a comparison of the MLRRD impossible in these situations. Where a comparison was possible, the MLRRD based on the regression with untransformed data was mostly larger than or equal to the MLRRD calculated with the three functions log transformed. We thus concluded that our statements regarding the diversitymultifunctionality relationship and the resource use-multifunctionality relationship were generally conservative.

Test of the effect of N fertilisation on multifunctionality (MLLRN)
The effect of N fertilisation on multifunctionality was tested by calculating the individual functions' LRR as defined by eqn. (8) followed by performing a bootstrap sampling taking among single LRRs, which was implemented in a generalised least square regression to test the MLRRN against zero. We performed three distinct comparisons related to N fertilisation: fourspecies equi-proportional mixture at N50 against i) the average of the two grass monocultures at N450, ii) the average of all monocultures at N450, iii) the four-species equi-proportional mixture at N450; therefore, three bootstrap runs were performed, one for each case. Table S1. Species proportions in monocultures and mixtures sown following a simplex design. All communities were sown at two overall densities and subjected to three levels of N fertiliser application (50, 150, 450 kg N ha -1 yr -1 ), except the binary and codominant mixtures, which were established only at 150 kg N ha -1 yr -1 .       *** P ≤ 0.001, ** P ≤ 0.01, * P ≤ 0.05, † P ≤ 0.1, ns: not significant ‡ Lp: L. perenne, Dg: D. glomerata, Tp: T. pratense, Tr: T. repens § The diversity effect is calculated as the difference between the performance of the equi-proportional four-species mixture and the average of monocultures and is split into the effects of mixing grass and legume species (DBGL effect), mixing the two grasses Lp and Dg (Lp•Dg effect), and mixing the two legumes Tp and Tr (Tp•Tr effect). ¶ Standard errors (s.e.) were calculated from the weighted average of variances at each N fertilisation level. # Weed biomass and NO3 in soil solution were natural log transformed and subsequently scaled to range between 0 and 100%; thus, scaling with the maximal functional performance is not reasonably possible. See Fig. S4 for back-transformed values on the linear scale. ‡ ‡ Lp•Dg and Tp•Tr effects and were omitted for Nsym and NO3 in soil solution to avoid the estimation of negative monoculture values (Nsym) and convergence problems (NO3). Figure S1. Percent of beneficial mixture effect (mixture performance greater than the average of monocultures) of the four-species equi-proportional mixture for the three forage quality functions scaled per hectare at three N fertilisation treatments (N50: 50 kg N ha -1 yr -1 , N150: 150 kg N ha -1 yr -1 , N450: 450 kg N ha -1 yr -1 ). Data were averaged across years. OM: organic matter. Further details are explained in Fig. 2, main text.    Figure S5. Percent of beneficial mixture effect (mixture performance greater than the average of monocultures) of the four-species equi-proportional mixture at three N fertilisation treatments for seven (a) and ten functions (b, c) over three experimental years (N50: 50 kg N ha -1 yr -1 , N150: 150 kg N ha -1 yr -1 , N450: 450 kg N ha -1 yr -1 ). Point estimates are based on multivariate linear mixed-effects regression (Tables S3, S4, S5) and error bars represent the 95% confidence intervals (CIs). Functions whose CI does not include 0 can be considered to reveal a significant mixture effect. Weed biomass and NO3 in soil solution were backtransformed to linear scale to calculate the beneficial mixture effect. Extreme CIs are truncated and reported as numbers, and X axes are log-scaled to equalise distances on both sides of parity. Note the different scale in (a).  Figure S6. Observed versus predicted values of the ten function based on multivariate regression analyses following eqn. 3, main text. Prior to analyses, stability, weed biomass, and NO3 in soil solution were natural log transformed to achieve a multivariate normal distribution of residuals, and all functions were standardised to range between 0 and 100%.