Soil Food Web Complexity Enhances the Link Between Soil Biodiversity and Ecosystem Multifunctionality in Agricultural Systems

Background: Belowground biodiversity supports multiple ecosystem functions and services that humans rely on. However, there is a dearth of studies conducted on a large spatial scale on the topic in intensely managed agricultural ecosystems. Existing studies have overlooked the fact that the functional diversity in other trophic groups within a food web could inuence function of an individual in another trophic group. Here, we report signicant and positive relationships between soil biodiversity (archaea, bacteria, fungi, protists, and invertebrates) and multiple ecosystem functions (nutrient provisioning, element cycling, and reduced pathogenicity potential) in 228 agricultural elds. Results: The relationships were inuenced by (I) the types of organisms with signicant relationships in archaea, bacteria, and fungi and not in protists and invertebrates, and (II) the connectedness of dominant phylotypes across soil food webs, which generate different ecological clusters within soil networks to maintain multiple functions. In addition, we highlight the role of soil food web complexity, reected by ecological networks comprising diverse organisms, which promote the multiple functions and enhance the link between soil biodiversity and ecosystem functions. Conclusions: Overall, our results represent a signicant advance in forecasting the impacts of belowground biodiversity within food webs on ecosystem functions in agricultural systems, and suggest that soil biodiversity, particularly soil food web complexity, should not be overlooked, but rather considered a key factor and integrated into policy and management activities aimed at enhancing and maintaining ecosystem productivity, stability, and sustainability under land-use intensication.

Such studies have focused on the biodiversity of a limited number of soil organism types and in natural terrestrial ecosystems. However, our understanding of the relationship between such multidiversity, including different groups of soil microbes (archaea, bacteria, fungi, protists) and invertebrates (nematodes, annelids, and arthropods), and multiple ecosystem function under human-managed conditions over large spatial scales, such agricultural elds, remains poor. Considering the importance of soil biodiversity for global food supply, understanding how soil biodiversity regulates multifunctionality in agricultural ecosystems is critical for the management of soil communities, which could be exploited for better provisioning of key ecosystem services that support sustainable crop yields.
Trophic interactions play important roles in stimulating ecosystem processes [1]. Soil organisms are structured, and form complex food webs, re ecting strong interrelationships within ecological networks [2,15]. Soil food web characteristics, such as architecture and connectedness, could strongly predict and in uence energy ux, food web stability, and carbon (C) and nitrogen (N) cycling processes [1,15,16]. In addition, the phylotypes within soil food webs with similar environmental and resource preferences, could form strongly connected ecological clusters in ecological networks, with major implications for ecosystem functioning [2,17]. Soil physico-chemical properties and ecosystem processes may vary between agricultural and natural ecosystems, with distinct biodiversity patterns and food web architecture [16,18,19]. For example, intensi ed land use induces shifts from fungi-dominated to bacteria-dominated communities associated with increased N losses and reduced C sequestration [20][21][22], in turn, in uencing consumers at higher trophic levels in the soil food web [16,23]. However, it remains unknown whether changes in complexity of food weds, indicated by interconnectivity among cooccurring species networks [24], in uence the role of biodiversity in ecosystem functioning in agricultural elds. Recent studies have demonstrated that the nature, variability, and intensity of competitive species interactions could shape biodiversity-ecosystem function relationships (BEF) [24,25]. Competition for resources could reduce community functioning if unproductive species become dominant or productive competitors are inhibited [26]. Conversely, diverse communities would be more productive due to their more e cient exploitation of limited resources, in turn, improving community performance [26,27].
However, the ndings above are based on controlled experimental conditions [24,26,28], while most natural communities host highly complex taxonomic diversity [6,29]. Currently, characterizing the in uence of complexity of food webs on BEF relationships in complex real-world ecosystems at large spatial scales is challenging.
Here, our aim was to (I) evaluate the relationships between soil biodiversity and ecosystem functions across complex environmental gradients in agro-ecosystems; and (II) investigate whether food web complexity in uences BEF relationships. We conducted a large-scale soil survey across 228 agricultural elds in China, and obtained diversity information on soil archaea, bacteria, fungi (mycorrhizal and saprotrophic fungi), protists, and invertebrates using high-throughput sequencing of 16S rRNA genes (for soil archaea and bacteria), ITS genes (for soil fungi), and 18S rRNA genes (for soil protists and invertebrates), comprising approximately 60,000 soil phylotypes. We also obtained data on a set of 15 ecosystem functions that are in uenced by soil organisms and correspond to key ecosystem services, such as nutrient provisioning, element cycling, and pathogen control (decreased relative abundance of potential fungal plant pathogens in soils), which represented core ecosystem functions that in uence crop growth and health in agro-ecosystems, and were both fundamental and quanti able. Food-web complexity was visualized based on ecological networks including phylotypes of all types of soil organisms. We applied the moving window approach to explore the consequences of changes in complexity in soil food webs on the strength of BEF relationships. Considering the intensi ed land-use activities in agricultural ecosystems, we explored the following two hypotheses: (I) soil biodiversity is positively correlated with multiple ecosystem functions, while larger soil phylotypes or phylotypes at higher trophic levels have weaker BEF relationships compared to smaller phylotypes or phylotypes at lower trophic levels; (II) an increase in soil food web complexity enhances BEF relationships, i.e., promotes the positive effects of soil biodiversity on ecosystem functions. By considering multiple aspects of soil biodiversity, our study could provide insights into the relationships between soil biodiversity and ecosystem functions in agricultural ecosystems and the potential impacts of shifts in food web complexity.

Results
In agricultural elds across eastern China, the diversity of most single groups of soil organisms exhibited signi cantly positive relationships with averaging multifunctionality (Fig. 1A) and most single ecosystem functions (Fig. 1D), excluding in the cases of protists and invertebrates. This relationship was also observed when considering Shannon (Additional le 1: Fig. 3) and phylogenetic diversity (Additional le 1: Fig. 4). We further explored the relationships between the diversity of major phylotypes in different groups and ecosystem functions. Most of the protist and invertebrate phylotypes showed weaker or no relationships with averaging multifunctionality (Additional le 1: Fig. 5) and single ecosystem functions (Additional le 1: Fig. 6), when compared to those of archaea, bacteria and fungi. Notably, protists, including Chlorophyta and Ochrophyta, which are common phototrophs, exhibited positive biodiversityecosystem function relationships (BEF), while other consumers or predators exhibited negative (e.g. Cercozoa) or non-signi cant (e.g. Ciliophora and Lobosa) BEF relationships. In addition, all the nematode functional groups had negative or non-signi cant BEF relationships. The results suggest that the relationships between soil biodiversity and multiple functions in agricultural ecosystems depend on the type of organisms and the identities of the dominant phylotypes across soil food webs.
The multidiversity of all groups of soil organisms exhibited greater positive associations with averaging multifunctionality (Fig. 1B), with steeper slopes (Fig. 1C and Additional le 1: Fig. 7), and more explained variation when compared to the associations with individual soil phylotypes. Positive association between soil multidiversity and multifunctionality was also observed in both maize and rice elds when examined separately (Additional le 1: Fig. 8). Moreover, such positive relationships were maintained when applying an alternative multifunctionality index that was weighted by three groups of ecosystem services (nutrient provisioning, element cycling, and pathogen control), such that functions from each ecosystem service contributed equally to multifunctionality (Additional le 1: Fig. 7), and when applying multidimensional functionality (Additional le 1: Fig. 9). The consistent results of multiple approaches validate the positive relationships between soil biodiversity and multiple ecosystem functions in agricultural soils.
Subsequently, we performed structural equation modelling (SEM) to test whether the relationship between soil biodiversity and multifunctionality was maintained when accounting for multiple multifunctionality drivers simultaneously (such as climate [MAT] and soil properties [soil pH, total C and percentage of clay]; Additional le 1: Fig. 10, a priori model). The positive effect of soil biodiversity on multifunctionality was maintained after accounting for other key factors ( Fig. 2A; Additional le 1: Table 1), veri ed based on the strong goodness-of-t of our model. The Random Forest models indicated that soil biodiversity was a signi cant and important predictor, although Soil total C appeared to be the most signi cant parameter in uencing multifunctionality (Fig. 2B), explaining 45% of the variation in multifunctionality (Fig. 2D). The spatial patterns of multifunctionality observed were consistent with the soil C spatial patterns ( Fig. 2C and Additional le 1: Fig. 1). However, other environmental factors only had weaker or non-signi cant relationships with multifunctionality ( Fig. 2D and Additional le 1: Fig. 11).
The soil correlation network was established to explore the strong associations (e.g. ecological clusters) among soil phylotypes of different organism groups, and their prediction of multifunctionality. Five major ecological clusters were identi ed, and they included > 91% of the soil phylotypes that strongly cooccurred within the network ( Fig. 3A and Additional le 1: Table 2). The richness of four ecological clusters had positive correlations with multifunctionality ( Fig. 3B and Additional le 1: Fig. 12) and most single ecosystem functions (Fig. 3C), indicating diversity of phylotypes with similar ecological preferences are also key predictors of multifunctionality. Subsequently, we calculated the topological features of the extracted sub-networks by preserving the phylotypes of individual soil samples to estimate the potential complexity of soil food webs. The number of nodes and edges, average degree, clustering coe cient, and graph density, re ecting the network complexity, were all positively correlated with multifunctionality; while average path length, denoting network sparsity, was negatively correlated with multifunctionality (Additional le 1: Fig. 13). The ndings suggest that soil food web complexity can promote the multiple ecosystem functions.
We then applied the moving window approach to explore the determinants of BEF relationships. We selected two window sizes of 20 and 30 samples to validate our ndings, generating 209 and 199 data subsets, respectively. The strength of BEF was re ected by the explained variation of soil biodiversity in relation to multifunctionality. The topological features of the soil networks exhibited signi cant correlations with BEF strength, with positive relationships of features re ecting network complexity, and negative relationships of features denoting network sparsity (Additional le 1: Fig. 14). Consequently, we calculated one index to re ect the potential soil food web complexity using the topological features of the soil networks via multidimensional scaling analysis (Additional le 1: Tables 3 and 4). We observed that soil food web complexity was positively correlated with BEF strength ( Fig. 4A and 4B, Additional le 1: Fig. 15 and Fig. 16). The SEM (Additional le 1: Fig. 17 Fig. 4C and 4D, Additional le 1: Tables 5 and 6). The Random Forest model results also indicated that soil biodiversity was a signi cant and important predictor (Additional le 1: Fig. 18). Soil pH was the most signi cant factor (Additional le 1: Fig. 17), and explained the most variation in BEF strength (Additional le 1: Fig. 15 and Fig. 16). The results based on 20 and 30 sample window sizes exhibited similar trends, validating the positive effect of soil food web complexity on BEF strength.

Discussion
Global change is causing biodiversity loss across many trophic groups [30,31], with potential effects on ecosystem services 6,29, [32] . Soil biodiversity is a key driver maintaining and promoting multiple ecosystem functions [2,[11][12][13][14]. Considering the high functional redundancy of soil microorganisms [33], most studies have focused on soil microbial diversity [12][13][14][15], while ignoring other soil biota, e.g., protists and invertebrates, which are key components of soil food webs in belowground ecosystems [1,34,35]. In addition, evidence for the link between soil biodiversity and multiple ecosystem functions is lacking in human-dominated agricultural ecosystems, and the factors in uencing BEF relationships remain largely unknown in complex real-world ecosystems over large spatial scales.
Here, we provide evidence that soil biodiversity is essential for the maintenance of multiple ecosystem functions in agricultural ecosystems. By considering averaging multidiversity, the diversity of speci c organism groups (archaea, bacteria, fungi, protists, and invertebrates), averaging multifunctionality, multiple individual functions, and the dominant phylotypes strongly co-occurring within soil ecological networks, we revealed the positive effect of soil biodiversity on ecosystem functions. Such relationships were in uenced by the type of organisms and the identities of dominant phylotypes across soil food webs, so that soil phylotypes with larger sizes or at higher trophic levels, e.g., invertebrates or protist predators, appeared to exhibit weaker or no BEF relationships when compared to those with smaller sizes or at lower trophic levels, e.g. archaea, bacteria, fungi, and protist phototrophs. In particular, our study highlights the role of soil food web complexity in promoting not only the multiple functions but also enhancing the positive effects of soil biodiversity on ecosystem functions. Our ndings suggest that management strategies should consider the value of soil biodiversity, particularly soil food web complexity, in the enhancement and maintenance of ecosystem productivity, stability, and sustainability, in agricultural systems under intensi cation.
The role of biodiversity in ecosystem functioning and the maintenance of multiple ecosystem functions could be attributed to a key feature, functional redundancy, which supports the 'insurance' [36] and 'rivetredundancy' [37]hypotheses, in addition to the 'portfolio' effect [38]. For soil microbial communities confronted with a potentially high level of functional redundancy [33,39], here, we demonstrate that greater soil biodiversity could also ensure greater performance in multiple ecosystem functions in agricultural systems. This was largely due to the effects of greater belowground biodiversity enhancing (1) the redundancy effect of the presence of more phylotypes present that support similar functions and (2) the presence of phylotypes that are predominantly associated with different ecosystem functions [15]. Nevertheless, we observed that not all the soil organisms participated equally in the maintenance of multifunctionality. Our results indicate that protist and invertebrate diversity exhibited non-signi cant BEF relationships in general. In addition, all the nematode functional groups, common predators or parasites, exhibited negative or non-signi cant BEF relationships. The major protist phylotypes, including Cercozoa, Ciliophora, and Lobosa, known as consumers or predators [35], exhibited negative or non-signi cant BEF relationships; whereas phototroph protists [35], e.g., Chlorophyta and Ochrophyta, exhibited positive BEF relationships. The results could be explained by the distinct role of the diversity of different soil organisms in supporting multiple ecosystem functions at different levels of functioning [2]. For example, larger soil organisms are essential for the maintenance of high levels of functioning [2], e.g., acting as predators to regulate the ow of resources to organisms (microorganisms) at lower trophic levels; whereas, small organisms are essential for the ne-tuning of multifunctionality (for example, by nutrient recycling) [2].
Considering nutrient cycling is the most important agricultural ecosystem process in uencing crop yield [1], the key functions in the major nutrient cycles could be performed by soil microorganisms, although the soil fauna have been demonstrated to have profound impacts on soil ecosystems and to regulate numerous essential soil processes [40]. Indeed, the intensi cation of human activity could reshape and simplify the architecture of soil food webs, leading to smaller-bodied organisms and fewer functional groups in agricultural ecosystems [16,18,19,40], as well as shifts from fungal-dominated to more bacteria-dominated communities [20][21][22], in turn, in uencing their consumers at higher trophic levels.
The ndings support our hypothesis (I) that soil biodiversity is positively associated with multiple ecosystem functions in agricultural ecosystems, while soil phylotypes with larger sizes or at higher trophic levels (e.g. protist predators and invertebrates) have weaker BEF relationships, when compared to smaller soil phylotypes or phylotypes at lower trophic levels (e.g., archaea, bacteria, fungi, and protist phototrophs).
Soil organisms live within complex soil food webs, forming ecological clusters of strongly co-occurring phylotypes within ecological networks [2,15]. Here, we revealed that strongly co-occurring dominant taxa within the food web could also predict multiple ecosystem functions, with positive associations between soil phylotype richness and multifunctionality within four ecological clusters. The nding indicated that soils that had larger numbers of phylotypes belonging to these four ecological clusters had greater levels of multifunctionality, while those belonging to other ecological clusters might not contribute substantially to multifunctionality. The observation could be supported by the ndings of a recent study [2]. Furthermore, the prediction of the network topological features with regard to the multifunctionality suggest that soil food web complexity could promote multiple ecosystem functions. Complex processes are driven by trophic interactions among soil organisms within food webs [1,15]. For example, large soil invertebrates comminute large amounts of animal and plant litter to generate resources for fungi and bacteria [1]; positive interactions, induced by niche partitioning or facilitation, could promote microbial community functioning [28,41]; microbial interkingdom associations could enhance ecosystem functioning [15] and promote plant health in the model plant Arabidopsis [42]. We also observed multiple potential associations among soil microbial biodiversity that could positively in uence ecosystem multifunctionality. For example, the protist diversity was positively associated with archaeal, bacterial, and fungal diversity, suggesting potential predator-prey associations (Additional le 1: Fig. 3) that could positively in uence multifunctionality [2,35]. The ndings suggest that greater soil biodiversity also facilitates greater association complexity within soil food webs that are jointly required to support multiple ecosystem functions simultaneously.
What was the most intriguing about our ndings was the positive effect of soil food web complexity, re ected by the ecological networks, including soil phylotypes from all types of organisms, based on BEF strength. Soil biodiversity would be a better predictor of multifunctionality when the associations among soil organisms are more complex. Although recent well-manipulated experimental studies have linked competitive species interactions to the direction of diversity-function relationship [24,26], it is still challenging to characterize the factors in uencing such relationships in complex natural ecosystems, largely due to the highly complex taxonomic diversity in natural communities [6,29]. To the best of our knowledge, this is one of the rst studies to link soil food web complexity to the strength of belowground biodiversity-function relationships. Previous studies have demonstrated that changes in soil food web structure could in uence C and N cycling process [16,22]. Species in highly competitive communities often grow less effectively due to intense competition for shared resources [26,28,43]. Weak-intransitive competition within communities made up of strong competitors exhibit negative diversity-function relationships [24]. Conversely, co-existing species, following niche partitioning based on distinct resources or facilitation, could positively interact, which can increase community functional performance [28,41,44]. Such cases support our ndings that the soil food web complexity not only promotes multiple functions but also enhances the link between soil biodiversity and ecosystem functions. Considering intensive agricultural practices are considered to lead to simpler soil food webs [19,40], and changes in the asymmetry of energy channels [22], our results reinforce the view that the elimination of complex belowground species associations can impair the biodiversity-driven ecosystem functioning in agricultural systems under intensi cation.
A few potential limitations should be considered within the context of the present study. First, the potential complexity of food weds was re ected by constructing ecological networks including soil phylotypes from all types of organisms, which was only based on correlation. Correlation network analyses are only a simplistic representation of a complex system. In addition, ecological networks that are based on correlations can yield spurious results, and associations between taxa within such networks cannot be automatically interpreted as interactions. Consequently, it may not be possible to comprehensively depict the architecture and connectedness of soil food webs in real-world conditions.
However, ecological network information is still essential for estimating potential species interrelationships within complex soil food webs, and, in turn, for revealing the in uence of soil food web complexity on biodiversity-driven ecosystem functioning [2]. Secondly, we highlighted the potential limitations of sequencing approaches for the quantifying of soil invertebrate biodiversity; larger soil organisms would be potentially underrepresented based on such an approach, although a few studies have applied the approach to estimate the soil invertebrate biodiversity [2,3,34].

Conclusions
Our ndings reveal the importance of soil biodiversity in supporting and maintaining ecosystem functioning in agricultural systems, by considering average multifunctionality, multiple individual functions, averaging multidiversity, the diversity of distinct groups (archaea, bacteria, fungi, protists and invertebrates), and the dominant phylotypes strongly co-occurring within the soil ecological network.
Furthermore, our results indicate that soil phylotypes with larger sizes or at higher trophic levels (e.g. protist predators and invertebrates) exhibit weaker BEF relationships compared to those with smaller sizes or at lower trophic levels (e.g. archaea, bacteria, fungi and protist phototrophs), suggesting that such relationships vary based on the type of organisms involved, and the identities of the dominant phylotypes within soil food webs. Particularly, we highlight the role of soil food web complexity in promoting multiple ecosystem functions and enhancing our capacity to predict the effect of soil biodiversity on ecosystem function and stability. Overall, our results represent a considerable advancement in the capacity to forecast the effects of belowground biodiversity within food webs on ecosystem functions in agricultural systems. Our work further demonstrates that the potential role of soil food web complexity in the enhancement and maintenance of ecosystem productivity, stability, and sustainability should not be overlooked under land-use intensi cation scenarios. Conversely, it is a key factor that should be take into account in policy development and soil management.

Soil sampling
Soil samples were collected between July and September 2017 from 114 locations in agricultural elds that had been under cultivation for at least 10 years. Maize and rice were selected as representative crops for dryland and wetland soils, respectively, and because they are globally important crops and are cultivated extensively across China. In each location, sampling sites were selected from two adjacent maize and rice elds less than 5 km apart, yielding 114 and 114 soil samples associated with maize and rice cultivation, respectively (see Fig. 1  from south to north. Three 100-m 2 plots were sampled at each site, and ve soil cores obtained per plot at a depth of 0-15 cm were combined. Standard testing methods were adopted for the measurement of soil pH, cation exchange capacity (CEC), total organic C, and percentage of clay, as previously described [45][46][47].

Soil biodiversity analysis
In total, we observed approximately 60,894 taxa in the 228 agricultural elds. The diversity of soil archaea, bacteria, fungi, protists, and invertebrates were analyzed using high-throughput sequencing, by targeting the 16S rRNA region in archaea, the 16S rRNA region in bacteria, the internal transcribed spacer 1 region for fungi, and the 18S rRNA region in protists and invertebrates. Corresponding polymerase chain reaction assays were performed using the Arch519F/Arch915R, 515F/907R, ITS5-1737F/ITS2-2043R, TAReuk454FWD1/TAReukREV3 primer pairs, respectively [48,49], and sequencing was performed on an Illumina HiSeq2500 platform (Illumina Inc., San Diego, CA, USA). The acquired sequences were ltered for quality control as previously described [50,51]. Any chimeric sequences were removed using the USEARCH tool based on the UCHIME algorithm [52]. Subsequently, the sequences were split into operational taxonomic units (OTUs) at a 3% dissimilarity level using the UPARSE pipeline [52]. OTUs with fewer than two sequences were removed, and their representative sequences were assigned to taxonomic lineages using the RDP classi er within the SILVA database (release 128) for bacteria and archaea, UNITE+INSD (UNITE and the International Nucleotide Sequence Databases) for fungi, and PR2 (Protist Ribosomal Reference Database) for protists and invertebrates.
Before we calculated the soil organism diversity, the OTU tables were resampled to a minimum number of sequences from each sample, at 36,880 for archaea, 27,712 for bacteria, 30,369 for fungi, 5393 for protists, and 207 for invertebrates, which showed even sampling depth within each belowground group of organisms (Additional le 1: Fig. 2). Protists were de ned as all eukaryotic taxa, excluding fungi, invertebrates (Metazoa), and vascular plants (Streptophyta), according previous studies [2,35]. Here, we used richness (that is, number of soil phylotypes) as a metric of soil biodiversity, which is the most extensively used, as well as the simplest, metric of biodiversity [2]. We use the term soil biodiversity to refer to different types of richness in general terms. Overall, we calculated the richness of the most prevalent prokaryotic and eukaryotic organisms in our soil samples, including archaea, bacteria, mycorrhizal and saprophytic fungi, protists, and invertebrates. The identities of saprophytic and mycorrhizal fungi were determined using FUNguild [53]. We used only highly probable and probable guilds with an identi ed single trophic mode for the analysis.
To obtain a multidiversity index, we combined soil biodiversity characteristics by averaging the standardized scores of richness of archaea, bacteria, mycorrhizal and saprophytic fungi, protists and invertebrates. The scores standardized based on a common scale ranging from 0 to 1 were calculated according to the following formula: STD= (X−X min )/(X max −X min ); where STD is the standardized variable and X, X min , and X max are the target variable, its minimum value, and its maximum value across all samples, respectively. The multidiversity index generally accepted and has been used extensively in the current biodiversity-function literature [9,14,30,54,55]. Particularly, the richness of each belowground group of organisms was highly correlated with their Shannon diversity (Pearson r=0.75-0.94, P<0.001, Additional le 1: Fig. 3). The multidiversity calculated based on Shannon diversity was highly correlated with that calculated based on richness (Pearson r=0.67, P<0.001), and also showed positive associations with average multifunctionality (Additional le 1: Fig. 3). This trend was also observed when considering phylogenetic diversity (Additional le 1: Fig. 4). In addition, the soil biodiversity indices included in the averaging index were not highly multicollinear (r<0.8 , and resistance to plant pathogens (reduced relative abundance of fungal plant pathogens in soils). Although we have included 15 ecosystem functions, some important functions are inevitably unmeasured, such as process rates including nitri cation rates, denitri cation rates, N mineralization, and decomposition rates, and future studies are encouraged to include more essential functions for comprehensive understanding of ecosystem functioning.
All of these are state variables rather than rates of changes.
In all the soil samples, the extracellular enzyme activities were measured using uorometry as described previously [56]. The soil dissolved organic C, available N, available phosphorus, microbial biomass C, microbial biomass N, available S, available Fe, available Cu, available Mn, and available Zn were all measured using standard soil testing procedures [45][46][47]. The relative abundance of potential fungal plant pathogens in soils was obtained from the sequencing analyses and was inferred by parsing the soil phylotypes using FUNguild [53]. Similarly, we applied only highly probable and probable guilds with an identi ed single trophic mode for the analysis. The inverse abundance (reduced relative abundance) of potential fungal plant pathogens was obtained by calculating the inverse of the variable (total relative abundance of fungal plant pathogens×−1). Because the values for ecosystem functions varied widely, we standardized all variables to a common scale ranging from 0 to 1 (as described above).
To derive a quantitative multifunctionality index for each sample, we used three independent multifunctionality approaches [57]: (1) multiple single functions [2], (2) averaging multifunctionality index [9], and (3) principal coordinate multifunctionality index [58]. To obtain an averaging ecosystem multifunctionality index, we averaged the standardized scores (a common scale ranging from 0 to 1) of all individual ecosystem functions. We also used principal coordinate analysis to examine the different dimensions of multifunctionality [58]. Note that we do not argue that one is better or more appropriate than the other.
Structural equation model and random forest analyses. We used structural equation modelling (SEM) to evaluate the direct link between soil biodiversity and multifunctionality (averaging), and between complexity of soil food webs and strength of BEF (explained below), after accounting for multiple drivers such as climate (mean annual temperature) and soil properties (soil pH, total C, and percentage of clay) simultaneously (Additional le 1: Fig. 9 and 13, priori models). We obtained climatic data including mean annual temperature (MAT) for all sampling sites from the Worldclim database (www.worldclim.org). Here, we did not consider mean annual precipitation since the source of water in the agricultural elds was partly from arti cial irrigation. All the variables were included as independent observable variables. The goodness of t of SEM models was checked using the following procedures: the χ 2 test, the root mean square error of approximation (RMSEA) and Comparative Fit Index (CFI), the model had a good t when the CFI value was close to 1 and the P values of the statistics were high (traditionally>0.05) [59]. With a good model t, we were able to interpret the path coe cients of the model and their associated P values. A path coe cient is analogous to the partial correlation coe cient, and describes the strength and sign of the relationship between two variables. SEM was conducted using the "lavaan" package in R environment (v3.5.1; http://www.r-project.org/) [60]. In addition, random forest (RF) analysis was performed to identify the major driving factors. To estimate the importance of the variables, we used percentage increases in MSE (mean squared error) of variables: higher MSE% values imply more important variables. The analysis was conducted using the "rfPermute" package in R [61]. Signi cance of the model was assessed with 5000 permutations of the response variable using the "A3" package in R [62].
Co-occurrence networks. To explore the species interrelationships within complex soil food webs, ecological co-occurrence networks consisting of dominant soil phylotypes from all types of organisms were constructed. We focused on the most dominant phylotypes-those that were both abundant (top 10% of all identi ed prokaryotes and eukaryotes in terms of relative abundance) and ubiquitous (>50% for prokaryotic and >10% for eukaryotic organisms of all locations) across all soil samples, and identi ed ecological clusters of strongly co-occurring soil phylotypes within the networks. Using such a ltering approach, our aim was to ensure that we identi ed dominant phylotypes in agricultural elds and minimize potential spurious correlations from the rare taxa [2,17,63]. We focused on the dominant soil phylotypes because they are expected to have a disproportionate functional importance in their ecosystems, and are distributed widely. Although many bacterial taxa are widely distributed, it is unlikely to be the case for eukaryotic organisms. Therefore, we applied a ubiquity threshold of >50% for prokaryotic and >10% for eukaryotic organisms in all locations.
Robust correlations with Spearman's correlation coe cients (ρ) >0.65 and false discovery rate (FDR)corrected P-values <0.001 were used to construct networks. We expected that the cut-off, which has been extensively used in literature and is comparable across studies [2], would have both a mathematical and biological meaning and reveal the organisms that are strongly correlated with each other. Nevertheless, we acknowledge that correlation network analyses can yield spurious results, and associations between taxa within such networks cannot be directly interpreted as interactions, and might not fully represent the architecture and connectedness of soil food weds in real-world conditions. However, the information derived from the networks is essential for generating novel hypothesis and ecological frameworks (to be tested in future experiments) on the role of strongly co-occurring phylotypes within food webs in controlling multifunctionality and the link between soil biodiversity and ecosystem functions.
Our network included 1,513 dominant and strongly co-occurring soil phylotypes. The soil phylotypes were dominated by 263 archaea, 1,077 bacteria, 114 fungi, 58 protists, and 1 invertebrate. We identi ed the ecological clusters within our ecological network using the "igraph" package, and calculated the richness of soil organisms within each ecological cluster across all samples. In addition, we extracted subnetworks by preserving the phylotypes of individual soil samples using the induced_subgraph function in "igraph" package in R [64]. The topological features of the sub-networks in each sample were calculated to estimate the potential complexity of soil food webs, including the number of nodes and edges, average degree, clustering coe cient, average path length, network diameter, and graph density. Average path length refers to the average network distance between all pairs of nodes; network diameter refers to the greatest distance between the nodes that exist in the network; average degree refers to the average connections of each node with another unique node in the network; clustering coe cient represents the degree to which the nodes tend to cluster together; and graph density refers to the intensity of connections among nodes [65,66] Therefore, higher numbers of nodes and edges, average degree, clustering coe cient, and graph density, and lower average path lengths and diameters suggest a more connected network, re ecting more potential complexity of soil food webs [45,67]. Networks were visualized using the interactive Gephi platform (https:// gephi.org).
Linking the complexity of soil food weds to the strength of BEF. We applied the moving window approach to explore the determinants of BEF relationships. The technique facilitates the analysis of multivariate data along an ecological gradient [68,69]. To ensure adequate data amounts, we selected two window sizes of 20 and 30 consecutive samples across sites, generating 209 (e.g. 1-20, 2-21 … 209-228) and 199 (e.g. 1-30, 2-31 … 199-228) data points, respectively. The window was advanced across the sampling sites after reordering along the latitude gradient, generating adjacent sampling subsets to better estimate the complexity of soil food weds and BEF relationships. Based on each sub-dataset, we estimated the amount of variance in multifunctionality explained (R 2 of the ordinary least squares linear regressions) by soil biodiversity to re ect the strength of BEF, and constructed soil networks as described above. We then calculated one index to re ect the potential complexity of soil food webs using the topological features of the soil networks via multidimensional scaling analysis [58] (Additional le 1: Table 3 and 4). Note that average path length and diameter, denoting the network sparsity, was calculated as the inverse of the variables (×−1) before the calculation of the index. Subsequently, SEM and RF analysis were performed to evaluate the effect of complexity of soil food webs on BEF strength (as explained above).

Declarations
Ethics approval and consent to participate Not applicable

Consent for publication
Not applicable Availability of data and materials The raw sequence data reported in this paper are available in the NCBI Sequence Read Archive under BioProject PRJNA544819

Competing interests
The authors declare that they have no competing interests.  The relationships between multifunctionality and biodiversity of organisms in agricultural ecosystems.
The linear relationships between multifunctionality and the biodiversity of selected groups of soil organisms (A; number of species, richness) or multidiversity (B; averaged standardized between 0 and 1). Statistical analysis was performed using ordinary least squares linear regressions; P values were indicated by asterisks: *P < 0.05, **P < 0.01 and ***P < 0.  Links between soil biodiversity and multifunctionality in agricultural ecosystems. (A) Structural equation model describing the direct relationship between the multidiversity of soil organisms and averaging ecosystem multifunctionality (EMF). We grouped the edaphic properties into the same box in the model for graphical simplicity, which did not represent latent variables. Numbers adjacent to arrows were indicative of the effect size of the relationship. R2 denotes the proportion of variance explained. Red arrows represented positive paths, and blue arrows represented negative paths. R2 denoted the proportion of variance explained. Signi cance levels were as follows: *P < 0.05, **P < 0.01 and ***P < 0.001. RMSEA: root mean square error of approximation; CFI: Comparative Fit Index; MAT: mean annual temperature; CEC: cation exchange capacity. Information about our a priori model was provided in Additional le 1: Fig. 10 and Additional le 1: averaged over the forest (5000 trees). Percentage increases in the MSE (mean squared error) of variables was used to estimate the importance of these predictors, and higher MSE% values implied more important predictors. (C) The predicted spatial distributions of multifunctionality. The cross-validation ("CV") of the maps was based on Pearson correlation between the predicted and observed values in each sampling site. (D) The linear relationships between multifunctionality and soil C and pH.