Study site

The study was conducted in the Gutianshan National Nature Reserve (29°140' N; 118°070' E) in Zhejiang Province, south-east China. The reserve, situated in mountainous terrain (250–1260 m above sea level), covers ca. 8000 ha of semi-evergreen, broadleaved forest. The climate is subtropical, with a mean annual temperature of 15.3 °C and mean annual precipitation of ca. 2000 mm.

In 2008, we established 27 study plots (30 m × 30 m) following a stratified random selection design based on the stand age ( < 20 to > 80 years) and woody plant species richness (25–69 species) typically encountered in the reserve48. We excluded one of the study plots from our analyses, because extreme values indicated processing errors of the soil samples used for the measurement of several belowground ecosystem functions that would have biased the statistical analyses (see Supplementary Fig. 10).

Ecosystem functions

We used 22 independent measurements of nine ecosystem functions that play important roles in energy and nutrient flow across trophic levels in our study system. These nine functions were: erosion control, microbial activity, primary productivity, nitrogen cycling, leaf decomposition, wood decomposition, herbivory resistance, predation, and parasitism (Supplementary Table 1).

Erosion control was quantified as the inverse of the estimated erosivity of rain drops. Microbial activity comprised measurements of enzyme activity and microbial biomass. We measured the activity of four enzymes that play important roles in decomposition processes and nutrient cycling by degrading cellulose (β-xylosidase), xylan (xylosidase), chitin (N-acetyl-glucosaminidase), and polyphosphates (acid phosphatase). Primary productivity was estimated as the increment in stem basal area of woody plants over four years. Soil nitrogen cycling was assessed by measuring net ammonification and net nitrification. The quantification of leaf decomposition comprised measurements of community-level and species-specific leaf decomposition rates. Wood decomposition was assessed by two complementary approaches: calculation of community wood decomposition rates and measurement of decomposition rates of standardized wood samples. Herbivory resistance was based on measurements of leaf damage by herbivorous insects of canopy trees and woody plant saplings. Predation was measured for two groups of common arthropod predators, ants and predatory wasps. Parasitism rates were assessed for parasitic Hymenoptera. See Supplementary Table 1 for full details on all measurements.

Our study provides a comprehensive analysis of these functions in a multifunctional and multitrophic context, which differs substantially in scope from previous analyses in our study system that looked at individual functions and their relationship with woody plant communities (Supplementary Table 11).

Biotic attributes

We conducted a complete inventory of the woody plant communities of each plot (abundance, species composition), with all tree and shrub individuals > 1 m height, in 2008. We estimated the stand age of the study plots from tree stem cores and diameter at breast height (DBH) measurements (to the closest 0.1 mm)48. To characterize the functional composition of the woody plant communities, we measured a range of leaf morphological, leaf chemical, and wood anatomical traits of the woody species: leaf area (LA), specific LA and LDMC, leaf C content, leaf carbon to nitrogen (C:N) ratio, leaf polyphenolics content, wood density, mean xylem vessel diameter, and wood fiber wall thickness. These traits determine key aspects of the life history strategies of plants44 and have previously been shown to strongly affect primary productivity, decomposition, herbivory, and other important ecosystem functions in subtropical and other forests46,52,53. Leaf traits were measured for ca. 80% of the 147 woody plant species recorded on the study plots and these species represented 95% of the total number of tree and shrub individuals at the study sites. Samples for trait measurements were taken from sun-exposed leaves of five to seven plant individuals in total54, collected from up to seven plots per species in the summer of 2008. Trait measurements followed standardized protocols55. Wood traits were available for 93 species that represented 83% of the total number of woody plants at the study site. Wood samples were taken from sun-exposed branches with a diameter of at least 3 cm (one to three samples per species), cut with a sliding microtome, permanently fixed, and analyzed under a microscope56.

Species richness of heterotrophic organisms was assessed with a range of different methods between 2008 and 2012. We considered 11 groups of arthropods in our analyses, which we assigned to different functional groups according to their feeding ecology and trophic rank: parasitoids (parasitic wasps [Hymenoptera: Braconidae, Chrysididae, Eurytomidae, Ichneumonidae, Leucospidae, Mutillidae, Pompilidae, Trigonalyidae]), predators (spiders [Arachnida: Araneae], centipedes [Chilopoda], cavity-nesting solitary wasps [Hymenoptera: Pompilidae, Sphecidae, Vespidae], and strictly predatory as well as omnivorous ants [Hymenoptera: Formicidae]), primary consumers/herbivores (moth and butterfly caterpillars [Lepidoptera], weevils [Coleoptera: Curculioninae], longhorn beetles [Coleoptera: Cerambycidae], and bark beetles [Coleoptera: Scolytinae]) and decomposers (millipedes [Diplopoda] and isopods [Isopoda]). Epigeic spiders, centipedes, epigeic ants, weevils, isopods, and diplopods were sampled with pitfall traps38 (four traps per plot from March to September 2009). Lepidopteran larvae, arboreal spiders, and ants were sampled from understory trees and shrubs by means of beating57 (25 saplings per plot on three sampling dates in 2011 and 2012). Cavity-nesting predatory wasps and associated parasitic wasps were sampled with reed-filled trap nests (Supplementary Table 1). Longhorn beetles, bark beetles, and canopy ants were sampled with flight interception traps (4 traps per plot from May to August 2010). In addition, ants were sampled with standardized protein and carbohydrate baits58 (36 baits per plot in May 2012). All arthropods were identified to species or morphospecies. Data on soil fungi were obtained from soil cores20 (eight cores of the upper 10 cm of soil per plot, taken in September 2012). The soil was sieved and freeze-dried for molecular analysis. Fungal DNA was extracted with the MoBio soil DNA extraction kit and analyzed by pyrotag amplicon sequencing of the fungal internal transcribed spacer (ITS). Sequence datasets were quality filtered, normalized and clustered into species-level operational taxonomic units (OTUs). Non-target taxa OTUs as well as singletons, doubletons, and tripletons were removed from the dataset. The fungal reference sequences were assigned to ecological functional groups (saprophytes, pathogens, ectomycorrhizae, arbuscular mycorrhizae) on the basis of sequence similarity using the default parameters of the GAST algorithm59 against the functional reference dataset60.

Abiotic attributes

We measured a range of spatial and environmental variables in the study plots that might directly or indirectly influence ecosystem functions. Elevation (m above sea level), slope (°), degree of northness and eastness (cosine- and sine-transformed radian values of aspect), latitude, and longitude were assessed during plot establishment in 2008. Soil pH (measured potentiometrically in a H 2 O suspension), soil nitrogen (N) and carbon (C) contents (measured with Vario ELIII elemental analyzer, Elementar, Hanau, Germany), and the soil C:N ratio were determined from a bulk sample of nine soil cores (0–10 cm) per plot, taken in summer 2009. Mean annual temperature and mean January and July temperatures per plot were obtained from continuous measurements with HOBO data loggers (one data logger in the center of each plot; 30 min time intervals from July 2011 to June 2012).

Calculation of multifunctionality

Each of the 22 measurements of ecosystem functions were scaled to range from 0 to 1 with the formula f(x) = (x i − x min )/(x max − x min ), where x is the variable of interest with its minimum (x min ) and maximum (x max ) values observed across all study plots (see Supplementary Table 12 for details on observed values). Data on herbivory and erosion were multiplied by −1 prior to scaling to represent herbivory resistance and erosion control. All scaled measurements of a given function were then averaged per plot to obtain an ecosystem function variable that represents the mean of the various independent measurements, giving each function the same weight in the multifunctionality analyses. This yielded nine variables corresponding to the functions described above, on which all analyses were conducted.

From the various methods available to calculate multifunctionality, we chose two of the most commonly used: the “averaging approach” and the “multiple threshold approach”5. The averaging approach takes the mean value across all standardized functions as an index of multifunctionality for each study plot5. Threshold approaches measure how many functions simultaneously exceed a predefined percentage of the maximum observed value of each individual function9,17,18. As the selection of a given threshold is arbitrary, analyzing multiple thresholds of maximal functioning is recommended5. Thus, we used thresholds of 30%, 60%, and 90% to analyze how diversity affects multifunctionality at low, medium, and high thresholds, respectively, of the observed maximum functioning (see refs. 7,19 for similar approaches). We used the mean of the three largest values of each function as the observed maximum to reduce the impact of potential outliers5. The number of functions surpassing a given threshold was calculated with the R-package “multifunc”.

Trophic-level species richness and community composition

We calculated a set of predictors that aggregated the information available for the diversity and community composition of woody plants and for the diversity of heterotrophic organisms at higher trophic levels (see Supplementary Table 12 for details on observed values).

From the woody plant inventory of each study plot, we calculated species richness and species composition of woody plants. Species composition was expressed as the first two principal components (PCs) of a principal components analysis (PCA) on relative abundance data, following the approach of ref. 39. PCAs were conducted with the R package “vegan”. Based on the relative abundance of the woody plant species, their morphological and chemical leaf traits, and their anatomical wood traits (see above), we quantified leaf morphological diversity, leaf chemical diversity, and wood anatomical diversity of the forest stands. We used Rao’s quadratic entropy Q61 for this, which describes the variance in pairwise trait dissimilarities among all individuals in a community. In addition, we calculated the community weighted mean values (CWMs) of the leaf and wood traits, i.e., the mean trait values across all individuals of a community. CWMs were subjected to a PCA for dimensionality reduction, yielding three PCs that accounted for the mean trait composition of the woody plant communities (Supplementary Table 6). Rao’s Q and CWMs were calculated with the R-package “FD”. To control for confounding effects of the stand structure on diversity and composition effects, we included stand age and the density of woody plants in the study plots as further predictors.

To obtain metrics of heterotrophic species richness that are comparable across trophic levels and that allow combining data assessed with different sampling methods, we followed the approach of ref. 62 and calculated a diversity index that merges the species richness patterns of individual taxa. This was done by averaging the scaled (scaled to the maximum observed value across the study plots) species richness of the taxa in each functional group (i.e., parasitoids, predators, herbivores, macrofaunal decomposers, saprophytic fungi, parasitic fungi, mycorrhizae).

Finally, we included information on the environmental plot conditions as predictors, as they might have an influence on the ecosystem functions considered in our study63. The environmental variables were subjected to a PCA and we retained the first two PCs as a characterization of elevational differences in aboveground temperature (PC1) and belowground soil conditions (PC2; Supplementary Table 13).

Effects of trophic-level diversity on multifunctionality

We tested for multicollinearity among the diversity and community composition variables and excluded the first PC of woody plant composition from the analyses (highly correlated with leaf chemical diversity, Pearson’s r = 0.78, P < 0.001), and the second PC of trait composition (highly correlated with leaf morphological diversity, Pearson’s r = 0.72, P < 0.001).

We used multiple linear regression on plot-level data with either the average multifunctionality index or the number of functions larger than the selected threshold as response variables. For the regression analyses, we used a model selection and averaging approach that calculated all possible subset models with up to four predictors and chose from this set those subset models with the lowest values (ΔAICc ≤ 2) of the Akaike Information Criterion corrected for small sample size (AICc). As potential predictors, we considered diversity, composition, stand age, and environmental variables described above. We also considered potential interactions between stand age and diversity metrics in those cases where an initial data check indicated the possibility of significant or close to significant interactions (decomposer species richness × stand age, predator average species richness × stand age) and in the case of our design variable (tree species richness × stand age). For the selected set of models we used model averaging to address model uncertainty64. Model selection and multimodel averaging were conducted with the R-package ”MuMIn”. Adjusted R²-values and partial R² values for individual predictors were averaged over all selected models. Model residuals were checked for normality and homogeneity of variances.

In addition to the analyses of multifunctionality based on the full set of functions and predictors across trophic levels, we conducted two further analyses to tease apart diversity effects within trophic levels vs. effects across trophic levels. As results for the above analyses on average and threshold multifunctionality showed the same general effects, we restricted the additional analyses to average multifunctionality for simplicity.

First, we calculated an overall index of total community richness as the average of the standardized diversity indices of (a) all trophic levels (average total community richness) and (b) all heterotrophic levels (average heterotrophic community richness) considered in our study. We re-ran the regression analyses as described above (AICc model selection and averaging based on a maximum of four predictors per model), but with one of the two community richness indices (and its interaction with stand age) instead of the individual diversity indices of the individual trophic levels as a predictor. Analogous to the regression analyses described above, environmental variables and plant-based predictors were considered as covariates. This analysis sheds light on the extent to which overall biodiversity (instead of the diversity of individual trophic levels) influences multifunctionality.

Second, we re-ran the analyses by trophic level (parasitoids, predators, herbivores, plants, decomposers, and fungi—the latter were considered together because functions such as microbial biomass and activity are the outcome of the combined effects of different soil fungi). For each trophic level, we excluded those functions from the multifunctionality index that were directly mediated by this trophic level (parasitoids: parasitism; predators: predation; herbivores: herbivory resistance; plants: primary productivity; decomposers: leaf and wood decomposition; fungi: microbial activity, nitrogen cycling, decomposition), as these functions might disproportionately influence the relationship between multifunctionality and diversity65,66. As predictors, we only considered the species richness of the respective trophic level, as well as environmental predictors and general plot characteristics (stand age, tree density) as covariables. These analyses allowed us to assess the extent to which diversity effects of a given trophic level propagate though the food web to affect ecosystem multifunctionality. We note that these comparisons across trophic levels are interdependent, which we cannot completely avoid with our relatively small sample size67,68. Results therefore need to be interpreted with care (see Results section for further details).

Finally, we used the same regression and model averaging approach with each of the nine individual ecosystem functions as a response variable and the set of environmental variables and the diversity metrics of all trophic levels as predictors. This allowed us to compare the multifunctionality results to the performance of individual ecosystem functions.

Path analysis

We used path analysis69 to shed light on potential causal direct and indirect pathways that determine the combined effects of multiple metrics of biodiversity across trophic levels. The path model was informed by the outcome of the regression analysis of average multifunctionality. We included all pathways from plant-based predictors (exogenous variables, the three variables retained in the model had correlations with Pearson’s r ≤ 0.2) to heterotrophic predictors and multifunctionality (endogenous variables), as well as pathways between heterotrophic predictors and from these predictors to multifunctionality. We did not perform any model simplification and assessed model fit based on χ2-statistics. Models were fitted with the “lavaan” R-package.

All analyses were conducted in R 3.3.1 (www.r-project.org).

Data availability

All taxon data is available on the BEF-China project database. The pyrosequencing dataset of soil fungi is deposited in the EMBL SRA database under study number PRJEB8979.