The links between ecosystem multifunctionality and above- and belowground biodiversity are mediated by climate. — TEI-tagged view

Xin JingpersName, Nathan J. SanderspersName, Yu ShipersName, Haiyan ChupersName, Aimée T. ClassenpersName, Ke ZhaopersName, Litong ChenpersName, Yue ShipersName, Youxu JiangpersName, Jin-Sheng HepersName (2015date). “The links between ecosystem multifunctionality and above- and belowground biodiversity are mediated by climate..” Nature Communicationstitle·j.

Object of study — ecosystem multifunctionalityterm·key and above- and belowground biodiversityterm·dom across 60 Tibetan Plateauterm·key grassland sites

Research question — how above- and belowground biodiversityterm·dom jointly affect EMFterm·key, and whether climateterm·dom mediates it

Contribution — biotic and abiotic variables explain up to 86% of EMFterm·key variation; above- and belowground biodiversityterm·dom together 45%

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

The links between ecosystem multifunctionalityterm·key and above- and belowground biodiversityterm·dom are mediated by climateterm·dom.

Abstract 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 biodiversityterm·dom on multiple ecosystem functionsterm·dom (for example, ecosystem multifunctionalityterm·key, EMFterm·key) or how climateterm·dom might mediate those relationships. Here we tease apart the effects of biotic and abiotic factors, both above- and belowground, on EMFterm·key on the Tibetan Plateauterm·key, China. We found that a suite of biotic and abiotic variables account for up to 86% of the variation in EMFterm·key, with the combined effects of above- and belowground biodiversityterm·dom accounting for 45% of the variation in EMFterm·key. Our results have two important implications: first, including belowground biodiversityterm·dom in models can improve the ability to explain and predict EMFterm·key. Second, regional-scale variation in climateterm·dom, and perhaps climateterm·dom change, can determine, or at least modify, the effects of biodiversity on EMFterm·key in natural ecosystems.

Ecologists have explored the links between biodiversity and ecosystem functioning for more than two decades 1 2 3 4 . Yet, most of our understanding of biodiversity–ecosystem functionterm·dom (BEFterm·dom) relationships comes from experimental or observational investigations of how a single ecosystem functionterm·dom, usually aboveground net primary productivity, depends on plant biodiversity in small plots at single sites 2 4 5 (but see refs 3 , 6 for multisite comparisons). Clearly, ecosystems perform multiple functions simultaneously (ecosystem multifunctionalityterm·key; EMFterm·key) 7 8 9 10 11 , and these functions can be mediated, if not controlled, by biodiversity in other components of ecosystems, such as the biodiversity of taxa belowground 11 12 13 14 15 16 17 . BEFterm·dom relationships might also vary geographically, perhaps because of spatial or temporal variation in climateterm·dom or other abiotic factors 3 5 9 18 .

Although above- and belowground communities and their interactions can shape multiple ecosystem functionsterm·dom, it is difficult to tease apart their relative effects 16 , perhaps because taxa operate at different spatial scales. For instance, belowground communities can be hyperdiverse over relatively small spatial scales 16 , and even Darwin recognized that these communities can be diverse and important 19 . However, while awareness of the immense biodiversity belowground is increasing, we are only beginning to elucidate the influence of belowground biodiversityterm·dom on ecosystem functionsterm·dom 16 20 . 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 provide 10 21 .

While biodiversity can influence ecosystem functionterm·dom, climateterm·dom can directly and indirectly affect it as well 9 22 (see Supplementary Fig. 1 ). For example, climateterm·dom directly influences ecosystem functionterm·dom by accelerating the activity and interactions among consumers, detritivores, decomposers and microbes 12 23 24 , and it indirectly affects ecosystem functionterm·dom by altering the composition of communities 16 . A growing number of experimental studies have crossed climateterm·dom manipulations with biodiversity treatments to examine the interplay of these two factors on ecosystem functionsterm·dom 3 8 25 . Observational studies have taken advantage of natural climatic variation across multiple sites to explore how climateterm·dom and aboveground (but rarely, if ever, belowground) biodiversity affect EMFterm·key. Together, these studies highlight the links among biodiversity and climateterm·dom on EMFterm·key, and to some extent, geographic variation in the relationship between EMFterm·key and biodiversity 10 26 . Yet, there are few studies that explore geographic variation in how above- and belowground biodiversityterm·dom, climateterm·dom and ecosystem functionterm·dom are related. Thus, to understand the individual and combined influences of biodiversity and climateterm·dom on EMFterm·key, a multifaceted approach where the links between biodiversity, both above- and belowground, and EMFterm·key across extensive climatic gradientsterm·key is needed. Such an approach is critical for understanding BEFterm·dom relationships, documenting whether these effects are large enough to rival the effects of the many other drivers of EMFterm·key, and scaling experimental results to real ecosystems at larger spatial scales 18 . Moreover, such an approach has important relevance for society: linking biodiversity above- and belowground to EMFterm·key across climatic gradientsterm·key can better inform predictions about the interacting effects of changing climateterm·dom and the loss and gain of species on the functioning of ecosystems 11 18 .

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

Results Relationships between biodiversity and EMFterm·key To examine whether there was a significant linear relationship between biodiversity and each component of EMFterm·key ( Supplementary Fig. 3–8 ), as well as whether there was a significant linear biodiversity–EMFterm·key 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 EMFterm·key varied among below- and aboveground components of the ecosystem. Biodiversity, both aboveground (plant species richnessterm·key) and belowground (soil biodiversity index, see Estimating biodiversity in the Methods section), was positively correlated with EMFterm·key, but the belowground effect varied among components: some groups were positively related (bacteria and fauna), while others were not related (archaea and AM fungiterm·key) to EMFterm·key ( Fig. 1a–d ). Aboveground biodiversity (that is, plant species richnessterm·key) explained more of the variation in EMFterm·key (42%) than did soil biodiversity (32%) ( Fig. 1e,f ).

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

The effects of biodiversity and abiotic factors on EMFterm·key To determine further whether the observed effects of above- and belowground biodiversityterm·dom influenced EMFterm·key to the extent that abiotic factors did, we conducted regression as well as correlation analyses with mean annual precipitationterm·dom (MAPterm·dom), mean annual temperatureterm·dom (MATterm·dom), soil moistureterm·dom, soil pH and soil CaCO 3 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 EMFterm·key, and the best-fit model with two components (PC1 and PC4) accounted for 84% of the variation in EMFterm·key. These results indicate that the components of PC1 and PC4 were the strongest drivers of EMFterm·key ( Supplementary Table 3 ). In addition, we found that PC1 represented MAPterm·dom, soil moistureterm·dom, soil pH, soil CaCO 3 , aboveground biodiversity and belowground biodiversityterm·dom. PC4 represented MATterm·dom, soil moistureterm·dom and soil CaCO 3 ( Supplementary Fig. 11 ).

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

Finally, we fitted a piecewise structural equation modelterm·key (SEM) to infer the direct and indirect effects of climateterm·dom (MAPterm·dom and MATterm·dom), soil properties (soil moistureterm·dom and soil pH), and above- and belowground biodiversityterm·dom on EMFterm·key. 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 EMFterm·key was mediated through soil moistureterm·dom, aboveground biodiversity, belowground biodiversityterm·dom and soil pH. Temperature directly ( β =0.16, standardized coefficient), and indirectly through soil pH ( β =−0.27, standardized coefficient), impacted EMFterm·key. Soil moistureterm·dom was the most significant parameter ( β =0.62, standardized coefficient) influencing EMFterm·key. Aboveground and belowground biodiversityterm·dom were interactively and approximately equal in magnitude in their influence on each other as well as on EMFterm·key ( Fig. 3 ). The strongest relationship observed in the SEM analysis was between MAPterm·dom and belowground biodiversityterm·dom ( β =0.79, standardized coefficient; Fig. 3b ). There was also a weaker and negative relationship between soil pH and EMFterm·key ( β =−0.11, standardized coefficient) ( Fig. 3 ).

Discussion Though debate exists about the mechanisms by which biodiversity is related to ecosystem functionterm·dom, numerous studies have demonstrated that aboveground biodiversity (that is, plant species richnessterm·key) tends to be positively correlated with ecosystem functionterm·dom 26 29 . Recent studies have expanded these analyses and explored the relationships among biodiversity and multiple ecosystem functionsterm·dom. Similar to the classical BEFterm·dom work, these studies tend to find that aboveground biodiversity is generally positively related to EMFterm·key 9 11 21 30 . Our study, across an extensive climatic gradientterm·key on the Tibetan Plateauterm·key, compliments the growing EMFterm·key literature by demonstrating that biodiversity and multiple ecosystem functionsterm·dom were positively correlated. However, our results are unique in that the combined effects of above- and belowground biodiversityterm·dom accounted for a relatively large fraction (45%) of the observed variation in EMFterm·key across sites. Moreover, these relationships were strongly modulated by climatic variation. The best-fitting models accounted for 86% of the variation in EMFterm·key and included soil moistureterm·dom, soil pH, above- and belowground biodiversityterm·dom as well as climatic variables (MATterm·dom). Thus, our results indicate that the abiotic environment dictates, or at least modifies, the effects of biodiversity on the functioning of ecosystems.

Combining belowground biodiversityterm·dom with measures of aboveground biodiversity increased the predictive powerterm·mth of biodiversity on EMFterm·key ( Supplementary Table 2 ; Fig. 3 ). Plant species richnessterm·key was always positively related to EMFterm·key. Belowground biodiversityterm·dom alone, however, accounted for only 4% of the variation in EMFterm·key ( Fig. 2 ). There are several explanations for the differences in above- and belowground impacts on EMFterm·key. First, the influence of belowground biodiversityterm·dom on EMFterm·key varied among taxa ( Fig. 1 ). Soil bacterial and soil faunal biodiversity were positively correlated with EMFterm·key, as was plant species richnessterm·key, but diversities of archaea and AM fungiterm·key were not correlated with EMFterm·key. Thus, the effects of biodiversity of multi-trophic organisms and ecosystem functionsterm·dom 31 32 may mediate the positive effects of plant, bacterial and faunal biodiversity on EMFterm·key 9 33 34 . Climatic and other environmental differences among sites, such as MAPterm·dom, may also be key drivers in shaping the positive relationship between biodiversity and EMFterm·key as well as the differences between above- and belowground communities 26 . Such a mechanism was supported by the SEMs ( Fig. 3 ) where the influence of precipitation on EMFterm·key was mediated through soil moistureterm·dom, soil pH, plant species richnessterm·key and soil biodiversity. Finally, soil bacteria and soil fauna are relatively broad taxonomic groupings with diverse associated traits and functions, while archaea and AM fungiterm·key 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 uptake 35 , which should increase EMFterm·key.

The relationship between EMFterm·key and the biodiversity of AM fungiterm·key and archaea may reflect differences in the life histories of these taxa. AM fungiterm·key are closely tied to plants and often provide limiting nutrients such as phosphorus to plants 36 . In the current study, there were no effects of AM fungal biodiversity on EMFterm·key. However, when we considered the effects of AM fungal biodiversity on single ecosystem functionsterm·dom, 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 EMFterm·key might not be a good metric to reflect the multiple and often interacting functions performed by AM fungiterm·key or the complex interactions between plant hosts and AM fungiterm·key. Similar to AMF, archaea showed a different pattern with EMFterm·key than did soil bacteria; however, it is premature to speculate about their influence on EMFterm·key given that their potential ecological impact and response to climateterm·dom is still being elucidated 36 37 .

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

In contrast to soil moistureterm·dom, MATterm·dom generally had no effect on among-site variation in ecosystem functionsterm·dom 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 variation 41 , suggesting that the communities are adapted to extensive variation in temperature. Thus, our results contrast with similar work from dryland ecosystems 9 , and suggest that an increase in temperature will not affect EMFterm·key in these alpine grasslandterm·key ecosystems. However, in combination with the other biotic (biodiversity both above- and belowground) and abiotic factors (MAPterm·dom, soil moistureterm·dom and soil pH) that were included in the GLM and SEM, MATterm·dom significantly increased the explanatory power of the models. Therefore, other environmental factors may covary with the temperature–EMFterm·key relationships in broad-scale studies across climatesterm·dom, soil types and plant communities. Again, teasing apart why temperature is related to EMFterm·key in some ecosystems but not others is an important challenge, given ongoing increases in global temperatures.

Soil CaCO 3 was highly correlated with soil pH, and was a better predictor of EMFterm·key than was soil pH, yet soil pH emerged as an important driver of EMFterm·key. Grassland soils across northern China, including those on the Tibetan Plateauterm·key, 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 grasslandterm·key ecosystems to perform multiple functions. Soil pH has been highlighted in other studies as being an important predictor of belowground biodiversityterm·dom 43 .

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

Belowground biodiversityterm·dom, as well as biodiversity aboveground, clearly influence EMFterm·key. Our works build on a growing body of work suggesting that the correlation between biodiversity and EMFterm·key may be a general pattern across natural ecosystems 9 . However, our work also demonstrates that the relationship between components of biodiversity, here above- and belowground components of ecosystems, and EMFterm·key can be positively or not related. Moreover, biodiversity–EMFterm·key relationships may be strongly mediated by climateterm·dom 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 EMFterm·key–biodiversity research globally. As climatesterm·dom 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 EMFterm·key.

Methods Sampling We sampled plant and soil communities at 60 study sites over an extensive area (>1,000,000 km 2 ) in the northeastern and central Tibetan Plateauterm·key 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 grasslandsterm·key on the Tibetan Plateauterm·key. The sites represent the three main vegetation types: alpine meadow, alpine steppe and desert steppe, and the 11 soil types (Genetic Soil Classification of China 44 ): 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 MATterm·dom and MAPterm·dom 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 climateterm·dom 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 MATterm·dom from −5.2 to 4.7 °C (mean, −0.17 °C), and in MAPterm·dom 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 m 2 ) 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 CaCO 3 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 moistureterm·dom, 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 processing 45 . 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 CaCO 3 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 CaCO 3 . The density of soil organic carbon, soil nitrogen, soil phosphorus and soil CaCO 3 were calculated in the top 0–5 cm. Soil moistureterm·dom 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 fungiterm·key 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 min 43 . 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 rRNA 48 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 primer 49 . 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 min 50 . 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 ChupersName 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 Uclust 53 and assigned to OTUs based on 97% similarity. A representative sequence was chosen from each phylotype by selecting the most highly connected sequence 51 . All representative sequences were aligned by Py NAST 54 . 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 length 49 were aligned to the reference SILVA database 58 . 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 methodology 59 , and to remove chimeric sequences using UCHIME 60 . 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 database 61 ( 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 apparatus 62 . 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-m 2 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-m 2 plot (0–5-cm deep, 7-cm diameter). A ratio of live:dead (56%) 63 was used to calibrate root biomass. Plant species richnessterm·key was assessed at three plots (1 × 1 m 2 ) located 10-m away from the soil samples. We listed all vascular plant species for each of the 180 plots.

Quantifying ecosystem multifunctionalityterm·key In total, we collected data that quantify key ecosystem functionsterm·dom 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 EMFterm·key 21 : the single functions approach, turnover approach, averaging approach and threshold approach. Because each approach has its strengths and weaknesses 21 , we used three distinct methods to calculate multifunctionality: single functions approach, an averaging approach, and a multiple thresholds approach. The relationships between EMFterm·key and biodiversity using the multiple thresholds approach were very similar to that obtained with the averaging approach (hereafter, EMFterm·key index), and hence we used the EMFterm·key index as described by Maestre et al . 9 in the main text and in further analyses. We acknowledge that the EMFterm·key index, calculated from functions that are often correlated with one other ( Supplementary Figs 13 and 14 ), does not provide a single best metric of EMFterm·key 21 . We use the EMFterm·key index for two reasons: (1) we are interested in biological or service-based outcomes 21 , and (2) this method is a straightforward and interpretable measure of a community's ability to sustain multiple functions simultaneously 9 . To obtain the EMFterm·key index for each site, we calculated the Z-scores for each of the eight ecosystem functionsterm·dom evaluated. The Z-scores for the measured variables were averaged to obtain an EMFterm·key for each site.

Estimating biodiversity We related EMFterm·key 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 fungiterm·key, 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 biodiversityterm·dom 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 functionsterm·dom and above- and belowground biodiversityterm·dom 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 EMFterm·key ( Supplementary Figs 3–8 ), as well as the relationships between biodiversity and EMFterm·key relationship ( Fig. 1 ) using OLS regressions. Then, we evaluated the combined effects of above- and belowground biodiversityterm·dom on EMFterm·key by fitting GLMs ( Supplementary Table 2 ). And, we fitted partial linear regression as described in Quinn and Keough 66 to identify the relative biodiversity effects of soil and plant biodiversity on EMFterm·key ( Fig. 2 ). The partial linear regression allowed us to estimate how much of the variation in EMFterm·key can be exclusively attributed to soil and plant biodiversity. To further determine whether the observed effects of soil and plant biodiversity influenced EMFterm·key as much as abiotic factors, we conducted both regression and correlation analyses for seven predictors including MAPterm·dom, MATterm·dom, soil moistureterm·dom, soil pH, soil CaCO 3 , soil biodiversity and plant species richnessterm·key ( Supplementary Figs 9 and 10 ). Then, we fitted GLMs to investigate the biotic and abiotic effects on EMFterm·key ( 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 climateterm·dom (MAPterm·dom and MATterm·dom), soil (soil moistureterm·dom and soil pH) and biodiversity (soil and plant) on EMFterm·key. 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 richnessterm·key drives belowground biodiversityterm·dom, the other one assumed that belowground biodiversityterm·dom drives plant species richnessterm·key ( 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 statistic 68 and AIC 69 in the R package ‘piecewiseSEM' ( https://github.com/jslefche/piecewiseSEM ). MAPterm·dom, 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 moistureterm·dom, soil CaCO 3 , 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 fungiterm·key 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: JingpersName, X. et al . The links between ecosystem multifunctionalityterm·key and above- and belowground biodiversityterm·dom are mediated by climateterm·dom. Nat. Commun. 6:8159 doi: 10.1038/ncomms9159 (2015).

Supplementary Material Supplementary Figures 1-14, Supplementary Tables 1-5, Supplementary Note and Supplementary References