The links between ecosystem multifunctionality and above- and belowground biodiversity are mediated by climate.

Plant biodiversity is often correlated with ecosystem functioning in terrestrial ecosystems. However, we know little about the relative and combined effects of above- and belowground biodiversity on multiple ecosystem functions (for example, ecosystem multifunctionality, EMF) or how climate might mediate those relationships. Here we tease apart the effects of biotic and abiotic factors, both above- and belowground, on EMF on the Tibetan Plateau, China. We found that a suite of biotic and abiotic variables account for up to 86% of the variation in EMF, with the combined effects of above- and belowground biodiversity accounting for 45% of the variation in EMF. Our results have two important implications: first, including belowground biodiversity in models can improve the ability to explain and predict EMF. Second, regional-scale variation in climate, and perhaps climate change, can determine, or at least modify, the effects of biodiversity on EMF in natural ecosystems.

Ecologists have explored the links between biodiversity and ecosystem functioning for more than two decades1234. Yet, most of our understanding of biodiversity–ecosystem function (BEF) relationships comes from experimental or observational investigations of how a single ecosystem function, usually aboveground net primary productivity, depends on plant biodiversity in small plots at single sites245 (but see refs 3, 6 for multisite comparisons). Clearly, ecosystems perform multiple functions simultaneously (ecosystem multifunctionality; EMF)7891011, and these functions can be mediated, if not controlled, by biodiversity in other components of ecosystems, such as the biodiversity of taxa belowground11121314151617. BEF relationships might also vary geographically, perhaps because of spatial or temporal variation in climate or other abiotic factors35918.

Although above- and belowground communities and their interactions can shape multiple ecosystem functions, it is difficult to tease apart their relative effects16, perhaps because taxa operate at different spatial scales. For instance, belowground communities can be hyperdiverse over relatively small spatial scales16, and even Darwin recognized that these communities can be diverse and important19. However, while awareness of the immense biodiversity belowground is increasing, we are only beginning to elucidate the influence of belowground biodiversity on ecosystem functions1620. Using integrative measures of biodiversity and multifunctionality should increase our ability to predict how biodiversity across taxa, both above- and belowground, shapes the suite of functions and services that ecosystems provide1021.

While biodiversity can influence ecosystem function, climate can directly and indirectly affect it as well922 (see Supplementary Fig. 1). For example, climate directly influences ecosystem function by accelerating the activity and interactions among consumers, detritivores, decomposers and microbes122324, and it indirectly affects ecosystem function by altering the composition of communities16. A growing number of experimental studies have crossed climate manipulations with biodiversity treatments to examine the interplay of these two factors on ecosystem functions3825. Observational studies have taken advantage of natural climatic variation across multiple sites to explore how climate and aboveground (but rarely, if ever, belowground) biodiversity affect EMF. Together, these studies highlight the links among biodiversity and climate on EMF, and to some extent, geographic variation in the relationship between EMF and biodiversity1026. Yet, there are few studies that explore geographic variation in how above- and belowground biodiversity, climate and ecosystem function are related. Thus, to understand the individual and combined influences of biodiversity and climate on EMF, a multifaceted approach where the links between biodiversity, both above- and belowground, and EMF across extensive climatic gradients is needed. Such an approach is critical for understanding BEF relationships, documenting whether these effects are large enough to rival the effects of the many other drivers of EMF, and scaling experimental results to real ecosystems at larger spatial scales18. Moreover, such an approach has important relevance for society: linking biodiversity above- and belowground to EMF across climatic gradients can better inform predictions about the interacting effects of changing climate and the loss and gain of species on the functioning of ecosystems1118.

Here we examine how biodiversity above- (plant species richness) and belowground (for example, bacterial operational taxonomic units (OTUs), faunal richness, archaeal OTUs and arbuscular mycorrhiza (AM) fungal richness), and climate and soil factors (Supplementary Table 1) influenced EMF at 60 sites arrayed along an extensive climatic gradient in alpine grasslands on the Tibetan Plateau, China (Supplementary Fig. 2). We predicted that (1) soil and plant biodiversity would be positively correlated with EMF, but plant biodiversity would better predict EMF911; (2) the combined effects of above- and belowground biodiversity would be a stronger predictor of EMF than either above- or belowground biodiversity alone711; (3) climate and other factors that vary geographically would, at least partially, mediate the effects of above- and belowground biodiversity on EMF (Supplementary Fig. 2)9202728. Our results indicate that EMF is positively associated with the biodiversity of plant, soil bacteria and soil fauna, but not related to the biodiversity of soil archaea or AM fungi. Moreover, the combined effects of above- and belowground biodiversity explained more of the variation in EMF among sites than did either factor alone, and models that included climatic and soil factors as well as above- and belowground biodiversity accounted for up to 86% of the variation in EMF among sites.

Results

Relationships between biodiversity and EMF

To examine whether there was a significant linear relationship between biodiversity and each component of EMF (Supplementary Fig. 3–8), as well as whether there was a significant linear biodiversity–EMF relationship, we conducted ordinary least squares (OLS) regressions using site means (n=60) for each above- and belowground component of the ecosystem (Fig. 1). As predicted, the relationship between biodiversity and EMF varied among below- and aboveground components of the ecosystem. Biodiversity, both aboveground (plant species richness) and belowground (soil biodiversity index, see Estimating biodiversity in the Methods section), was positively correlated with EMF, but the belowground effect varied among components: some groups were positively related (bacteria and fauna), while others were not related (archaea and AM fungi) to EMF (Fig. 1a–d). Aboveground biodiversity (that is, plant species richness) explained more of the variation in EMF (42%) than did soil biodiversity (32%) (Fig. 1e,f).

Figure 1

Relationships between below- and aboveground biodiversity and ecosystem multifunctionality.

Relationships between soil bacterial (a) soil faunal (b) soil archaeal (c) soil AM fungal (d) plant (e) and soil (f) biodiversity and ecosystem multifunctionality (EMF). There were no significant linear relationships between EMF and soil archaeal biodiversity or soil AM fungal biodiversity. The red fitted lines are from OLS regression. Only significant fitted lines are displayed on the graphs. Shaded areas show 95% confidence interval of the fit.

The effects of above- and belowground biodiversity on EMF

Most previous BEF studies have considered only plant species richness, while largely ignoring biodiversity belowground. Here we examined the combined effects of above- and belowground biodiversity on EMF using general linear models (GLMs). We found that 45% of the variation in EMF was accounted for by above- and belowground biodiversity (Supplementary Table 2). Next, we fitted partial linear models to tease apart the relative effects of above- and belowground biodiversity on EMF (Fig. 2). Plant species richness accounted for 19% of the variation in EMF. Soil biodiversity, however, accounted for 4% of the variation in EMF, although this result was not statistically significant (P=0.107, partial linear regression). As predicted, the combined effects of above- and belowground biodiversity explained more of the variation (22%) in EMF than did either factor alone. Overall, plant species richness alone was the best single predictor of EMF (Figs 1e, 2).

Figure 2

The relative effects of below- and aboveground biodiversity on ecosystem multifunctionality.

The relative effects of soil biodiversity (a) and plant species richness (b) on ecosystem multifuncitonality (EMF).The red fitted line is from a partial linear regression. Only significant fitted lines are shown on the graphs. The results of variance partitioning for soil biodiversity and plant species richness are shown by the R2 on this graph. Shaded areas show 95% confidence interval of the fit.

The effects of biodiversity and abiotic factors on EMF

To determine further whether the observed effects of above- and belowground biodiversity influenced EMF to the extent that abiotic factors did, we conducted regression as well as correlation analyses with mean annual precipitation (MAP), mean annual temperature (MAT), soil moisture, soil pH and soil CaCO3 as predictors in the models (Supplementary Figs 9 and 10). Then, we derived seven components from a principal component analysis (PCA) with the seven predictors mentioned above (Supplementary Fig. 11). The full model, with the seven components, accounted for 86% of the variation in EMF, and the best-fit model with two components (PC1 and PC4) accounted for 84% of the variation in EMF. These results indicate that the components of PC1 and PC4 were the strongest drivers of EMF (Supplementary Table 3). In addition, we found that PC1 represented MAP, soil moisture, soil pH, soil CaCO3, aboveground biodiversity and belowground biodiversity. PC4 represented MAT, soil moisture and soil CaCO3 (Supplementary Fig. 11).

Next, we fitted GLMs to investigate the biotic and abiotic effects on EMF (Table 1). We found that soil moisture was the most important driver, accounting for 65% of the variation in EMF (Supplementary Fig. 9). In total, 86% of the variation in EMF was accounted for by a model that included soil pH, aboveground biodiversity, belowground biodiversity and MAT (Table 1). MAP and soil CaCO3 had no influence on the model predictions. However, we found that MAP explained 43%, soil pH explained 30% and soil CaCO3 explained 33% of the variation in EMF (Supplementary Fig. 9). MAP and soil moisture were positively, and soil pH and soil CaCO3 were negatively, correlated with EMF. Somewhat surprisingly, we did not detect any relationship between temperature and EMF (Supplementary Fig. 9).

Table 1

Summary of the general linear models (GLMs) for the effects of abiotic factors (soil moisture, soil pH and mean annual temperature) and biodiversity (plant species richness and soil biodiversity) on ecosystem multifunctionality (EMF).

Source	Estimate	SE	t-value	Significance Pr (>|t|)	MS	F-value	Significance Pr (>F)	VIF
Multiple R2 0.86; residual SE 0.265 on 54 df
SM	1.492	0.138	10.818	< 0.001	17.17	243.753	< 0.001	1.286
pH	−1.111	0.529	−2.100	0.04038	1.82	25.833	< 0.001	1.466
Plant SR	0.178	0.063	2.832	0.00648	2.59	36.735	< 0.001	2.343
Soil biodiversity	0.307	0.103	2.995	0.00414	0.64	9.048	0.003988	1.939
MAT	0.037	0.013	2.820	0.00671	0.56	7.951	0.006705	1.081

df, degree of freedom; MS, mean square; MAT, mean annual temperature; pH, soil pH; plant SR, plant species richness; SE, standard errors; SM, soil moisture; VIF, variance inflation factor.

The best-fit model of the GLMs can be expressed as: EMF=0.982+1.492 × SM−1.111 × pH+0.178 × plant SR+0.307 × soil biodiversity+0.037 × MAT.

Finally, we fitted a piecewise structural equation model (SEM) to infer the direct and indirect effects of climate (MAP and MAT), soil properties (soil moisture and soil pH), and above- and belowground biodiversity on EMF. Overall, two models with Fisher's C statistic P>0.05 and ΔAICc=0.224 were selected (Fig. 3). The SEM demonstrated that the influence of precipitation on EMF was mediated through soil moisture, aboveground biodiversity, belowground biodiversity and soil pH. Temperature directly (β=0.16, standardized coefficient), and indirectly through soil pH (β=−0.27, standardized coefficient), impacted EMF. Soil moisture was the most significant parameter (β=0.62, standardized coefficient) influencing EMF. Aboveground and belowground biodiversity were interactively and approximately equal in magnitude in their influence on each other as well as on EMF (Fig. 3). The strongest relationship observed in the SEM analysis was between MAP and belowground biodiversity (β=0.79, standardized coefficient; Fig. 3b). There was also a weaker and negative relationship between soil pH and EMF (β=−0.11, standardized coefficient) (Fig. 3).

Figure 3

Structural equation models of climate, soil and biodiversity as predictors of ecosystem multifuncitonality (EMF).

Solid red arrows represent positive paths (P<0.05, piecewise s.e.m.), solid black arrows represent negative paths (P<0.05, piecewise s.e.m.) and dotted grey arrows represent non-significant paths (P>0.05, piecewise s.e.m.). We report the path coefficients as standardized effect sizes. Overall fit of piecewise s.e.m. was evaluated using Shipley's test of d-separation: Fisher's C statistic (if P>0.05, then no paths are missing and the model is a good fit) and Akaike information criterion (AIC). MAP, mean annual precipitation; MAT, mean annual temperature; SM, soil moisture, plant SR, plant species richness.

Discussion

Though debate exists about the mechanisms by which biodiversity is related to ecosystem function, numerous studies have demonstrated that aboveground biodiversity (that is, plant species richness) tends to be positively correlated with ecosystem function2629. Recent studies have expanded these analyses and explored the relationships among biodiversity and multiple ecosystem functions. Similar to the classical BEF work, these studies tend to find that aboveground biodiversity is generally positively related to EMF9112130. Our study, across an extensive climatic gradient on the Tibetan Plateau, compliments the growing EMF literature by demonstrating that biodiversity and multiple ecosystem functions were positively correlated. However, our results are unique in that the combined effects of above- and belowground biodiversity accounted for a relatively large fraction (45%) of the observed variation in EMF across sites. Moreover, these relationships were strongly modulated by climatic variation. The best-fitting models accounted for 86% of the variation in EMF and included soil moisture, soil pH, above- and belowground biodiversity as well as climatic variables (MAT). Thus, our results indicate that the abiotic environment dictates, or at least modifies, the effects of biodiversity on the functioning of ecosystems.

Combining belowground biodiversity with measures of aboveground biodiversity increased the predictive power of biodiversity on EMF (Supplementary Table 2; Fig. 3). Plant species richness was always positively related to EMF. Belowground biodiversity alone, however, accounted for only 4% of the variation in EMF (Fig. 2). There are several explanations for the differences in above- and belowground impacts on EMF. First, the influence of belowground biodiversity on EMF varied among taxa (Fig. 1). Soil bacterial and soil faunal biodiversity were positively correlated with EMF, as was plant species richness, but diversities of archaea and AM fungi were not correlated with EMF. Thus, the effects of biodiversity of multi-trophic organisms and ecosystem functions3132 may mediate the positive effects of plant, bacterial and faunal biodiversity on EMF93334. Climatic and other environmental differences among sites, such as MAP, may also be key drivers in shaping the positive relationship between biodiversity and EMF as well as the differences between above- and belowground communities26. Such a mechanism was supported by the SEMs (Fig. 3) where the influence of precipitation on EMF was mediated through soil moisture, soil pH, plant species richness and soil biodiversity. Finally, soil bacteria and soil fauna are relatively broad taxonomic groupings with diverse associated traits and functions, while archaea and AM fungi are more narrow taxonomic groupings with a more limited suite of functions. Thus, increasing biodiversity of broad taxonomic groups of free-living soil organisms, such as bacteria, should increase the diversity of substrates that are decomposed in the soil and returned to the plant community for uptake35, which should increase EMF.

The relationship between EMF and the biodiversity of AM fungi and archaea may reflect differences in the life histories of these taxa. AM fungi are closely tied to plants and often provide limiting nutrients such as phosphorus to plants36. In the current study, there were no effects of AM fungal biodiversity on EMF. However, when we considered the effects of AM fungal biodiversity on single ecosystem functions, we found that some relationships were negative, some were positive and some were non-significant (Supplementary Fig. 4). For example, AM fungal biodiversity was negatively correlated with soil phosphorus concentration, which could offset the positive effects of AM fungal biodiversity on root biomass and phosphorus content in aboveground biomass. Thus, the average of multiple functions encompassed in the index of EMF might not be a good metric to reflect the multiple and often interacting functions performed by AM fungi or the complex interactions between plant hosts and AM fungi. Similar to AMF, archaea showed a different pattern with EMF than did soil bacteria; however, it is premature to speculate about their influence on EMF given that their potential ecological impact and response to climate is still being elucidated3637.

Our results underscore the obvious fact that the abiotic environment influences EMF. At least on the Tibetan Plateau (and probably elsewhere), abiotic factors, such as precipitation and soil moisture, are key drivers shaping the biodiversity–EMF relationship (Fig. 3; Table 1; Supplementary Fig. 9). While our estimates of soil moisture come from a single measurement and are not integrated across the entire season (see Supplementary Note 1 for an extensive discussion of soil moisture), they were collected at the same time and were comparable to one another among the sites. Further, in alpine grassland ecosystems, higher water availability is generally associated with higher precipitation, higher aboveground net primary productivity and higher nutrient availability263839, which could facilitate resource use by plants840 and ultimately support more species.

In contrast to soil moisture, MAT generally had no effect on among-site variation in ecosystem functions when we considered it in isolation of other potential factors (Supplementary Fig. 9a). This pattern may emerge because, in this ecosystem, diurnal variation in temperature is large and often approximates seasonal variation41, suggesting that the communities are adapted to extensive variation in temperature. Thus, our results contrast with similar work from dryland ecosystems9, and suggest that an increase in temperature will not affect EMF in these alpine grassland ecosystems. However, in combination with the other biotic (biodiversity both above- and belowground) and abiotic factors (MAP, soil moisture and soil pH) that were included in the GLM and SEM, MAT significantly increased the explanatory power of the models. Therefore, other environmental factors may covary with the temperature–EMF relationships in broad-scale studies across climates, soil types and plant communities. Again, teasing apart why temperature is related to EMF in some ecosystems but not others is an important challenge, given ongoing increases in global temperatures.

Soil CaCO3 was highly correlated with soil pH, and was a better predictor of EMF than was soil pH, yet soil pH emerged as an important driver of EMF. Grassland soils across northern China, including those on the Tibetan Plateau, have experienced significant acidification since the 1980s (ref. 42). Soil acidification should inhibit decomposition of soil organic C, which might enhance the ability of alpine grassland ecosystems to perform multiple functions. Soil pH has been highlighted in other studies as being an important predictor of belowground biodiversity43.

Given that soils sustain life aboveground, a better understanding of the relative and combined effects of above- and belowground biodiversity is needed to predict the potential consequence of biodiversity loss for the future maintenance of ecosystem functions and services20. We note, however, that our approach to characterizing biodiversity belowground provides only a rough estimate of the actual belowground biodiversity because the sequencing methods we used did not detect all taxa (Supplementary Fig. 12). Hence, while unlikely, the relative effect of soil biodiversity on EMF may have been different had we used another suite of molecular methods to characterize soil biodiversity. Future studies linking above- and belowground biodiversity to EMF should use higher resolution sequencing profiling methods as well as techniques that focus on the active community to characterize the active component of belowground biodiversity. Despite this caveat, our work nevertheless demonstrates that abiotic factors (climate and soil) mediate the effects of above- and belowground biodiversity on EMF, which will depend strongly on the climatic gradients considered.

Belowground biodiversity, as well as biodiversity aboveground, clearly influence EMF. Our works build on a growing body of work suggesting that the correlation between biodiversity and EMF may be a general pattern across natural ecosystems9. However, our work also demonstrates that the relationship between components of biodiversity, here above- and belowground components of ecosystems, and EMF can be positively or not related. Moreover, biodiversity–EMF relationships may be strongly mediated by climate and be contingent on the environmental constraints of different ecosystems. Obviously, experimental work would help disentangle these factors further, and such experiments are a clear next step in EMF–biodiversity research globally. As climates change and species are lost and gained from ecosystems, predicting how ecosystems will function in the future will require experiments and observations that link biodiversity above- and belowground to EMF.

Methods

Sampling

We sampled plant and soil communities at 60 study sites over an extensive area (>1,000,000 km2) in the northeastern and central Tibetan Plateau in Qinghai Province and Tibetan Autonomous Region, China (Supplementary Fig. 2) during the peak-growing season (July–August) of 2011. Our survey captured a substantial range of the vegetation types, soil classes and climatic conditions found in the alpine grasslands on the Tibetan Plateau. The sites represent the three main vegetation types: alpine meadow, alpine steppe and desert steppe, and the 11 soil types (Genetic Soil Classification of China44): brown pedocals, castanozems, chernozems, cold calcic soils, dark felty soils, felty soils, frigid calcic soils, frigid frozen soils, grey–brown desert soils, grey-cinnamon soils and meadow soils on the Plateau. We selected sites to minimize the potential effects of grazers and other disturbances on soil and plant community structure. For each of the sites, we compiled MAT and MAP from the National Meteorological Bureau of China database. Data were compiled by interpolating data of monthly mean temperature and monthly precipitation records (1951–2010) from 716 climate stations across China (http://cdc.cma.gov.cn). The sites ranged in elevation from 2,918 to 5,228 m (mean, 4,064 m), in MAT from −5.2 to 4.7 °C (mean, −0.17 °C), and in MAP from 66 to 560 mm (mean, 365 mm) (Supplementary Table 4).

We collected 180 soil samples (0–5 cm; 3 samples per site) from the 60 sites. At each site, we established a 100-m transect and randomly placed three plots (1 × 1 m2) on each transect, with the stipulation that the plots were at least 40-m apart. Within each plot, 5–7 soil cores (5 cm in diameter) were collected, bulked, and homogenized in the field. Soil samples for total soil carbon, nitrogen, phosphorus and CaCO3 analysis were air-dried, sieved (2-mm mesh) and ground to a fine powder using a ball mill. Soil samples for molecular and physiochemical (soil moisture, soil available nitrogen and soil pH) analysis were packed in polyethylene bags, immediately stored in portable refrigerator powered with a car battery and then stored in the lab at −20 °C until processing45. Subsamples were in the freezer for no >1 week before measurements.

Physiochemical data measurements

Total soil carbon and nitrogen were determined by combustion on a CHN elemental analyzer (2400 II CHN elemental analyzer, PerkinElmer, Boston, MA, USA). Total soil phosphorus was determined by the molybdenum blue method with a ultraviolet–visible spectrophotometer (UV-2550, Shimadzu, Kyoto, Japan). Soil CaCO3 was analysed volumetrically on ground subsamples using a Calcimeter (Eijkelkamp, Netherland). Soil organic carbon was calculated as the difference between total soil carbon and carbon bound in soil CaCO3. The density of soil organic carbon, soil nitrogen, soil phosphorus and soil CaCO3 were calculated in the top 0–5 cm. Soil moisture was measured gravimetrically after ∼10-h desiccation at 105 °C. Soil available nitrogen (sum of ammonium, nitrate and dissolved organic nitrogen) were determined using a TOC-TN analyzer (Shimadzu, Kyoto, Japan). Soil pH was determined in a 1:5 ratio of fresh soil to deionised water slurry on a pH meter (Thermo 0rion-868).

DNA extraction and pyrosequencing

The biodiversity of soil bacteria, AM fungi and archaea were assessed using a suite of molecular techniques. The nucleic acids from each plot (n=180) were extracted from 0.5 g of soil using a FastDNA Spin kit (Bio 101, Carlsbad, CA, USA), according to the manufacturer's instructions, and then stored at −40 °C. Extracted DNA was diluted to ∼25 ng μl−1 with distilled water and stored at −20 °C until use. For bacterial, archaeal and AM fungal amplification 2 μm diluted DNA samples from each plot were used as a template.

To assess bacterial communities, we amplified the V4–V5 hypervariable regions of bacterial 16S ribosomal RNAs (rRNAs; Escherichia coli positions 515–907) using the primer set: F515 (ref. 46): 5′-GTGCCAGCMGCCGCGG-3′ with the Roche 454 ‘A' pyrosequencing adapter and a unique 7-bp barcode sequence, and primer R907 (ref. 47): 5′-CCGTCAATTCMTTTRAGTTT-3′ with the Roche 454 `B' sequencing adapter at the 5′-end of each primer. We conducted PCR amplification with 25-μl 2 × premix (TaKaRa), 0.5-μl 20 mM each forward and reverse primer and 50 ng of DNA, and the volume was completed to 50 μl with double-distilled water. Each sample was amplified in triplicate with the 50-μl reaction under the following conditions: 30 cycles of denaturation at 94 °C for 30 s, annealing at 55 °C for 30 s and extension at 72 °C for 30 s; with a final extension at 72 °C for 10 min43. PCR products from each sample were pooled together, purified with an agarose gel DNA purification kit (TaKaRa), quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, USA), then combined in equimolar ratios in a single tube and run on a Roche FLX454 pyrosequencing machine (Roche Diagnostics Corp., Branford, CT, USA), producing reads from the forward direction F515.

To assess archaeal communities, the V3–V5 hypervariable regions of archaeal 16S rRNA48 were amplified using the primer set: Arch344F: 5′-ACGGGGYGCAGCAGGCGCGA-3′ with the Roche 454 ‘A' pyrosequencing adapter and a unique 7-bp barcode sequence, and primer Arch915R: 5′-GTGCTCCCCCGCCAATTCCT-3′ with the Roche 454 ‘B' sequencing adapter at the 5′-end of each primer. We conducted PCR amplification with 25-μl 2 × premix (TaKaRa), 0.5-μl 20 mM each forward and reverse primer and 50 ng of DNA, and the volume was completed to 50-μl with double-distilled water. Each sample was amplified in triplicate with the 50-μl reaction under the following conditions: 94 °C for 5 min, 10 cycles of touchdown PCR (denaturation at 94 °C for 30 s, annealing for 30 s with a 0.5 °C per cycle decrement at 61 °C above the respective annealing temperatures and extension at 72 °C for 1 min), followed by 25 cycles of regular PCR (94 °C for 30 s, 30 s at the respective annealing temperature, 72 °C for 1 min and a final extension step for 7 min at 72 °C)48. PCR products from each sample were pooled together and purified with an Agarose Gel DNA purification kit (TaKaRa), quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, USA), then combined in equimolar ratios in a single tube and run on the Roche FLX454 pyrosequencing machine, producing reads from the forward direction Arch344F.

To assess AM fungal communities, the 18S rRNA gene fragment for the 454 GS-FLX pyrosequencing platform was amplified using the primer set: AMV4.5NF (5′-AAGCTCGTAGTTGAATTTCG-3′) with the Roche 454 ‘A' pyrosequencing adapter and a unique 7-bp barcode, and primer AMDGR (5′-CCCAACTATCCCTATTAATCAT-3′) with the Roche 454 ‘B' sequencing adapter at the 5′-end of each primer49. We conducted PCR amplification with 25-μl 2 × premix (TaKaRa), 0.5-μl 20 mM each forward and reverse primer and 50 ng of DNA, and the volume was completed to 50-μl with double-distilled water. Each sample was amplified in triplicate with the 50- μl reaction mixtures under the following conditions: 30 cycles of denaturation at 94 °C for 30 s, annealing at 58 °C for 30 s and extension at 72 °C for 30 s; with a final extension at 72 °C for 10 min50. PCR products from each sample were pooled together, purified with a Agarose Gel DNA purification kit (TaKaRa), quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, USA) and then combined in equimolar ratios in a single tube. 454-sequencing was performed by BGI (Shenzhen, China) on a Roche FLX454 pyrosequencing machine, producing reads from the forward direction AMV4.5NF.

Pyrosequencing data analyses

We processed and analysed bacterial and archaeal biodiversity data as described by Chu et al.45 and Hamady et al.51 with the Quantitative Insights into Microbial Ecology pipeline (http://qiime.sourceforge.net/)52. In brief, sequences >200-bp long with an average quality score >25 and no ambiguous characters were included in the analyses. Bacterial and archaeal phylotypes were identified using Uclust53 and assigned to OTUs based on 97% similarity. A representative sequence was chosen from each phylotype by selecting the most highly connected sequence51. All representative sequences were aligned by Py NAST54. The taxonomic identity (Supplementary Table 5) of each phylotype was determined using the Greengenes database (http://greengenes.lbl.gov/). Filtering of the sequences to remove erroneous, OTUs due to sequence errors and chimeras was conducted using the USEARCH tool in QIIME, version 1.8.0. Because we relied on filtering rather than denoising of the data, we also deleted singletons and set a threshold for a high-quality score (that is, 30) when running the command split_libraries.py, similar to approaches from other studies (for example, refs 51, 55, 56). Total sequence counts were normalized to 4,000 and 2,000 sequences per sample for bacterial and archaeal samples, respectively (Supplementary Fig. 12). AM fungal sequences were processed and analysed using Mothur v.1.30.1 (ref. 57). The fasta, quality and flow files were extracted from the standard flowgram format file of Roche, and were trimmed and de-noised (pdiffs=0, bdiffs=0, minflows=360, maxflows=550, maxhomop=8, minlength=230 and flip=T). Remaining sequences >230 bp in length49 were aligned to the reference SILVA database58. Sequences not aligning within the optimized alignment region were removed from the analysis with the screening function. The remaining sequencing were processed to reduce sequencing noise by a pre-clustering methodology59, and to remove chimeric sequences using UCHIME60. The OTUs were assigned by an average neighbour algorithm with a 0.03 dissimilarity cutoff, based on the sequences and/or OTUs obtained within the phylum of Glomeromycota. In spite of the specificity of the AM fungal primer, Alveolata, Metazoa, Viridiplantae, Stramenopiles, Dikarya and unclassified consensus taxonomy were detected; they were removed as environmental contaminants. AM fungal OTUs corresponding to unique AM virtual taxa were defined on the basis of a BLAST search of representative sequences from each OTU against the MaarjAM database61 (Supplementary Table 5; Supplementary Fig. 12).

Soil fauna extraction and identification

We assessed soil faunal biodiversity by collecting and homogenizing three soil cores within each plot (3.5 cm in diameter, 0–15 cm in depth). Cores were returned to the lab and soil fauna were extracted using a modified Berlese–Tullgren apparatus62. To be more specific, soil samples for microarthopods who prefer dry environments were extracted through Tullgren funnels for 48 h (dry funnel method). Soil samples for nematodes and enchytraeids, who prefer wet environments, were wrapped in nylon cloth and extracted through Berlese funnels for 48 h (wet funnel method). They were counted, identified at order level, and preserved in 75% ethyl alcohol.

Plant biomass measurement and identification

Before soil sampling, the plant communities were surveyed and harvested in each 1-m2 plot to measure standing aboveground biomass. The plant material was dried at 60 °C for ∼12 h and weighed. Then, we ground the aboveground plant material to a fine powder on a ball mill and analysed for plant nitrogen, and phosphorus using the techniques outlined for soils above. We estimated root biomass by collecting 4–8 soil cores per 1-m2 plot (0–5-cm deep, 7-cm diameter). A ratio of live:dead (56%)63 was used to calibrate root biomass. Plant species richness was assessed at three plots (1 × 1 m2) located 10-m away from the soil samples. We listed all vascular plant species for each of the 180 plots.

Quantifying ecosystem multifunctionality

In total, we collected data that quantify key ecosystem functions and related variables: (1) aboveground biomass, (2) root biomass, (3) soil organic carbon, (4) soil nitrogen, (5) soil available nitrogen, (6) soil phosphorus, (7) plant nitrogen (nitrogen pools in aboveground biomass), and (8) plant phosphorus (phosphorus pools in aboveground biomass).

There are several ways to estimate the relationship between biodiversity and EMF21: the single functions approach, turnover approach, averaging approach and threshold approach. Because each approach has its strengths and weaknesses21, we used three distinct methods to calculate multifunctionality: single functions approach, an averaging approach, and a multiple thresholds approach. The relationships between EMF and biodiversity using the multiple thresholds approach were very similar to that obtained with the averaging approach (hereafter, EMF index), and hence we used the EMF index as described by Maestre et al.9 in the main text and in further analyses. We acknowledge that the EMF index, calculated from functions that are often correlated with one other (Supplementary Figs 13 and 14), does not provide a single best metric of EMF21. We use the EMF index for two reasons: (1) we are interested in biological or service-based outcomes21, and (2) this method is a straightforward and interpretable measure of a community's ability to sustain multiple functions simultaneously9. To obtain the EMF index for each site, we calculated the Z-scores for each of the eight ecosystem functions evaluated. The Z-scores for the measured variables were averaged to obtain an EMF for each site.

Estimating biodiversity

We related EMF to biodiversity for each site, where biodiversity is the complete tally of all species at each sampled site. More specifically, plant biodiversity (species richness) is the number of plant species, bacterial and archaeal biodiversity is the number of OTUs, AM fungal richness is the number of virtual taxa, and faunal richness is the number of taxonomic units at the order level tallied at each site. We were conservative in our estimates and used a similarity of 97% (ref. 64) to delineate OTUs to represent soil bacteria and soil archaea species. We used virtual taxa for AM fungi, which are defined with bootstrap support and a sequence similarity of ≥97% (ref. 65). Note that we use the term ‘biodiversity' throughout, although the microbial and archaeal groups are phylotypes, AM fungal groups are virtual taxa, faunal groups are taxa at the order level. In addition, we used an approach of averaging standardized values of belowground biodiversity into a single index (soil biodiversity) as described by Wagg et al.11. An index was used from the average of all standardized biodiversity indices (bacterial OTUs, archaeal OTUs, faunal richness and AM fungal richness) to represent changes in soil biodiversity.

Statistical analyses

The soil characteristics, ecosystem functions and above- and belowground biodiversity were assessed for a total of 180 plots. However, we use site means for the data analyses, except for the biodiversity assessments that were based on number of species or OTUs present in each site (n=60). First, we evaluated the relationships between biodiversity and each component of EMF (Supplementary Figs 3–8), as well as the relationships between biodiversity and EMF relationship (Fig. 1) using OLS regressions. Then, we evaluated the combined effects of above- and belowground biodiversity on EMF by fitting GLMs (Supplementary Table 2). And, we fitted partial linear regression as described in Quinn and Keough66 to identify the relative biodiversity effects of soil and plant biodiversity on EMF (Fig. 2). The partial linear regression allowed us to estimate how much of the variation in EMF can be exclusively attributed to soil and plant biodiversity. To further determine whether the observed effects of soil and plant biodiversity influenced EMF as much as abiotic factors, we conducted both regression and correlation analyses for seven predictors including MAP, MAT, soil moisture, soil pH, soil CaCO3, soil biodiversity and plant species richness (Supplementary Figs 9 and 10). Then, we fitted GLMs to investigate the biotic and abiotic effects on EMF (Table 1). To address multicollinearity, we first conducted a PCA using a correlation matrix with the seven predictors mentioned above. Seven components of this PCA were retained for further analyses. Two components that had eigenvalues >1 (PC1=1.871, PC2=1.120), explaining 68% variation in the data. We fitted two models of GLMs, one is a full model with the seven components, the other one is a reduced model with two components (PC1 and PC2) (Supplementary Fig. 11). Then, we used variables with VIF (variance inflation factor)<3 to select variables for the GLMs analyses (Table 1).

In addition, because the partial linear regression and the GLMs assume that the biodiversity and/or abiotic effects are additive, but not interactive, we further fitted a piecewise SEM (piecewise SEM)67 to infer relative importance of climate (MAP and MAT), soil (soil moisture and soil pH) and biodiversity (soil and plant) on EMF. Compared with the traditional variance covariance-based SEM, the piecewise SEM could (1) piece multiple separate (generalized) linear models together to a single causal network, (2) use Shipley's test of d-separation to test whether any paths are missing from the model and (3) use Akaike information criterion (AIC) to compare nested models and for small sample size (http://jonlefcheck.net/2014/07/06/piecewise-structural-equation-modeling-in-ecological-research/). We constructed the piecewise SEM based on the schematic diagram in the Supplementary Fig. 1. Because piecewise SEM is only for recursive models, we selected two models for the piecewise SEM; one assumed that plant species richness drives belowground biodiversity, the other one assumed that belowground biodiversity drives plant species richness (Fig. 3). We fitted the component models of the piecewise SEM as linear models. We reported the standardized coefficient for each path from each component models. Overall fit of the piecewise SEM was evaluated using Shipley's test of d-separation: Fisher's C statistic68 and AIC69 in the R package ‘piecewiseSEM' (https://github.com/jslefche/piecewiseSEM). MAP, soil pH, aboveground biomass, belowground biomass, soil organic carbon, soil nitrogen, soil available nitrogen, soil phosphorus, plant nitrogen and plant phosphorus were log-transformed. Soil moisture, soil CaCO3, bacterial, faunal, archaeal, AM fungal and plant biodiversity were sqrt-transformed. All statistical analyses were performed using R version 3.0.2 (R Foundation for Statistical Computing, Vienna, Austria, 2013).

Additional information

Accession codes: The 454 pyrosequencing data set of soil bacteria are deposited in the DDBJ Sequence Read Archive (http://trace.ddbj.nig.ac.jp/DRASearch) with accession number: DRA001226, the 454 pyrosequencing data set of soil archaea are deposited in the European Nucleotide Archive (http://www.ebi.ac.uk/ena) with accession number: PRJEB8007, and the 454 pyrosequencing data set of soil AM fungi are deposited in the NCBI Bio Project Archive (http://www.ncbi.nlm.nih.gov/bioproject/269515) with accession number: PRJNA269515.

How to cite this article: Jing, X. et al. The links between ecosystem multifunctionality and above- and belowground biodiversity are mediated by climate. Nat. Commun. 6:8159 doi: 10.1038/ncomms9159 (2015).