A long-term field experiment, i.e., a rice-wheat cropping system, was established at Ezhou (114.7°E, 30.4°N), Hubei Province, China, in 1981. Wheat is sown in late November and harvested in late May, and rice is sown in early June and harvested in early October. According to the National Data Center for Meteorological Sciences,² from 2019 to 2020, the average precipitation was 3.2 and 2.7 mm, and the average temperature was 27.3 and 12.5^°C in rice and wheat season, respectively, per day. The soil is derived from yellow-brown paddy soil, classified as Cambisols (Nachtergaele, 2017). The soil before cultivation in 1981 had the following properties: pH 6.3, SOC 15.9 g kg^–1, TN 1.8 g kg^–1, available phosphate (AP) 5.0 mg kg^–1, available potassium (AK), 98.5 mg kg^–1. The experiment had a complete randomized design with different fertilization regimes and three replicates. Each plot measured 5 m × 8 m. The fertilization regimes in this study contained: (1) no fertilizer, NoF, (2) inorganic nitrogen, phosphate and potassium fertilizers, NPK, (3) manure, M, and (4) the application of NPK combined with manure, NPKM. We applied urea, superphosphate, and potassium chloride as N, P, and K fertilizers, respectively. The detailed application rates of fertilizers are displayed in Table 1. The organic manure was mainly composed of swine manure (69% moisture) and had the following properties: 262.2 g kg^–1 SOC, 15.1 g kg^–1 TN, 20.8 g kg^–1 phosphorous pentoxide, and 13.6 g kg^–1 potassium dioxide.

Five intact soil cores were collected in October 2019 and May 2020, using stainless steel cores (5 cm inner diameter × 10 cm height) from a depth of 0–10 and 10–20 cm in each plot, which finally contained 48 samples (4 fertilization regimes × 2 crop season × 2 soil depth × 3 replicates). The fresh soil was transported to the laboratory immediately, sieved through a 2 mm mesh, and divided into three subsamples. One subsample was prepared for microbial DNA extraction, one was stored at 4^°C for the analysis of extracellular enzyme activity, and the third was air-dried to analyze the soil physicochemical properties. SOC and TN were determined by an Elemental Analyzer (VARIO MACRO cube, Germany).

The extraction of soil amino sugars, including GlcN, GalN, and MurN, was determined following Zhang and Amelung (1996). Briefly, a weighted air-dried soil sample (contained about 0.3 mg N) was hydrolyzed with 6 mol l^–1 HCl at 105^°C for 8 h, and then cooled to room temperature. The solution was mixed with 100 μl inositol and dried by a rotatory evaporator. The precipitate was dissolved in distilled water, adjusted to pH 6–6.8, centrifuged at 3,000 r min^–1, and freeze-dried. Then, methanol was added to freeze-dried supernatants, centrifuged at 3,000 r min^–1, transferred to vials, and dried by N₂ at 45^°C. Next, the residue was dissolved in 1 ml distilled water and 100 μl N-methylglucamine, mixed, and freeze-dried. The purified amino sugars were reacted with ethyl acetate and n-hexane. The amino sugar derivatives were separated on a gas chromatograph equipped with an HP-5 column (30 m × 0.25 mm × 0.25 mm) and quantified using a flame ionization detector (GC-7890B, Agilent, United States). The temperature program was set as follows: 120^°C for 1 min, increased by 10^°C min^–1 to 250^°C and held for 2.5 min, increased by 20^°C min^–1 to 270^°C and held for 2 min. The injector temperature was 250^°C and the split ratio was 10: 1. The peaks were identified by comparing sample retention times with those of pure amino sugar standards using the Chemstation software. Myo-inositol was added as the first standard after hydrolysis and N-methylglucamine was added as the second standard before the derivatization to assess the recovery of amino sugars.

The contents of amino sugars, microbial residue-C and N were calculated using the analytical equations developed by Appuhn and Joergensen (2006) and Joergensen (2018).

where m_i represents the quality of inositol, and A_i and A_x present the peak area of inositol and amino sugars, respectively. In addition, 45 is the conversion coefficient of MurN to BRC, 6.67 is the conversion coefficient of MurN to bacterial residues-N (BRN), 179.17 is the molecular weight of GlcN, 9 is the conversion coefficient of fungal GlcN to FRC, and 1.4 is the conversion coefficient of GlcN to fungal residues-N (FRN).

DNA was extracted from 500 mg soil using the FASTDNA Spin Kit (MP Biomedical, Santa Ana, CA, United States) according to the manufacturer’s protocols. The DNA concentration and quality were checked by NanoDrop (Thermo NanoDrop, 2000, United States). For bacteria, the V3-V4 hypervariable regions of the 16S rRNA gene were targeted using primers 515F (GTGCCAGCMGCCGCGGTAA) and 806R (GGACTACVSGGGTATCTAAT); for fungi, the ITS1 region was amplified using ITS1F (CTTGGTCATTTAGAGGAAGTAA) and ITS2R (GCTGCGTTCTTCATCGATGC). Sequencing was performed on the Illumina MiSeq platform. Sequenced paired-end reads were joined using FLASH³ and were filtered according to quality using Fastp tools.⁴ Briefly, sequences were joined (overlapping pairs) and grouped by samples following the barcodes and then the barcodes were removed. Then, sequences <50 bp or with ambiguous base calls were removed. Reads were denoised by unoise3, and any chimeric sequences were removed using the USEARCH tool based on silva_16s_v123. Sequences were then split into operational taxonomic units (OTUs) at 97% similarity using the UPARSE pipeline in USEARCH. Plastid and non-bacteria were removed by USEARCH based on silva_16s_v123. This involved quality filtering and 97% clustering of the ITS1 region as indicated above for the 16S processing, using the UNITE database (utax_reference_dataset_all_02.02.2019), for chimera removal, taxonomic identification, and non-fungal removal of representative OTUs. The bacterial and fungal OTU abundance tables were rarefied to an even depth of 30,000 reads for 16S and 25,000 for ITS sequences. The raw sequencing data were deposited in the NCBI Sequence Read Archive (SRA)⁵ with accession No. PRJNA796782.

The data were processed using Microsoft Excel 2019 and the analysis of variance (ANOVA) was conducted by using SPSS 18.0 (IBM, United States), followed by a Fisher’s least significant difference (LSD) test at p < 0.05 to determine the effects of fertilization regimes, crop season, and soil depth on the contents of C and N in soil and microbial residues. PERMANOVA was conducted using the adonis () function with 9,999 permutations, variance portioning analysis (VPA) was determined using the varpart () function, and differences within-group centroids were calculated using the betadisper () function. All these functions were from the ‘‘vegan’’ package⁶ in R (version 3.6.1). Principal co-ordinates analysis (PCoA) was conducted to calculate the differences within microbial community compositions using the BetaDiv () function in “phyloseq” package, which was based on the Bray-Curits distance (Adonis test) with 999 permutations.

Partial least squares path model (PLS-PM) was performed using the “plspm” package (Sanchez et al., 2015) to further investigate the potential direct and indirect effects of inorganic and organic fertilization, crop season, and soil depth on soil C and N contents by influencing the compositions of microbial communities and their residues, which finally influence soil C pools. Bacterial community was represented by the sample values of PC1 and PC2 derived from PCoA analysis. Bacterial residues contained BRC and BRN, while fungal residues contained FRC and FRN. The quality of the PLS-PM was evaluated by examining the goodness of fit (GoF) index, in which >0.7 means an acceptable overall prediction performance of the model, and by calculating the coefficients of determination of the variables (R²). Then, we estimated the standardized total effects of all variables upon soil SOC and TN by using the “pls” package.

To explore the complex relationships between microbial communities and their residues, we further constructed the networks using the Python language package Networkx (see text footnote 1). Only OTUs with top 500 abundances were considered. Edges were represented as links between microbial OTUs and the soil variables, which were tested by the Pearson correlation analysis. Only edges with p < 0.05 and |ρ| > 0.5 were included in the final networks. The network properties, including node degree, degree centrality, closeness centrality, and betweenness centrality, were also calculated. A node with the highest score of closeness centrality, betweenness centrality, and degree centrality are indicated as keystone taxa (Chen et al., 2020). And the networks were visualized using the Matplotlib Python package.⁷

Some figures were plotted using OriginPro 9 (OriginLab, United States) and modified using Adobe Illustrator CS6 (Adobe, United States).

Fertilization regimes and crop season exerted significant impacts on the contents of SOC and TN (p < 0.05), while soil depth did not exert any impact (Table 2, p > 0.05). For fertilization regimes, compared with the original soil, the long-term fertilization enlarged SOC and TN pools, with an increase of 3.8–31.4% and 1.5–19.2%, respectively, while the long-term non-fertilization significantly decreased, with an average decrease of 15.8 and 6.6%, respectively. Compared with soils in NPK, the long-term application of organic manure, especially NPKM, significantly increased SOC and TN contents by 42.2 and 24.6%, respectively. Compared with soils in rice season, a relatively higher TN content was observed in wheat fields, which increased by 12.7%.

For the contents of soil microbial residues, fertilization regimes and crop season provided significant contributions to their variations, especially crop rotation. Meanwhile, soil depth only exerted significant impacts on the contents of FRC and FRN. In detail, when compared with original soil, the long-term fertilization significantly increased the contents of BRC and BRN by 29.0–220.5% and 28.9–221.0%, respectively, while NPKM significantly increased the FRC and FRN by an average of 13.4 and 13.5%, respectively. Compared with NoF, NPKM significantly increased the accumulations of BRC, BRN, FRC, and FRN by 12.7–52.6%, 12.8–53.3%, 13.7–59.2%, and 13.6–59.3%, respectively. For crop season, the accumulations of microbial residues were relatively higher in the wheat fields, and the contents of BRC, BRN, FRC, and FRN were on average 24.2, 24.4, 33.7, and 33.8% higher than those in the rice fields. Compared with the upper soil layer, the contents of FRC and FRN were significantly higher in the lower layer with an increase of 33.7 and 33.8%, respectively.

Moreover, the ratios of microbial residues to SOC and TN can reflect the contributions of the bacterial or fungal residues to soil C and N pools. The ratios of FRC/SOC and FRN/TN accounted for 23.7–51.1% and 30.0–55.0%, respectively, which were significantly higher than the ratios of BRC/SOC (9.1–15.7%) and BRN/TN (10.9–19.9%).

Differences in the overall microbial community compositions based on the PCoA revealed that both bacterial and fungal community compositions demonstrated a significant structuring in response to fertilization regimes, crop season, and soil depth (Figure 1A and Table 3). In detail, soil depth made the highest contribution to the variation of bacterial community compositions, and organic manure contributed more to the variation of fungi (Figure 1B). Major bacterial and fungal phyla were observed in all treatments (Figure 1C). Within fertilization regimes, compared with NoF, NPKM significantly increased the relative abundances of Planctomycetota, Verrucomicrobia, and Ascomycota, with an increase of 22.6, 15.1, and 13.3%, respectively. On the contrary, NPKM significantly increased the relative abundances of Chloroflexi, Gemmatimonadetes, Basidiomycota, Glomeromycota, and Chytridiomycota, with a decrease of 24.2, 13.1, 108.9, and 10.7%, respectively. For crop season, only the relative abundance of Bacteroidetes was 20.7% higher in the wheat fields than that in the rice fields. Compared with the lower soil layer, the relative abundances of Actinobacteria, Firmicutes, Alphaproteobacteria, Gammaproteobacteria, Rozellomycota, Mortierellomycota, and Chytridiomycota were significantly higher in the upper layer, while the relative abundances of Chloroflexi, Bacteroidetes, Nitrospirae, and Glomeromycota were significantly lower in the upper layer.

Regression of Bray-Curtis dissimilarity with soil variables revealed the influence of soil C and N on microbial community compositions (Table 4). For soil microbial communities in rice fields, the contents of soil TN, BRC, BRN, FRC, and FRN significantly correlated with bacterial and fungal community compositions, and SOC only affected soil fungal community compositions (R² = 0.66, p < 0.001). In the wheat fields, FRC and FRN (bacteria R² = 0.23, p < 0.05; fungi R² = 0.26, p < 0.01) were significantly correlated with microbial community compositions, and fungi was also affected by SOC (R² = 0.53, p < 0.001) and TN (R² = 0.56, p < 0.001).

Considering the significant correlations between microbial community compositions and microbial residues, we constructed a PLS-PM to access the direct and indirect effects of fertilization regimes, crop season, and soil depth on the contents of SOC and TN by affecting living and dead microorganisms in the soil (Figure 2). We found that inorganic fertilization increased soil SOC contents directly and indirectly by positively correlating with soil bacterial community compositions and increasing the contents of microbial residues C and N. Both organic fertilization and microbial communities exerted indirect positive impacts on SOC and TN contents. Relative to bacterial residues, fungal residues made more contributions to the contents of SOC and TN.

To further investigate the correlations between soil-specific microorganisms and microbial residues, we conducted a network analysis (Figure 3). Results revealed that both C and N in soil and microbial residues showed significant links to specific microbial OTUs. For the bacterial community in the rice fields, most OTUs affiliated to the phyla Actinobacteria, Acidobacteria, Chloroflexi, and Gammaproteobacteria showed significantly positive correlations with microbial residues, while the OTUs affiliated to the phyla Nitrospirae and Myxococcota showed significantly negative correlations with the microbial residues. On the contrary, in the wheat fields, there were fewer correlations between bacterial OTUs and bacterial residues. OTUs affiliated to Acidobacteria and Methylomirabilota were positively correlated with the fungal residues, while OTUs affiliated with Actinobacteria, Gammaproteobacteria, Alphaproteobacteria, and Gemmatimonadetes were negatively correlated. Considering their relative abundances, only Actinobacteria, which correlated positively with fungal residues, showed relatively higher abundance (1.85–6.42%, Supplementary Figure 1). Among all the identified bacterial genera, Gaiella (0.28–0.93%) and Marmoricola (0.17–0.86%) showed relatively higher abundances. Specifically, fertilization regimes and soil depth affected the genus Gaiella, which showed relatively higher abundances in the upper soil layer with M and NPKM treatments (Supplementary Figure 2). For the fungal community, most OTUs affiliated with the class Eurotiomycetes presented negative correlations with bacterial residues, and OTUs from the class Sordariomycetes showed negative links to the fungal residues, including the genera Cosmospora, Nigrospora, and Zopfiella. All these genera showed the lowest relative abundances in NPKM treatment than other fertilization treatments (Supplementary Figure 3). However, based on their relative abundances, only the class Pezizomycetes, which correlated positively with the fungal residues, presented a relatively higher abundance (0.22–61.09%), which was primarily attributed to the genus Pseudaleuria.

Compared with fertilization regimes, crop season, and soil depth, only fertilization regimes significantly affected the relative abundance of Psedaleuria (Figures 4B, P < 0.001), which showed higher abundances in M and NPKM with an average abundance of 34.1%. Considering its high relative abundance and significant correlations with fungal residues, we constructed the PLS-PM (Figure 4A) to further investigate the possible effects on microbial residues caused by Psedaleuria. Results revealed that the relative abundance of Psedaleuria was significantly positively correlated with fungal residues, which then indirectly increased the stocks of SOC and TN (Figures 4C,D, p < 0.01).