Microbial diversity and the abundance of keystone species drive the response of soil multifunctionality to organic substitution and biochar amendment in a tea plantation

High nitrogen (N) fertilizer inputs accelerate soil acidification and degradation in tea plantations, thus posing a threat to soil microbial diversity, species composition, and ecosystem service functions. The effects of organic fertilizer and biochar applications on improving soil fertility have been extensively studied on cropland; however, little is known about their effectiveness in promoting soil multifunctionality on rapidly expanding acidic soils in tea plantations. In this study, we conducted a two‐year field experiment in a subtropical tea plantation to investigate the effects of organic fertilizer substitution and biochar amendment on soil microbial communities and multifunctionality. The results showed that soil multifunctionality was enhanced in plots amended with organic fertilizer and biochar. Soil multifunctionality was significantly and positively correlated with alpha‐diversity of bacteria but not fungi. We also found that organic fertilizer substitution and biochar amendment improved soil multifunctionality by altering the abundance of keystone species. The abundance of keystone species classified as module hubs in the bacterial co‐occurrence network contributed significantly and positively to soil multifunctionality. In contrast, the keystone species categorized as module hubs in the fungal co‐occurrence network negatively affected soil multifunctionality. Soil pH was a key driver of soil microbial community composition, indicating that the increase in soil pH under organic fertilizer and biochar amendment had a crucial role in biological processes. These results suggest that organic substitution and biochar amendment are beneficial in preventing soil degradation and maintaining soil multifunctionality in subtropical tea plantations.


| INTRODUCTION
Complex and diverse soil microbial communities are fundamental to maintaining multiple ecosystem functions (Delgado-Baquerizo et al., 2017;Wagg et al., 2014). Any changes in microbial communities would have crucial impacts on soil multifunctionality by influencing nutrient cycling, litter decomposition, primary production, and climate change (Delgado-Baquerizo et al., 2016;Philippot et al., 2013). Recent investigations have shown that agricultural management practices, particularly fertilization, can severely positively or negatively affect soil multifunctionality (Ding & Wang, 2021;Guo et al., 2021). However, the ecological role of microbial communities in regulating soil multifunctionality under different fertilization practices is not well understood Li, Zhang, et al., 2021;Yang et al., 2021). More empirical evidence concerning how fertilizer management influences the soil microorganisms and multifunctionality and their linkages are essential for the development of sustainable agriculture.
Excessive loading of inorganic N fertilizer in the tea plantation ecosystem has resulted in severe soil acidification, structural damage, and nutrient leaching, leading to negative impacts on soil microorganisms (Ji et al., 2018;Yan et al., 2020). Tea plants (Camellia sinensis L.) are widely distributed throughout tropical and subtropical regions in China, and the area of tea cultivation makes up 70% of the world total (FAOSTAT, 2020). The popularity of tea has recently been increasing steadily because of the abundance of beneficial vitamins, antioxidants, and amino acids in tea (Yang et al., 2021). To obtain rich nutritional ingredients and desirable yield, Chinese primary tea-producing areas receive approximately 500 kg N ha −1 yr −1 (Ni et al., 2019). Several recent studies have shown that bacterial and fungal diversities were drastically reduced in tea-planted soil after the long-term application of chemical fertilizers (Gu et al., 2019;Ji et al., 2020). In turn, all of these changes may negatively impact the availability and turnover of soil nutrients and soil multifunctionality. Therefore, there is an urgent need to prevent soil degradation of tea plantations resulting from the long-term application of inorganic fertilizer using current knowledge-based strategies.
Mounting evidence suggests that organic fertilizer substitution may contribute to improving soil structure, neutralizing soil pH, and adding nutrients simultaneously (Lavallee et al., 2020;Oldfield et al., 2019). Organic fertilizers made from animal manure can also improve agricultural soil health in general by increasing essential micronutrients and biological activity . A vital advantage of the organic fertilizer substitution option is to manipulate soil enzyme activities and microbial community structures (Bowles et al., 2014;García-Ruiz et al., 2008;Kotroczó et al., 2014), which would provide the maintenance or promotion of soil ecosystem functionality (Chaparro et al., 2012). Therefore, manure application may be effective agricultural management to restore soil environments of tea plantations caused by acidification and degradation. A better understanding of how organic fertilizer substitution would affect microorganisms and soil multifunctionality is essential to understand the mechanisms of biotic and abiotic changes caused by organic fertilization in tea-planted soil.
Biochar amendment has attracted significant attention because it can sequestrate organic carbon and enhance soil fertility (Cooper et al., 2020;Wang et al., 2016). The presence of biochar in the soil can also change microbial community diversity and composition in cropland ecosystems (Wang et al., 2015;Yang et al., 2021). However, the impacts of biochar amendment on microorganisms hinge on the initial soil properties, environmental conditions, and the type and application rate of biochar. Biochar addition is beneficial to microbial diversity through increasing soil pH, aeration, and essential nutrient content (Palansooriya et al., 2019). Conversely, the use of biochar greatly reduced the number of bacteria in a Fluvo-aquic soil due to the decreased content of total phospholipid fatty acid (Wang et al., 2015). There was no apparent change in the microbial abundance and biomass in Aridsol soils under biochar amendment, whereas in the same soil when biochar was combined with cattle manure, remarkable differences were found (Elzobair et al., 2016). Thus, it is necessary to probe into the overall effect and complex mechanism of biochar amendment on microbial community diversity and composition in the tea plantation ecosystem toward better use of biochar.
Network analysis has been applied in soil microbial ecology to explore patterns of microbial communities assembly, and it has revealed why some microbial groups consistently occur together or whether specific microbial taxa are more critical for maintaining network structure, and visualize response patterns of taxonomic groups to agronomic practices (Hartman et al., 2018;de Vries et al., 2018). For this reason, investigating the alteration of soil microbial co-occurrence networks driven by organic fertilizer substitution and biochar amendment in tea-planted soil can provide insights into the relationships between soil microbiome and multifunctionality.
In the present study, we aimed to evaluate the response of soil multifunctionality and microbial community structure to organic fertilizer substitution and biochar amendment in a subtropical tea plantation. We also determined whether organic fertilizer substitution and biochar amendment-induced changes in soil multifunctionality are associated with microbial community diversity and composition. We hypothesized that organic fertilizer substitution and biochar amendment would positively affect soil microbial community structures and multiple functionalities. We performed a two-year field experiment with five fertilization treatments to measure a range of soil physical and chemical properties and generated a soil multifunctionality index to test this hypothesis. The 16S and ITS rRNA gene amplicon sequencing were used to quantify the diversity and compositions of soil bacteria and fungi and then the co-occurrences networks were constructed to investigate the associations among microbes. These data were used to identify the relationship between microbial community structures and soil multifunctionality.

| Site description and soil sampling
A two-year field study was established in the Nonglin tea planting farm in Zhenjiang City, China (31.97N, 119.14E, elevation 106 m). The region is characterized by a subtropical humid monsoon climate, with the mean annual air temperature and precipitation being 16.2°C and 1192 mm, respectively . The soil type of experimental site is planosols and contains 35.9 ± 2.1% clay (<0.002 mm), 28.0 ± 4.0% silt (0.002-0.05 mm), and 36.1 ± 2.0% sand (0.05-2 mm). Other soil properties include total N content of 0.32 ± 0.01%, organic C content of 2.14 ± 0.03%, pH of 4.5 (1:2.5 H 2 O), and bulk density of 1.21 g cm −3 .
The two-year field experiment was initiated in July 2018. We set up five treatments with three replicates: conventional chemical fertilizer (N: P 2 O 5 : K 2 O = 16:16:16) (CF, local farm practice), cattle manure (organic fertilizer substitution; OF), CF combined with biochar (CF+Bc), OF combined with biochar (OF+Bc), and the control treatment with no fertilization. Field plots were arranged in a randomized complete block design, and the size of each field plot is 18 m 2 . The distance between the canopies was approximately 0.3 m, and the breadth of the canopy was about 0.8 m (Han, Wang, Xu, Sun, et al., 2021). According to the local practice, N fertilizer was applied at a rate of 450 kg N ha −1 yr −1 , divided into basal fertilizer and top-dressing. The basal fertilizer was applied on October 12 each year at a rate of 150 kg N ha −1 yr −1 , while top-dressing was only used on February 28, 2019, at a rate of 300 N ha −1 yr −1 due to the COVID-19 pandemic lockdown in 2020. Fertilizers and biochar were applied in the inter-row spaces with a width of 0.2 m using the band application method and incorporated into soils at a depth of 0.1 m according to local practice. Biochar was only added together with basal fertilizer on October 12, 2018, at a rate of 20 t ha −1 . Organic fertilizer has a C/N ratio of 7.6, a total N of 2.76%, and a pH of 7.85. Biochar had a total N of 0.63%, an ash content of 14.6%, and a pH of 9.61.
Soil samples were collected at three points from the soil between rows and under the canopy simultaneously within each plot using stainless steel corers on August 21, 2020. Fresh soil for each plot was pooled and passed through a 2-mm sieve to remove visible gravel and plant residues. One portion was stored at 4°C to measure soil functional metrics, and the other was kept at −80°C to analyze soil microbiota.

| Soil analyses
To examine the net N transformation rates, about 30 g of fresh soil samples were placed in 50 ml plastic specimen cups and incubated aerobically for a month in the dark at 25°C (Hart et al., 2018). The rates of net N mineralization (Nm, mg N kg −1 d −1 ) and net nitrification (Nn) were calculated according to changes in the concentrations of mineral N (ammonium (NH 4 + ) and nitrate (NO 3 − )) before and after incubation. We extracted the exchangeable NH 4 + and NO 3 − of fresh soils using 1 M KCl (1:5 w/v). NH 4 + and NO 3 − concentrations in the extracts were measured through a UV-VIS spectrophotometer (U-2900, Hitachi, Japan) using the colorimetric salicylate procedure and vanadate methods, respectively. Soil moisture was determined gravimetrically by drying soil in an oven (24 h, 105°C). Soil pH was measured with a soil-to-water ratio of 1:2.5. We extracted dissolved organic carbon (DOC) using ultrapure water, and a TOC analyzer (TOC-L, Shimadzu, Japan) was used to determine the concentrations of DOC. Soil organic C (SOC) and total N (TN) were measured with a EuroEA 3000 elemental analyzer (EuroVector, Milano, Italy). Soil enzymes are the primary biological mechanism of nutrient cycling and organic matter decomposition, which can be used as metrics of biological diversity, soil fertility, and ecosystem functioning (DeForest, 2009 Saiya-Cork et al. (2002). Activities of hydrolytic enzymes were fluorometrically measured in black 96-well plates using 4-methylumbelliferone and 7-amino-4-methylcoumarin as the substrates. Activities of oxidase enzymes were determined in clear 96-well plates using L-3, 4-dihydroxyphenylalanine as the substrate. Urease activity was measured by the indophenol blue colorimetric method using 10% urea solution as the substrate, and the soil mixture was incubated for 24 h at 37°C (Guan et al., 1986).
Soil multifunctionality was assessed as described previously (Fanin et al., 2017;Lefcheck et al., 2015). In brief, data on soil physicochemical properties and enzyme activities were tested for normality by the Shapiro-Wilk test before analyses. The non-normally distributed data were square root-transformed or their logarithm was considered when necessary. To fairly combine response measured on different scales, we transformed each variable by dividing the absolute value of the maximum observed variable level. We then standardized the variables to the scale of zero to one and derived an average multifunctionality index by taking a mean of the scaled responses for each treatment.

Illumina sequencing
Approximately 0.25 g soil sample was taken for DNA extraction using the DNeasy PowerSoil Kit (Qiagen. Inc.). The integrity of genomic DNA was determined through agarose gel electrophoresis, and the Nanodrop 2000 and Qubit 3.0 spectrophotometer was used for measuring the concentration and purity of genomic DNA. Samples were sequenced via Illumina ® MiSeq (Genesky Biotechnologies Inc.). The V4-V5 hypervariable regions of the 16S rRNA gene were amplified with the primers 907R (5'-CCGTCAATTCMTTTRAGTTT-3') and 515F (5'-GTGCCAGCMGCCGCGG-3'). The ITS1 hypervariable regions of the ITS rRNA gene were amplified with the primers ITS1 (CTTGGTCATTTAGAGGAAGTAA) and ITS2 (GCTGCGTTCTTCATCGATGC). The raw read sequences were processed in QIIME2 (Bolyen et al., 2019). The adaptor and primer sequences were trimmed using the cutadapt plugin. The quality control and identification of amplicon sequence variants were performed using DADA2 plugin (Callahan et al., 2016). Taxonomic assignments of representative sequences were performed with a confidence threshold of 0.8 by a pre-trained Naive Bayes classifier. Illumina raw sequences for bacterial and fungal communities were deposited in the NCBI SRA under bioproject number PRJNA780408.

| Statistical analyses
Statistical analyses were performed in R version 4.0.3 (R Core Team, 2020). Shapiro-Wilk tests and histograms tested the normality of data in the present study, and parameters with unequal variances or non-normal distributions were either square-root or logarithmically transformed. Analysis of variance (ANOVA) was used to determine the differences in soil physical and chemical properties, soil multifunctionality, alpha-diversity metrics, and relative abundance of dominant bacterial and fungal phyla among five treatments. The 'randomForest' package was used for random forest analysis. Principal coordinates analysis (PCoA) analyses were conducted using the 'vegan' package to explore the differences in the bacterial and fungal community structures among the different treatments. Relationships of alpha-diversity metrics and the relative abundance of generalists in the microbial network to soil multifunctionality were established by linear regression analysis using the 'stats' package.
Based on Illumina sequencing data according to Molecular Ecological Network Analyses Pipeline (http://129.15.40.240/mena/), the random matrix theory (RMT) was used for constructing co-occurrence networks. Only the relative abundance of OTUs >0.01% was adopted in the analyses. RMT was achieved to identify 0.81 and 0.80 as the thresholds for bacterial and fungal networks, respectively. All pair-wise Spearman correlations between OTUs were calculated, and the correlations with a Spearman's coefficient of less than 0.65 and a p-value of more than 0.01 were removed (Fan et al., 2019). Visualization of the co-occurrence networks was performed through Gephi version 0.9.2. The abundance of each generalist was calculated by adding the relative abundance of species that belong to it. The topological roles of individual nodes in the network were determined by the threshold values of Zi and Pi (Ling et al., 2016). Nodes were classified into four categories: module hubs (Zi > 2.5), network hubs (Zi > 2.5 and Pi > 0.62), connectors, and peripherals (Zi < 2.5 and Pi < 0.62).

| Soil properties and multifunctionality
Note that we mainly focused on the data from the soil between rows because the fertilizer types and biochar amendment had no influence on multifunctionality or microbial community structure from the soil under the canopy (Table S1; Figures S1 and S2). On the contrary, there was a significant change in soil physical and chemical properties with organic fertilizer substitution and biochar amendment (Table 1). Soil pH under organic fertilizer plots was generally higher than that under synthetic fertilizer plots (p < 0.001). Soil NH 4 + concentrations were more abundant under the CF plots than under the OF plots, but there was no detectable variation on NO 3 − concentrations among the five treatments (p > 0.05). Organic fertilizer application increased the DOC concentrations significantly compared with the control plots (p < 0.05), while there were no differences among fertilization treatments. Regardless of fertilizer type, biochar amendment had greater SOC content and higher urease activities. The activities of CBH, BG, NAG, LAP, and PER under organic fertilizer treatments were generally higher than those under synthetic fertilizer treatments, irrespective of biochar addition. On the contrary, organic fertilization decreased PPO activities by a mean of 104% relative to synthetic fertilizer treatments. Soil multifunctionality was higher under organic fertilizer treatments, while biochar amendment had no noticeable effect on soil multifunctionality regardless of fertilizer type.

| Soil microbial community diversity and composition
The alpha-diversity of soil bacteria was highly variable among the five treatments, while there was no variation in fungal alpha-diversity (Figure 1a,b). Specifically, compared with synthetic fertilizer treatments, the richness and Shannon index of bacterial communities were remarkably higher under organic fertilizer treatments irrespective of biochar addition (p < 0.05). Relative to the control, the richness and Shannon index of bacterial communities were remarkably lower for the plots under synthetic fertilizer application (p < 0.05). Regardless of fertilizer type, biochar amendment generally increased the richness and Shannon index of bacterial communities by an average of 33% and 10%, respectively. Soil multifunctionality was significantly and positively correlated with the richness (r = 0.52, p = 0.045) and Shannon index of bacterial communities (r = 0.54, p = 0.037; Figure 1c). Therefore, the enhancement of soil multifunctionality was associated with increased bacterial diversity. The community structures of bacteria and fungi were significantly affected by organic fertilizer substitution and biochar amendment (Figure 2). The PCoA showed that microbial communities strongly clustered according to treatments, which explained about 49% and 36% of the total variation for bacteria and fungi, respectively ( Figure  2a,b). The bacterial community primarily consisted of members of the phyla Acidbacteria, Proteobacteria, Actinobacteria, Chloroflexi, and Bacteroidetes, accounting for ~80% of the relative abundance of all the phyla (Figure 2c). Relative to the control, the abundance of Acidbacteria was lower for the plots under fertilization treatments (p < 0.05; Table S2). Both organic fertilizer substitution and biochar amendment significantly increased the relative abundance of Proteobacteria T A B L E 1 Effects of full organic substitution for synthetic fertilizer and biochar amendment on the soil properties between rows (p < 0.05). Nevertheless, manure application and biochar amendment did not remarkably change the relative abundances of Actinobacteria, Chloroflexi, and Bacteroidetes (p > 0.05). The fungal communities mainly consisted of Ascomycota, Mortierellomycota, and Basidiomycota at the phyla level ( Figure 2d). Fertilization had only a significant effect on the abundance of Mortierellomycota. When organic fertilizer was combined with biochar, the relative abundance of Mortierellomycota was increased compared with the OF treatment, while there was a decrease in the Mortierellomycota abundance when synthetic fertilizer was combined with biochar relative to the CF treatment. We used the random forest model to determine crucial soil properties that influenced microbial community composition (Figure 3). In the model of bacterial community composition (R 2 = 0.703, p < 0.001), the significant importance of soil properties in decreasing order were soil pH > PPO > AP > SOC > NH 4 + > BG > Nn > LAP > C/N ratio (Figure 3a). Similarly, in the model of fungal community composition (R 2 = 0.345, p = 0.01), fungal community composition was greatly affected by seven soil properties (Figure 3b). The relative importance was as follows: soil pH > AP > DOC > LAP > PPO > Nn > NH 4 + .

| Bacterial and fungal cooccurrence networks
Co-occurrence networks were constructed to assess the different co-occurrence patterns of soil bacterial and fungal communities following organic fertilizer and biochar application. In general, there were 256 nodes and 1015 links in the bacterial network (Figure 4a), while the fungal network contained 113 nodes and 252 links ( Figure  4c). This suggested that bacterial community composition was more complex than fungi in tea-planted soil. Zi-Pi plots were constructed to confirm the topological roles of bacterial and fungal network nodes (Figure 4b,d). The nodes allocated into the network hub, module hub, and connector were the generalists, which might be analogous to key organisms to the microbial community as predicted from network theory (Montoya et al., 2006). Overall, modules hubs and connectors were observed in bacterial and fungal networks, but network hubs were only observed in the bacterial network. Six nodes, OTU1378, OTU45, OTU31, OTU14, OTU132, and OTU782, were categorized as module hubs in the bacterial network, which were highly connected to many nodes in their modules ( Figure  4a,b; Table S3). Eight nodes (OTU166, OTU98, OTU133, OTU583, OTU54, OTU577, OTU20, and OTU99) were classified as connectors that were highly connected to several modules, and two nodes were allocated into network hubs (both module hubs and connectors). Similar to the key OTUs in the bacterial network, there were one node (OTU49) and ten nodes (OTU227, OTU20, OTU196, OTU503, OTU43, OTU256, OTU126, OTU95, OTU188, and OTU689) categorized as module hub and connectors in the fungal network, respectively (Figure 4c,d; Table S4).

| The relative abundance of keystone species
Module hub bacterial OTUs belonged to Opitutus, Gp2, Gp1, Mizugakiibacter, Actinomycelates, and Rhodospirillales ( Figure 5a; Table S3). Relative to the control, synthetic fertilizer treatments increased the relative abundance of module hub bacterial OTUs while organic fertilization decreased it. Regardless of biochar amendment, organic fertilizer substitution increased the relative abundance of Gp1 but decreased the relative abundance of Mizugakiibacter. Network hub bacterial OTUs belonging to GP3 and Sporosarcina, mainly distributed in synthetic fertilization treatments (Figure 5b; Table S3). Biochar amendment combined with synthetic fertilizer increased the number of Sporosarcina but decreased the number of Gp3 compared to unamended treatment. The relative abundance of connector bacterial OTUs belonged to Subdivision3_genera_incertae_sedis, WPS-1_genera_in-certae_sedis, Gp13, Gp3, and Gp1, which increased in the plots under organic fertilization treatments relative to those under synthetic fertilization treatments (Figure 5c; Table S3). However, biochar amendment decreased the relative abundance of connector bacterial OTUs compared to unamended treatments, irrespective of fertilizer type.

F I G U R E 2
Principal coordinate analysis (PCoA) plot of Bray-Curtis distance of bacteria (a) and fungi (b), and the relative abundance of the dominant bacterial (c) and fungal (d) phyla from the soil between rows in a subtropical tea plantation. Control, unfertilized plots; CF, plots with chemical fertilizer; OF, plots with composted cattle manure; CF+Bc, CF plots with wheat straw-derived biochar amendment; OF+Bc, OF plots with wheat straw-derived biochar amendment There was only one module hub fungal OTU for the keystone species in the fungal community, which was classified as Microascaceae and mainly distributed in organic fertilization treatments (Figure 5d; Table S4). There was a nearly two-fold increase in the relative abundance of Microascaceae in the OF+Bc treatment compared to the OF plots. Organic fertilizer substitution decreased the abundance of connector fungal OTUs belonging to Filobasidiales, Xylariales, GS13, Mollisia, Myxocephala, Chytridiomycetes, Agaricomycetes, Sordariales, and Fusarium, regardless of biochar addition (Figure 5e; Table  S4). Pearson correlation showed that soil multifunctionality was mainly influenced by the relative abundance of species categorized as module hubs (Table 2). Overall, soil multifunctionality was negatively and significantly related to the relative abundance of species in the module hub of the bacterial network, while correlated positively and significantly with the relative abundance of species in the module hub of the fungal network.

| Treatment effects on soil multifunctionality
Our results demonstrated that organic fertilizer substitution significantly promoted soil multifunctionality, consistent with previous findings that both biological traits and ecological functions can be improved with manure application (Chen et al., 2020;Li, Zhang, et al., 2021). There may be several reasons for this. First, we found that organic fertilization enhanced soil pH from 4.79 to 5.97, directly or indirectly involving soil quality and ecosystem functioning (Table 1). Among the abiotic factors, soil pH was the only selected environmental predictor that influenced the resistance and resilience of soil multifunctionality . Soil pH also drives the response of soil multifunctionality to environmental disturbances in managed ecosystems, probably through its effect on soil microbial communities (Delgado-Baquerizo et al., 2017). Therefore, the neutralization of soil pH in the plots under manure treatments may be vital for improving soil multifunctionality. Second, the added nutrients under organic fertilizer treatments may be the mechanism responsible for the increased multifunctionality . For example, the DOC and TN concentrations were higher in the soil under the OF treatment compared to the CF treatment. The increase in soil nutrients could positively affect functionality because increased resource availability was beneficial to the metabolism and composition of microorganisms (Tiemann & Billings, 2011). Similar results were also reported that high soil multifunctionality is generally accompanied by high values of soil nutrients (Liang et al., 2017;Qiu et al., 2021). Finally, our study demonstrated that organic fertilizer substitution increased soil enzyme activities (Table 1), which are involved in the soil C-and N-cycling (García-Ruiz et al., 2008). Ros et al. (2003) and Marinari et al. (2000) also described a generalized short-to-medium term increase in C-, N-, P-, and S-cycling-related enzyme activities after organic fertilization practices. Soil enzymes can be regarded as early indicators of soil quality change under contrasting agricultural management practices (García-Ruiz et al., 2008). The previous study has demonstrated that soil multifunctionality was significantly and positively associated with enzyme activity in the rhizosphere soil (Mahmoudi et al., 2021). It is worthy to note that despite the absence of biochar amendment effect on soil multifunctionality, more research is needed to understand how biochar addition would affect soil multifunctionality under various conditions. F I G U R E 3 Main factors affecting bacterial (a) and fungal (b) composition characterized by random forest modeling analysis. The figure shows the random forest analysis results with primary predictor importance (i.e., the percent increase in mean square error [%IncMSE]) of environmental drivers on microbial communities' composition from the soil between rows in a subtropical tea plantation. ** p < 0.01, * p < 0.05 ** ** ** ** ** ** ** ** *

| Linkages between soil multifunctionality and microbial community diversity and composition
Fertilization management significantly shifted the profiles of bacterial alpha-diversity in tea-planted soil, while there was no evident effect on fungal alpha-diversity (Figure 1a,b). This result provided an important implication that the diversity of fungi exhibits a lower sensitivity than bacteria to fertilization in the tea-planted soil. Specifically, we found that organic fertilizer substitution and biochar amendment improved the alpha-diversity of bacteria, which had a positive and significant correlation with soil multifunctionality (Figure 1c). Consistent with our finding, a recent study also found that microbial diversity has vital roles in maintaining ecosystem functioning through a wide variety of ways (Delgado-Baquerizo et al., 2016). For example, microorganisms perform many critical ecological processes such as nutrient cycling and decomposition, thereby can support the essential mechanisms linking belowground and aboveground communities in cropland ecosystems (Cardinale et al., 2011;Van Der Heijden et al., 2008). A great diversity of bacteria can improve the rapid breakdown of litter derived from tea F I G U R E 4 Co-occurrence networks of bacteria (a) and fungi (b) from the soil between rows in a subtropical tea plantation based on RMT analysis from OTU profiles, and Zi-Pi plots of bacterial (c) and fungal (d) based on OTUs topological roles in bacterial and fungal networks. The threshold values of Zi and Pi for categorizing OTUs were 2.5 and 0.62, respectively plantations, enhancing soil organic matter content and thereby fostering the activity of soil microorganisms (Van Der Heijden et al., 1998). Besides, the cooperation of a large and diverse group of microorganisms can promote the process of organic matter degrading from complex and recalcitrant polymers into simpler and more labile monomers (Wardle et al., 2004). Overall, we speculated that bacterial diversity supports soil multifunctionality by altering the nutrient supply and resource distribution, enabling high material processing rates in tea plantation ecosystems. In these respects, our study recommended the use of organic fertilizer substitution and biochar amendment in the tea plantation ecosystem to maintain or improve soil multifunctionality.
Fertilizer management resulted in a dramatic shift in bacterial community composition, consistent with Xiang et al. (2020), who found that bacterial community F I G U R E 5 The relative abundance of keystone species in module hub (a), network hub (b), and connector (c) of the bacterial network, and the relative abundance of keystone species in module hub (d) and connector (e) of the fungal network. Control, unfertilized plots; CF, plots with chemical fertilizer; OF, plots with composted cattle manure; CF+Bc, CF plots with wheat straw-derived biochar amendment; OF+Bc, OF plots with wheat straw-derived biochar amendment composition under long-term inorganic and organic fertilization treatments showed a noticeable difference in the tea-planted soils. At the phylum level, we found the major phylum of bacteria in the tea-planted soil includes Acidobacteria, Proteobacteria, Actinobacteria, Chloroflexi, and Bacteroidetes (Figure 2c). Acidobacteria is oligotrophs that appear in a high abundance in poor-nutrient environments, while copiotrophic such as Proteobacteria favor rich nutrient conditions (Sun et al., 2013;Yang et al., 2021). Consequently, it was not surprising to find that fertilization reduced the abundance of Acidobacteria in our study (Figure 2c; Table S2). Qiu et al. (2021) also found that Acidobacteria are astricted by soil moisture, having a negative correlation with most soil nutrients, but a positive correlation with soil moisture. Regardless of biochar addition, the application of organic fertilizer increased the abundance of Proteobacteria when compared to synthetic fertilizer treatment ( Figure 2c; Table S2). Many Proteobacteria are considered copiotrophic, having a relatively fast growth rate and using various substrates (Spain et al., 2009;Zhang et al., 2016). Increases in the relative abundance of Proteobacteria in the organic fertilizer treatment may be related to organic fertilizer substitutioninduced accumulation of available substrates and nutrients. As mentioned above, the high soil multifunctionality can also explain the increase of Proteobacteria in the organic fertilization plots. Actinobacteria are demonstrated to play a significant role in soil C-cycling due to the high capacity in the decomposition of soil organic matter (Araujo et al., 2020), which is the possible reason for the higher abundance of Actinobacteria under organic fertilizer substitution and biochar amendment treatments. Similarly, the fertilizer type and biochar amendment also significantly changed the composition of the fungal community. The predominant phyla of fungi in the tea-planted soil were Ascomycota, Mortierellomycota, and Basidiomycota (Figure 2d). It has been reported that Ascomycota is the most abundant phyla in tea-planted soil (Ji et al., 2020;, which can tolerate stressful conditions such as low nutrient availability to achieve more efficient resource use in a challenging environment (Chen et al., 2017). Regardless of biochar amendment, organic fertilizer substitution enhanced the relative abundance of Mortierellomycota (Table S2). Mortierellomycota can metabolize organic matter rapidly in soil, and nutrient substances stimulate their abundances under the organic fertilizer treatment . Another reason to explain the higher abundance of Mortierellomycota under the organic fertilizer treatment is that Mortierellomycota has a positive correlation with soil organic carbon content and the recalcitrant C of organic fertilizer contributions a large proportion of the carbon pool in soils of this study (Yao et al., 2017). Besides, the relative abundance of Basidiomycota decreased with organic fertilizer application. Basidiomycota plays a vital role in degrading high-lignin content plant litter . With the organic fertilizer substitution, the proportion of lignin was decreased in the tea-planted soil, probably restricting the growth and reproduction of Basidiomycota.

| Linkages between soil multifunctionality and the abundance of keystone species
Organic fertilizer substitution and biochar amendment altered the abundance of keystone species in microbial co-occurrence networks, which have crucial ecological functions in the microbial community structure. Regardless of biochar amendment, synthetic fertilization increased the relative abundance of keystone species that were classified into module hub in the bacterial co-occurrence network, while organic fertilizer substitution decreased them as compared to the control ( Figure  5a), which meant that synthetic fertilizer appeared to be beneficial for the survival of this ecological cluster. Specifically, synthetic fertilizer application improved the relative abundance of Actinobacteria and Acidobacteria_ Gp2 compared with the control. An interesting finding was that the relative abundance of keystone species in the module hub of the bacterial co-occurrence network negatively correlated with soil multifunctionality (Table  2). Fertilization and biochar addition also affected the relative abundance of keystone species in the network hub and connector of bacterial co-occurrence network (Figure 5b and c; Table S3), although the dominant bacterial taxa showed little control on multifunctionality. This phenomenon may be ascribed to the fact that soil multifunctionality is only associated with specific microbial taxonomic groups (Chen et al., 2019), rare microbial taxa, or specific microbial functional groups in certain regions (Schuldt et al., 2018). The disproportionately role of keystone species in multifunctionality seems counterintuitive given the expected high functional redundancy in microbial communities (Chen et al., 2020). In other words, many species share similar functions, and keystone species might not be necessary for maintaining functions. Overall, our results indicated that keystone taxa in the module hub of bacterial co-occurrence network might have an over-proportional role in soil multifunctionality, which to some extent could explain why fertilization significantly changed the soil multifunctionality.
Our study also demonstrated that organic fertilizer substitution and biochar amendment dramatically changed the relative abundance of keystone species in the fungal co-occurrence network (Figure 5d,e; Table S4). The module hub fungal OTU was classified as Microascaceae, and the relative abundance of Microascaceae was significantly and positively correlated with soil multifunctionality (Table 2). This result implies that the effects of microbial attributes on soil multifunctionality are most limited to those from microbial composition via key microbial taxa. Information on the ecological role of the Microascaceae in natural habitats is scarce, while fungal populations within the Microascaceae were shown to be surprised and positively reactive to degradable labile carbon in soil (Lueders et al., 2006). Therefore, the net activity of soil Microascaceae is essential for ecosystem processes, and the specific contribution of Microascaceae in terrestrial carbon transformations should be explored deeply in future research. In the same way, the fertilization treatment also changed the relative abundance of keystone species in the connector of the fungal co-occurrence network. However, these keystone species have no significant effect on soil multifunctionality. Microbial species considered functionally less relevant under a given environmental condition may be less critical in providing unique functional traits under a favorable condition (Chen et al., 2020). Our results suggested that preserving or even improving the abundance of particular species may be crucial to sustaining the provision of tea plantation ecosystem services in the future.

| CONCLUSIONS
In summary, our experimental findings demonstrated that organic fertilizer substitution and biochar amendment increased microbial alpha-diversity and soil multifunctionality and changed microbial community composition and the abundance of keystone species in the tea plantation ecosystem. The primary driving factor of microbial community structure was soil pH, influenced by fertilizer type. The significant correlations between soil multifunctionality and bacterial diversity and the abundance of keystone microbial species indicated that the fertilization-induced shifts in the structure of microbial communities were the primary drives in shaping the multifunctionality in the tea-planted soil. Altogether, adopting the replacement of synthetic fertilizer by manure or biochar amendment is a potential way to alleviate soil degradation and improve soil multifunctionality in the tea plantation ecosystem.