The pot experiment was carried out in Jiujiang City (29°68′N, 115°98′E), Jiangxi Province, South-Eastern China, which features a subtropical warm and humid monsoon climate, exhibiting four distinct seasons. The paddy soil (classified as Stagnic Anthrosols) was collected from local farmland in Jiujiang before the 2020 rice planting season. There are four treatments with three replicates. Maize straw biochar, which pyrolyzed at 400°C, was used with the dosage of 2.5 g kg^–1 (B1), 5.0 g kg^–1 (B2), and 10.0 g kg^–1 (B3) in biochar treatments, respectively. Pots (20 cm × 20 cm × 25 cm) without biochar amendment were treated as control (CK). In the Spring of 2020, rice seedlings (Taifengyou3301) were planted in plastic pots with 10.0 kg soil. The total number of pots was 4 biochar levels × 2 sampling times (heading stage and maturation stage of rice) × 3 replicates = 24 pots. The pot experiment was initiated in May 2020, with sample collection during the heading stage (n = 36) conducted in August 2020, followed by the maturity stage (n = 36) sampling in October 2020. The experimental duration lasted for 4 months. During the maturity stage sampling, the average seedling height measured 90 cm. After sampling, the soil was separated into three aggregate classes: >2000 μm (large macroaggregate), 250–2000 μm (small macroaggregate), and <250 μm (microaggregate) at each sampling time following the previous study (Elliott, 1986).

In the present study, we quantified 12 metrics related to enzyme activities, nutrient provisioning, and element cycling: soil organic carbon (SOC), dissolved organic C (DOC), total N (TN), total P (TP), ammonium (NH₄⁺), nitrate (NO₃^–), saccharase, β-glucosidase, urease, phosphatase, available S (AS) and available Fe. The contents of SOC, total N, total P, ammonium, and nitrate were measured using standard soil testing procedures (Hu et al., 2021). The concentrations of DOC, available S, and available Fe were detected following the previous approach (Jiao et al., 2022). The extracellular enzyme activities were determined using the model substrates methods as described previously (Yao et al., 2019; Han et al., 2021). Ecosystem multifunctionality measures reflect the ability of an ecosystem to supply multiple functions or services simultaneously. To deliver a quantitative multifunctionality index for each sample, the 12 metrics measured were normalized and standardized using the Z-score transformation. After standardizing, the ecosystem functions were averaged to calculate the multifunctionality index (Delgado-Baquerizo et al., 2016).

Soil microbial DNA was extracted from fresh soil in three aggregate levels by using the Fast DNA SPIN Kit (MP Biomedicals, Santa Ana, CA, USA) according to the manufacturer’s instructions. The diversity and composition of bacteria, fungi, protists, and nematodes were analyzed using sequencing (Illumina HiSeq-PE250 platform). Bacterial 16S rRNA was amplified the V4-V5 region with the primer pairs: 5′-GTGCCAGCMGCCGCGGTAA-3′ and 5′-CCGTCAATTCMTTTRAGTTT-3′. Eukaryota 18S rRNA was amplified using the paired primers: 5′-GGTGGTGCATGGCC GTTCTTAGTT-3′ and 5′-TACAAAGGGCAGGGACGTAAT-3′ (Fan et al., 2021). Bacterial and eukaryotic sequences have been submitted to the NCBI Sequence Read Archive (SRA) with the accession number PRJNA984952.

Raw sequences were processed using the QIIME pipeline and then clustered into operational taxonomic units (OTUs) at a 97% similarity threshold after removing the low-quality sequence. The SILVA database was used to classify bacteria and nematodes, while the UNITE database and the Protist Ribosomal Reference (PR2) database were employed for the classification of fungi and protists, respectively, with all taxonomic assignments adhering to a 97% similarity threshold. The even sampling depths for each group were 12049 for bacteria, 2006 for fungi, 7524 for protists, and 119 for nematodes. Soil biodiversity was characterized as richness (i.e., number of soil phylotypes) and Shannon diversity. The multidiversity index was calculated by averaging the standardized values (0–1) of richness of all types of organisms: soil bacteria, fungi, protists, and nematodes (Wang et al., 2019).

We chose the abundant (top 10% of all identified bacteria and eukaryotes) phylotypes to construct the co-occurrence networks. All pairwise Spearman correlations were calculated in the package “WGCNA” in R and the correlation threshold was above 0.65 and P-value < 0.01. Additionally, we employed the analysis of variance (ANOVA) to eliminate the variance. The data for network analysis is derived from the OTU data of all 72 samples, obtained at various aggregation levels of microbial OTUs. The networks were visualized using Gephi software. The modules were identified from soil microbial groups that were closely interacting with each other. Structural equation models (SEMs) were used to assess the direct and indirect effects of soil pH, biochar, soil biodiversity, and modules on multifunctionality. The data used for the structural equation model (SEM) is based on 72 samples. For bacteria, fungi, protists, and nematode diversity data, we utilized richness data. Regarding module data, we performed principal component analysis (PCA) on the relative abundances of modules 1, 2, 3, and 4 and subsequently extracted the data from the first principal component (PCA1). SEM was conducted in IBM SPSS Amos 21. Linear regressions (based on Spearman correlations) were performed to establish the relationships between soil biodiversity, multidiversity, and multifunctionality. Heatmaps were conducted to estimate the Spearman correlations between biodiversity, modules, and single functions. We also evaluated the role of main protists and nematodes in multifunctionality by linear regressions.

We examined the influence of biochar addition on soil ecosystem multifunctionality at different aggregate sizes. The results suggested that the addition of biochar significantly enhanced the soil multifunctionality compared to the control treatment in the early stage (P < 0.01) (Figure 1A). Moreover, soil multifunctionality increased along with the increasing doses of biochar application in the early stage. For example, the soil ecosystem multifunctionality in B1 and B3 increased by 20.1 and 35.7%, respectively, compared to that in CK (Supplementary Figure 1A). Nevertheless, the impacts of biochar amendment on soil multifunctionality were inconsistent at different stages. The biochar amendment had a significantly positive effect (P < 0.01) on soil multifunctionality in the early stage, but the effect was diminished in the late stag (P > 0.05) (Figures 1A, B and Supplementary Figure 1). These results suggested that the biochar amendment exhibited a stronger positive association with soil multifunctionality in the early stage, and such an effect was enhanced by the increase of biochar application. At the late stage, however, biochar application exhibited weaker or no relationships with soil ecosystem multifunctionality compared to the control treatment. At the aggregate level, no significant correlations were observed for biochar amendment on multifunctionality (Figure 1A). Regarding aggregate size, in the early stage, the highest biodiversity is observed in the microaggregates (Supplementary Figure 2), with a gradual decrease noted in larger aggregate sizes. The large aggregates exhibit the lowest biodiversity during this stage. In the late stage, there is a notable increase in biodiversity for small aggregates, which exhibit the highest biodiversity, followed by microaggregates, with large aggregates displaying the lowest biodiversity among the three categories (Supplementary Figure 2).

The relationship between microbial diversity and soil multifunctionality was explored using the standardized average of twelve variables: SOC, TN, TP, nitrate, ammonium, Fe, DOC, phosphatase, urease, β-glucosidase, saccharase, and AS (see section “2 Materials and methods”). The results indicated that the diversity of soil organisms shifted the soil ecosystem multifunctionality (Figures 2A–D and Supplementary Figure 3). The least-squares regression models revealed a negative linear correlation between soil microbial richness and soil multifunctionality for both bacteria (r = −0.28, P = 0.018) and nematodes (r = −0.35, P = 0.003) (Figures 2A, D). In contrast, the richness of protists was positively related to soil ecosystem multifunctionality along with the soil diversity gradient (r = 0.25, P < 0.05) (Figure 2C). However, the fungal richness showed weaker or no multifunctional relationships than the richness of bacteria, protists, and nematodes (Figure 2B). Interestingly, this relationship was found to be different when considering the Shannon index (Supplementary Figure 3). Unlike in Figure 2B, the fungal Shannon diversity significantly decreased soil ecosystem multifunctionality (P = 0.01), but weaker or no significant associations between the Shannon diversity of bacteria and nematodes and soil multifunctionality were observed (P > 0.05) (Supplementary Figure 3). Notably, the relationship between the diversity of protists and soil multifunctionality was maintained when using the Shannon index (Supplementary Figure 3). The diversity of protists was the only exception presenting a significantly positive effect when both richness and Shannon index were employed as an indicator to evaluate soil microbial diversity (Figure 2C and Supplementary Figure 3).

By identifying the correlation between the diversity of single groups of organisms and single ecosystem functions, it was observed that various soil organisms had distinct impacts on the different aspects of ecosystem functions and services (enzyme activities, nutrient provisioning, and element cycling) (Lefcheck et al., 2015; Jiao et al., 2022). Among them, soil organic carbon (SOC), dissolved organic C (DOC), total N (TN), total P (TP), ammonium (NH₄⁺), and nitrate (NO₃^–) represent nutrient provisioning. Saccharase, β-glucosidase, urease, and phosphatase are indicative of enzyme activities, while available S (AS) and available Fe are associated with element cycling. Activities of P-, N-, and C- cycle enzymes, including phosphatase, urease, β-glucosidase, and saccharase, showed significant correlations with different soil microbes (bacteria, fungi, protists, and nematodes) (Figure 2E). Among the enzymes assessed, the activities of phosphatase, urease, and β-glucosidase were negatively related to the diversity of bacteria and nematodes, while the activities of phosphatase and urease positively correlated with the diversity of protists (Figure 2E). In particular, saccharase was the only enzyme displaying positive association with most of the microorganisms, especially on nematodes (P < 0.001) (Figure 2E). In addition, the accumulation of NH₄⁺ and TP was positively associated with the diversity of bacteria and nematodes but were not positively related to that of fungi and protists. No correlation, however, was found between some single ecosystem functions (including the accumulation of TN, SOC, NO₃^–, and DOC) and the diversity of four soil groups. Fe-element cycling and AS in soil were negatively associated with nematodes but positively with protists (Figure 2E).

To explore whether the biochar addition and multidiversity showed a significant correlation with multifunctionality, we further averaged the standardized scores (z scores) of all ecosystem functions to obtain a single index of soil ecosystem multifunctionality (Supplementary Figure 4). Similarly, the multidiversity index reflects the overall community compositional changes in concert with changes in soil microbial biodiversity. The biochar amendment exhibited highly significant positive effects on multifunctionality (P < 0.001), whereas the changes in ecosystem multifunctionality had negative associations with the average of the indicators of soil biodiversity (P < 0.05) (Supplementary Figure 4). These results suggested that the biochar amendment and changes in soil biodiversity significantly impact ecosystem multifunctionality in paddy soil.

To reveal how soil phylotypes of different modules impact the soil ecosystem multifunctionality, the soil co-occurrence networks were constructed for four soil groups (bacteria, fungi, protists, and nematodes) (Figure 3). The integrated co-occurrence networks were composed of 440 nodes and 2582 edges. As shown in Figure 3A, positive correlations dominated the co-occurrence networks. The proportion of positive associations occupied 99.3% for the total edges of the network. Based on the co-occurrence networks, the multitrophic network was clustered into four modules (subunits with highly inter-connected nodes) (Figure 3A). Co-occurrence network analysis showed that a higher number of nodes and edges among nodes was observed in Module 1, indicating that Module 1 was more closely and complicatedly with higher connectedness than other modules (Figure 3A). We further measured the OTUs proportion of dominant taxa in Modules 1–4 and found that protist OTUs accounted for a considerable proportion of all soil microbial organisms within modules, with the highest relative abundance reaching 48.08% (Figure 3B).

To characterize microbial interactions within each module, we next structured the microbial networks of four ecological assemblies based on the co-occurrence network of total microbial OTUs. Consistent with the results presented in Figure 3A, the network of Module 1 contained a higher number of significantly co-occurring OTUs (nodes = 169) than the other networks. The network complexity within Module 1, as indicated by the number of edges (1557), was also higher than in other modules. The keystone taxa, represented by nodes with both high degree and low betweenness centrality values in the networks, differed among various modules (Figures 4A–D). The keystone species at the phylum level, found in Module 1, were Bacteroidetes and Nitrospirae, and Dikarya, which belong to bacteria and fungi, respectively (Figure 4A). Within Module 2, these were Actinobacteria, Firmicutes, Cryptomycota, and Archaeplastida, which belong to bacteria, fungi, and protists, respectively (Figure 4B). For Module 3, Proteobacteria, Bacteroidetes, and Alveolata which belong to bacteria and protists, respectively, were found to be the keystone phyla (Figure 4C). And the keystone phyla with Module 4 were Bacteroidetes, Acidobacteria, Firmicutes, and Proteobacteria (all belonging to bacteria) (Figure 4D).

It was found that protists displayed a high connectedness in the networks. Likewise, at the genus level, the affiliation of keystone taxa also varied significantly among the modules. Major keystone protistan species within Module 1 included members of the phylum Stramenopiles, comprising the genus Navicula, Sellaphora, and Nitzschia; the phylum Alveolata, consisting of the genus Furgasonia, and the phylum Opisthokonta, comprising the genus Chaetonotida, Monogononta, Gastropoda, and Crustacea (Figure 4A). Within Module 2, the majority of keystone protistan taxa were the phylum Stramenopiles, including the genus Stauroneis and Pinnularia; the phylum Alveolata, including the genus Stylonychia, Cyrtolophosidida, and Didinium, and the phylum Opisthokonta, comprising the genus Monogononta (Figure 4B). Most keystone protistan genera within Module 3 were composed of Neidium and Cyrtolophosidida, while within Module 4, they were Rhabditida and Chromadorea (Figures 4C, D).

In this work, emphasis was also centered on the associations between the biodiversity of each module and single ecosystem functions. The biodiversity within Modules 2–4 was positively correlated with the activity of phosphatase, urease, and β-glucosidase, whereas it was negatively correlated with the activity of saccharase (P < 0.01) (Figure 4E). Besides, the accumulation of TP and NH₄⁺ showed significantly negative correlations (P < 0.01) with the biodiversity of Modules 2–4. The accumulation of TN, NO₃^– and Fe-element cycling was positively related to the biodiversity of Module 2 but not significantly associated with that of most other modules (Figure 4E). No correlation was observed between the accumulation of SOC, DOC, and AS and the diversity of Modules 1–4 (Figure 4E). Interestingly, the significant relationship between the effects of Module 1 on single ecosystem functions appeared to be the exact reverse of Modules 2–4 (Figure 4E).

The relationships between each module and multifunctionality were evaluated to further investigate the mechanism of the modules on adjusting soil multifunctionality. Significant positive correlations (P < 0.05) between the diversity of Modules 2–4 and soil multifunctionality were detected (Figure 5). The richness of Module 1 was the only exception, presenting a significantly negative relationship (P < 0.01) (Figure 5). The specificity of Module 1 (Figure 5) and the characteristics of the microbial co-occurrence networks (Figures 3, 4) indicated that soil microbial taxa in modules regulate multifunctionality. We next explored the importance of dominant taxa within ecological modules as manipulators of soil ecosystem multifunctionality.

Given that the effects of protists on the majority of single functions were consistent with most modules (Figures 2E, 4E), and protists showed specificity in their impacts on soil multifunctionality (Figure 2C and Supplementary Figure 3), we speculated that protists in the modules were the dominant taxa regulating soil multifunctionality. The associations between protistan diversity within each module and soil multifunctionality were then tested. The richness of protists within Modules 2, 3, and 4 was positively related to the multifunctionality (Figure 6). By contrast, the protistan diversity in Module 1 had a negative association with the multifunctionality (Figure 6). The impacts of protists within modules on soil multifunctionality were coherent with those of each module (Figures 5, 6). Besides, we assigned the protists into various functional groups based on their nutrient-uptake mode (Xiong et al., 2020), including phagotrophs, parasites, phototrophs, and plant pathogens, and tested their associations with single functions. Protists, especially microbial-feeding phagotrophs, are significantly correlated with soil multifunctionality (Figure 6). The phagotrophic protists exerted a significant impact on both enzyme activity and the supply of nitrogen and phosphorus nutrients (Figure 6). These results provided independent evidence for the hypothesis of protists regulating multifunctionality.

Structural equation models (SEMs) were then conducted to infer the multiple direct and indirect impacts of biochar amendment, microbial diversity, and modules in explaining soil ecosystem multifunctionality. To test whether various indirect pathways drive the biochar-biodiversity-multifunctionality relationships, we divided the microbial diversity into multidiversity, modules and different soil microbes, including bacteria, fungi, protists, and nematodes, and performed SEMs for both cases separately. Our results from both cases suggested that the direct positive impacts (P < 0.001) of biochar amendment on soil ecosystem multifunctionality were maintained even after accounting for the effects of other ecosystem factors simultaneously (Figure 7). The biochar amendment was the predominant driver of soil multifunctionality and positively related to soil pH which had limited contribution to multifunctionality indirectly (Figure 7). Consistent with the results shown in Figure 2A and Supplementary Table 1, SEMs for soil taxonomic microorganisms revealed that the diversity of nematodes was negatively correlated with soil multifunctionality, whereas that of protists was positively correlated with soil multifunctionality (P < 0.05) (Figure 7A). More importantly, although the bacterial and fungal diversity did not affect soil multifunctionality directly, they did affect the multifunctionality indirectly by altering the diversity of nematodes and protists, respectively (Figure 7A). In addition, modules had direct and significant positive effects (P < 0.001) on soil multifunctionality (Figure 7B). Multidiversity, however, indirectly influenced soil ecosystem multifunctionality primarily by negatively affecting the modules (R = −0.42, P < 0.001) (Figure 7B).

In summary, these results demonstrated that the effects of multidiversity on multifunctionality were indirectly driven by changes in modules. The association of the diversity of soil higher trophic level microorganisms with soil ecosystem multifunctionality was stronger and more significant than that of the lower ones.

Previous studies have confirmed that biochar addition actively enhanced the soil ecosystem functions (Liang et al., 2017; Dong et al., 2022; Han et al., 2022; Ren et al., 2023). In this study, we also found that the different gradients of biochar amendment improved the soil multifunctionality in the early stage (P < 0.01) (Figure 1A). This was largely attributed to the biochar addition increasing soil pH, thus hindering soil acidification (Dong et al., 2022; Figure 7). Soil acidification decelerated nutrient provisioning, litter decomposition, and element cycling, thereby negatively affecting soil ecosystem multifunctionality (Eldridge et al., 2020; Shen et al., 2021; Wei et al., 2022). In addition, with its highly porous and loose characteristics, biochar can promote the growth and proliferation of soil microorganisms by offering more room (Jaafar et al., 2014, 2015; Ye et al., 2017). Furthermore, biochar may enhance soil multifunctionality by improving soil biological properties, such as boosting enzyme activities. Soil enzyme activities are crucial indicators that should be taken into consideration in evaluating the soil ecosystem multifunctionality. Wang et al. (2015) reported that the activities of some extracellular enzymes involved in nutrient cycling were increased in biochar-amended fluvo-aquic soil (Wang et al., 2015). However, the impact of biochar addition on soil multifunctionality was weakened in the late stage, which mainly accounted for the buffering capacity of paddy soils. The soil buffering capacity alleviated the pH enhancement induced by the biochar amendments (Li et al., 2019). Simultaneously, Xu et al. (2012) revealed that the biochar application led to an improvement in the soil buffering capacity (Xu et al., 2012). These results reflected that the biochar amendment on paddy soil exhibited the buffering capacity, which in turn increased with adding biochar. The buffering capacity, nevertheless, was not immediately noticeable in the early stage, while only became evident in the later stage. This supports our first hypothesis that biochar positively affects soil multifunctionality by increasing soil pH, but this effect faded owing to the soil buffering capacity at the late stage. Moreover, aggregate size was closely related to biodiversity levels in biochar-amended paddy soils, with smaller aggregates supporting a higher diversity of soil organisms (Supplementary Figure 2). This indicated that smaller soil aggregates exhibited higher biodiversity compared to larger aggregates. The microscale spatial heterogeneity and complexity of soil aggregates impacted the distribution of microorganisms, concurrently influencing microbial diversity in biochar-amended paddy soils (Han et al., 2021).

The biodiversity of soil protists was positively related to soil multifunctionality, whereas that of bacteria, fungi, and nematodes exhibited negative correlations with multifunctionality (Figures 2A–D and Supplementary Figure 3). SEMs further revealed that the direct associations between the diversity of protists and nematodes and soil multifunctionality were robust after accounting for four soil microorganisms simultaneously. By contrast, bacterial and fungal diversity was indirectly correlated with multifunctionality (Figure 7B). This is corresponding to our second hypothesis that the relationship between the biodiversity of soil higher trophic level microorganisms and soil ecosystem multifunctionality was stronger and more significant than that of the lower ones. Moreover, multiple potential associations between soil microbial biodiversity impacted ecosystem multifunctionality (Delgado-Baquerizo et al., 2020). Specifically, the diversity of protists and nematodes showed a positive association with that of fungi and bacteria, respectively (P < 0.01) (Figure 7B), indicating the presence of potential predator-prey relationships. Bacteria and fungi may indirectly impact multifunctionality by being preyed upon by protists and nematodes (Han et al., 2022). For instance, the negative relationships (P < 0.01) between the typical nematodes (Aracolaimida and Chromadorida) and multifunctionality were observed (Supplementary Table 1). According to the SEM, there was a positive correlation between the diversity of nematodes and that of bacteria (Figure 7B). Considering these associations and interactions among different trophic levels, it could indicate that nematodes preying on bacteria involved in soil functions lead to a significant negative correlation between nematodes and multifunctionality. This is aligned with the earlier finding of Wagg et al. (2019), revealing that the interactions among the multitrophic members of soil microbiome also have contributed to soil ecosystem functions directly or indirectly (Wagg et al., 2019).

Overall, the alterations in soil multifunctionality exhibited a negative association with the average of soil biodiversity (multidiversity) (Supplementary Figure 4), indicating that changes in soil microbial diversity negatively impact soil ecosystem multifunctionality. Currently, the effects of soil microbial diversity on ecosystem processes remain largely controversial because the positive, negative, and neutral relationships between soil multidiversity and multifunctionality have been reported (Nielsen et al., 2011; Bradford et al., 2014; Wagg et al., 2014; Delgado-Baquerizo et al., 2016, 2020; Chen et al., 2020; Dong et al., 2022; Jiao et al., 2022; Fan et al., 2023). These contradictory findings reported in prior studies may be attributed to two main reasons: (1) environmental change factors, including biotic and abiotic factors, can alter the effects of soil diversity on soil ecosystem functions (Hu et al., 2021; Yang et al., 2022); and (2) the soil biodiversity-multifunctionality relationships may be contingent upon the number of functions and their combinations (Bradford et al., 2014; Lefcheck et al., 2015; Soliveres et al., 2016; Chen et al., 2020). All soil ecosystem functions have the potential to be either positively or negatively influenced by soil diversity or the interplay among the functions (Kareiva et al., 2007; Lefcheck et al., 2015). As more functions and trophic levels are considered, the biodiversity-multifunctionality relationships become more comprehensive and persuasive (Soliveres et al., 2016; Wang et al., 2019). Thus, the soil biodiversity-multifunctionality relationships shift as the spatial-temporal environmental conditions change. Different ecosystems have distinct purposes and provide various services; hence, the results of measuring the relationship above vary in different ecosystems (Jiao et al., 2022; Fan et al., 2023).

The correlations of module diversity with soil single functions were stronger and more significant than those of soil microbial taxa (Figures 2E, 4E). Enzyme activities and the accumulation of NH₄⁺ and TP were significantly associated with the diversity of modules (Figure 4E). The stronger relationships between modules and single functions are probably explained by the fact that the ecological clusters can support similar functions (Jiao et al., 2022). Theoretically, species with similar ecological preferences are more likely to cluster and form modules that promote similar functions (Delgado-Baquerizo et al., 2020; Wang et al., 2022). Interactions among soil horizontal and vertical species lead to altering ecosystem functions (Duffy et al., 2007). The microbial co-occurrence networks were dominated by positive correlations (Figure 3A), suggesting potential cooperative interactions among taxa within modules. The cooperative interactions among the species in modules result in their synergistic and complementary relationships that consistently enhance the impact of modules on various ecosystem functions (Chen et al., 2022). No correlation, however, was found between some single ecosystem functions (including the accumulation of TN, SOC, AS, and DOC) and the diversity of modules (Figure 4E). This discrepancy is largely triggered by functional redundancy as soil microbial communities may have overlapping functions (Allison and Martiny, 2008), implying that multiple taxonomic groups can perform the same ecosystem functions (Li et al., 2021). However, the existence of functional redundancy in ecosystems affects evaluating the relationships between biodiversity and soil ecosystem functions. Microbial communities display highly functional redundancy for some basic functions, such as microbial respiration, whereas redundancy may not be observed for specific or specialized functions (Langenheder et al., 2006; Reich et al., 2012; Chen et al., 2020). The functional redundancy tends to fade as more functions are measured (Hector and Bagchi, 2007; Isbell et al., 2011). It also explains why multitrophic levels and more appropriate functions are needed to be considered in exploring the mechanism of biodiversity-multifunctionality relationships.

Four primary modules were detected based on the soil co-occurrence network of soil organisms to reveal how soil modules regulate the soil ecosystem multifunctionality. Soil microbial community diversity and multifunctionality exhibited a positive correlation in Module 2–4 (P < 0.05), yet a negative correlation in Module 1 (P < 0.01) (Figure 5). The specificity of Module 1 has been primarily attributed to closer and more complicated connectedness and more robust interactions between species of Module 1 than other modules (Figure 3A). This result revealed that trade-offs exist when measuring soil biodiversity-multifunctionality relationships (Byrnes et al., 2014; Hu et al., 2021). The negative relationships and complicated interactions within Module 1 reflected an increasing frequency of trade-offs among ecosystem functions as more species and their interactions were considered (Lefcheck et al., 2015).

According to Fan et al. (2021), keystone taxa within modules played a pivotal role in regulating the ecosystem functions of terrestrial ecosystems subjected to 40 years of fertilization (Fan et al., 2021). Keystone taxa, integral to the flow of energy and materials, substantially contribute to the productivity of the soil ecosystem in ecological networks (Toju et al., 2018; Shi et al., 2020). Protists in Module 1 exhibited a negative correlation with multifunctionality, suggesting a unique interaction dynamic within this module that may involve complex trade-offs. Conversely, protists in Modules 2, 3, and 4 showed positive correlations with multifunctionality, highlighting their fundamental contribution to promoting a range of soil functions (Figure 6). Therefore, we identified the protistan keystone taxa in the modules based on the microbial co-occurrence network analysis (Figures 4A–D). Moreover, the major phyla of predatory protists, such as Amoebozoa, Stramenopiles, Excavata, and Alveolata, were found to positively correlate with ecosystem multifunctionality (Supplementary Table 1). The predation or grazing of protists was essential for directing carbon and energy flow toward higher trophic levels while facilitating the release of dissolved nutrients throughout the food webs (Bjorbaekmo et al., 2020). Guo et al. (2022) demonstrated that predatory protists interacted with bacteria in the plant-soil microbiome, ultimately augmenting plant health and crop yields (Guo et al., 2022). In our study, the biochar amendment of paddy soil enhanced the abundance of predatory protists (Asiloglu et al., 2021) and played a critical role in maintaining ecosystem functioning.

The SEMs demonstrated that the modules exerted a significant effect on soil multifunctionality (P < 0.001), while multidiversity indirectly affected multifunctionality through its impact on modules (P < 0.001) (Figure 7A). Simultaneously, the influence of modules on soil multifunctionality surpassed that of individual microorganisms (Figure 7B). Recent research also revealed that modules regulated soil multifunctionality more remarkably than individual microorganisms and multidiversity (Delgado-Baquerizo et al., 2020; Jiao et al., 2022; Yang et al., 2022). These findings convincingly illustrated the crucial role of modules formed by soil microbes with similar ecological preferences as critical drivers of ecosystem functioning. Together, our study highlights the more prominent role of modules in determining soil multifunctionality.