The long-term fertilization experiment commenced at the Yingtan National Agroecosystem Field Experiment Station of the Chinese Academy of Sciences (28°15′20′′N, 116°55′30′′E) in Jiangxi Province, China. The experiment site has a typical subtropical climate with a mean annual temperature of 17.6°C and precipitation of 1,795 mm. The soil is derived from Quaternary red clay, which is classified as Ferric Acrisol according to the FAO classification system and Udic Ferralsol according to Chinese Soil Taxonomy. The long-term field experiment followed a completely randomized design with three replicates. The experiment was conducted since 2002, which consisted of 12 concrete plots with the following size: 2-m long, 2-m wide and 1.5-m deep. The four manure treatments were (1) no manure (M0); (2) low manure with 150 kg N ha⁻¹ y⁻¹ (M1); (3) high manure with 600 kg N ha⁻¹ y⁻¹ (M2); and (4) high manure with 600 kg N ha⁻¹ y⁻¹ and lime applied at 3,000 kg Ca (OH)₂ ha⁻¹ 3 y⁻¹ (M3). Notably, the experiment was not a full factorial design as lime was only applied to the high manure treatment. Such an experimental design allowed us to determine the substantial effects of manure treatments on soil nematodes and microbial communities in the macroaggregates, while the impact of lime was just estimated under the high manure application. Pig manure had an average total carbon content of 386.5 g kg⁻¹ and a total nitrogen content of 32.2 g kg⁻¹ on a dry matter basis. The field was planted annually with corn monoculture (cultivar No. 11 from Denghai) from April to July. No management measures were taken with the exception of weeding by hand.

The topsoil samples (0–20 cm) were collected after the harvest of maize in late July 2013. Ten soil cores were collected from each plot using an auger and pooled together to form a composite sample. Fresh samples were chilled on ice immediately after collection in the field and transported in a cooler to the laboratory. Soils were forced through a 4-mm sieve to remove visible residue and then homogenized. Soil was manually fractionated into three aggregate sizes, including large macroaggregates (>2 mm), small macroaggregates (0.25–2 mm), and microaggregates (<0.25 mm) (Jiang et al., 2013). In this study, only large macroaggregates were used for the subsequent stability analysis and micro-computed tomography scanning and thereafter referred to as macroaggregates. The macroaggregates were further subdivided into three subsamples for analyzing soil physicochemical properties, nematode assemblages, and microbial communities.

Soil pH was determined using a glass electrode in a soil: water ratio of 1: 2.5 (w/v). Soil organic carbon and total nitrogen were determined by the Walkley-Black wet digestion method and the Kjeldahl method, respectively (Jackson, 1956; Nelson and Sommers, 1982). Cation exchange capacity was measured by the ammonium acetate method (Sumner and Miller, 1996). Dithionite-citrate-bicarbonate soluble Fe and oxalate-soluble Fe were detected using inductively coupled plasma optical emission spectrometry (PerkinEimer's, Optima 8000, USA). Water stability of soil macroaggregates was determined using a modification of the wet sieving method (Nimmo and Perkins, 2002) and was expressed as mean weight diameter. Tensile strength of 10 macroaggregates was detected by a crushing test using a uniaxial load frame (SANS CMT0104, Shenzhen, China) for evaluating mechanical stability (Dexter and Watts, 2001). The measurement and calculation of tensile strength were described in detail by Zhou et al. (2012).

Macroaggregates (n = 5) were scanned using synchrotron-based micro-computed tomography at beam line BL13W1 at the Shanghai Synchrotron Radiation Facility. Detailed scanning settings and reconstruction procedures can be found in Zhou et al.'s (2013). The final 8-bit grayscale slices were stored in tiff format with a size of 1,024 × 1,024 × 1,024 voxels and a resolution of 3.7 μm. A cubic region of interest with a size of 600 × 600 × 600 voxels was selected from the central part of the soil macroaggregates for further analysis. Image preprocessing included removing the ring artifacts, enhancing the contrast, normalizing the brightness, and reducing noise with a 3D median filter (Zhou et al., 2013). The indicator kriging method was used to segment the voxels to pores and solids based on both the grayscale values and the spatial locations (Wang et al., 2011). The image-based porosity and pore size distribution of the segmented images were determined using the “opening” operation in the Quantim software (Vogel et al., 2010).

To investigate soil microbial community, soil samples were characterized by PLFA analysis following a modified method described by Frostegård et al. (1993). Briefly, total lipids were extracted from 2 g of freeze-dried soil samples with a chloroform-methanol-citrate buffer solvent (1:2:0.8, v/v/v) and partitioned into neutral, glyco- and phospho-lipids by a silica acid column. Phospholipids were subjected to mild alkaline methanolysis and converted into fatty acid methyl esters prior to analysis. The fatty acid methyl esters were quantified based on the addition of the internal standard methyl nonadecanoate (19:0) by a HP 6890 Series gas chromatograph instrument equipped with a 7683 Series injector using helium as a carrier gas (Hewlett Packard, Wilmington, Delaware, USA).

Identification was performed using bacterial fatty acid standards and a software from the Microbial Identification System (MIS; Microbial ID Inc., Newark, DE, USA). Individual PLFAs were named using the standard nomenclature. Total microbial biomass was calculated by summing the mass of all the detected fatty acids (Bossio et al., 1998) and expressed as nmol PLFA g⁻¹ dry soil. The sum of Gram-positive bacteria (GP, sum of all branched PLFAs) and Gram-negative bacteria (GN, sum of all monounsaturated PLFAs) was expressed as the total biomass (Frostegård and Bååth, 1996). From this value, the ratio of Gram-positive to Gram-negative bacteria (GP/GN) was calculated. The ratio of bacteria to fungi (B/F) was calculated by dividing the sum of bacterial PLFA markers by that of the fungal PLFA markers (18:1ω9 and 18:2ω6,9) (Zelles, 1997; Bååth and Anderson, 2003). The sum of PLFAs was further classified as actinomycetes (10Me PLFAs) (Frostegård et al., 1993) and AMF (16:1ω5) (Bach et al., 2010). Furthermore, two ratios of cyclopropyl PLFAs to their monoenoic precursors (cy17:0/16:1ω7 and cy19:0/18:1ω7) were used to indicate physiological or nutrient stress in the microbial communities (Kaur et al., 2005).

The capability of soil microbial communities to utilize a variety of carbon sources was measured with Biolog EcoPlate (Biolog Inc., USA) (Zak et al., 1994). The Biolog EcoPlate system consisted of 31 different carbon sources plus a blank well in three replicates. The carbon sources were subdivided into six substrate groups including carbohydrates, carboxylic acids, amino acids, polymers, phenolic acids, and amines (Weber and Legge, 2009). Briefly, soil microorganisms were extracted as follows: 5 g soil (dry weight equivalent) was added to 45 ml sterile 0.85% (w/v) saline solution, and the mixture was shaken for 30 min at 90 rpm. Each well of the Biolog EcoPlate was inoculated with 200 μl of the mixed suspension and incubated at 25°C in the dark for 7 days. The average well color development (AWCD) reflecting carbon utilization was recorded as the optical density value at 590 nm with a plate reader at regular 24-h intervals. As previously described, the absorbance of a single time point at 96 h was used to compare the wells.

Nematodes were extracted from soil macroaggregates using a modified Baermann funnel method (Barker, 1985). Briefly, a handful (100 g) of soil was placed on a double cotton filter (Hygia rapid, Hartmann AG, Heidenheim, Germany) over a wire screen. Soil samples were placed on top of a Baermann funnel, which was filled with distilled water up to 1/4 of the core height. The nematodes were allowed to move through the filter into the distilled water in the funnel for 2 days at room temperature, and then clean suspensions were ready for the counting and identification of nematodes. Four functional groups of nematode assemblages, including bacterivores, fungivores, plant parasites, and omnivores and predators, were identified based on known feeding habits or stoma and esophageal morphology (Yeates et al., 1993). Nematode populations were counted and expressed as the number of nematodes per 100 g of dry weight soil.

One-way analysis of variance (ANOVA) was used to test the effects of fertilization treatments on soil physicochemical properties, abundances of nematodes and microbial communities, and microbial carbon metabolism using Bonferroni's post-hoc test in the SPSS software (SPSS, Chicago, IL, USA). Principal coordinate analysis (PCoA) was performed to estimate the influence of fertilization on the Bray-Curtis distances of nematode and microbial community compositions (Anderson and Willis, 2003). We employed the “capscale” function for PCoA and “permutes” permutation-based testing for the calculation of the significance values.

Random forest modeling was conducted to quantitatively assess the important predictors of carbon metabolic capacity including soil properties, nematodes, and microbial community. The importance of each factor was evaluated by the increase in the mean square error between the observations and predictions (that is, the decrease in prediction accuracy) when the predictor was randomly permuted (Breiman, 2001). Random forest modeling was conducted using the randomForest package (Liaw and Wiener, 2002), and the significance of the model and predictor importance were determined using the A3R and rfPermute packages in the R software, respectively (Fortmann-Roe, 2013; Archer, 2016). Structural equation modeling was developed to uncover the direct and indirect contributions of soil physicochemical properties and nematodes to the microbial communities and carbon metabolic capacity. Structural equation modeling analysis was performed by the robust maximum likelihood estimation using AMOS 20.0 (Amos, Development Corporation, Meadville, PA, USA). The fitness of structural equation model was examined on the basis of nonsignificant chi-square test (P > 0.05), goodness of fit index, and the root mean square error of approximation (Hooper et al., 2008).

The results of one-way ANOVA revealed that fertilization treatments significantly affected the physical and chemical properties of soil [F_{(3, 8)} = 7.49–73.87, P < 0.01]. The proportion of soil macroaggregate fraction was significantly increased by 30.6–33.7% and 37.8–40.6% under the high level manure treatments (M2 and M3) compared with that under low manure (M1) and no manure (M0) treatments. The M2 and M3 treatments were characterized by significantly (P < 0.05) higher pH, SOC, total nitrogen, cation exchange capacity, and oxalate-soluble Fe but lower dithionite-citrate-bicarbonate soluble Fe than the M0 and M1 treatments (Supplementary Table 1). Tensile strength and mean weight diameter were significantly (P < 0.05) elevated by high manure application. In addition, distinguishable differences in the 3D pore structure were found under the four fertilization treatments (Figure 1). The image-based intra-aggregate porosity substantially declined under the high manure treatments (Supplementary Table 1, P < 0.05). The pore size distribution was found to vary greatly under manure treatments. The proportion of the pore classes >100 μm in soil macroaggregates was significantly higher under the M2 (42.2%) and M3 (43.7%) treatments compared with that under the M0 (15.7%) and M1 (26.9%) treatments, while the proportions of the 0–30, 30–60, and 60–100 μm pore classes followed the opposite trend (Figure 2). Microbial carbon metabolism determined by the AWCD was significantly (P < 0.001) distinguished by fertilization regimes, following a general trend of M2 > M3 > M1 > M0 (Figure 3). The utilization of six substrate groups, including carbohydrates, carboxylic acids, amino acids, polymers, phenolic acids, and amines, followed a trend similar to that shown over the whole plate (Figure 3, P < 0.01).

We determined the microbial biomass and community structure in soil macroaggregates by PLFA analysis. The total biomass varied considerably in response to manure treatments and was significantly higher under manure treatments (M1, M2, and M3) than that under the M0 treatment (Figure 4, P < 0.05). The general pattern of M2 > M3 > M1 > M0 was observed in specific PLFA signatures, such as bacteria, fungi, Gram-negative bacteria, Gram-positive bacteria, and actinomycetes (Figure 4). The distribution of microbial PLFA biomarkers as an index of community structure exhibited a substantial response to manure treatments. Principal coordinate analysis revealed that the microbial community structure under the high level of manure addition was distinguished clearly from that under the M0 and M1 treatments (Figure 4, P < 0.001). The ratio of B/F increased significantly, but the ratio of GP/GN decreased with the increase in manure addition (Supplementary Figure 1, P < 0.05). The ratios of cyclopropyl fatty acids to their precursors (17:0 cyclopropyl per precursor and 19:0 cyclopropyl per precursor) calculated as two stress indexes were in the microbial communities. The ratios of 17:0 cyclopropyl to 16:1ω7 (cy17:0/16:1ω7c) and 19:0 cyclopropyl to 18:1ω7 (cy19:0/18:1ω7c) under the high manure treatments were approximately 4.0–5.9-fold lower than those under the M0 treatment and 1.6–2.1-fold lower than those under the M1 treatment (Supplementary Figure 1).

The total number of nematodes in soil macroaggregates was counted under fertilization treatments, as well as four functional groups including bacterivores, plant parasites, fungivores, as well as omnivores and predators. Manure application considerably (P < 0.01) affected the number and composition of nematode assemblages. Bacterivores (46.6%) and plant parasites (30.5%) were the two most abundant trophic groups in the nematode assemblages (Supplementary Figure 2). Manure treatments significantly increased the total number of nematodes and the four functional groups in soil macroaggregates (Figure 5, P < 0.05). However, the elevated soil pH by lime reduced the total number of nematodes as well as bacterivores and plant parasites (Figure 5). Principal coordinate analysis revealed that the nematode assemblages were well segregated by fertilization treatments in regard to the abundances of bacterivores and plant parasites (Figure 5, P < 0.001). Plant parasites were enriched under the M1 treatment, comprising 73.5% of the total number, while bacterivores were dominated under the M2 and M3 treatments, comprising 57.1 and 69.9% of the total number, respectively (Supplementary Figure 2).

Bacterial and fungal biomass was positively associated with soil pH (r = 0.885, P < 0.001 and r = 0.665, P = 0.018) and SOC (r = 0.876, P < 0.001 and r = 0.598, P = 0.040), as well as microbial carbon metabolic function (r = 0.892, P < 0.001 and r = 0.771, P = 0.003). The proportion of pores with size >100 μm was significantly associated with bacterivores (r = 0.915, P < 0.001) and fungivores (r = 0.841, P = 0.001), as well as omnivores and predators (r = 0.881, P < 0.001). The abundance of AMF was significantly related to mean weight diameter (r = 0.641, P = 0.025). The abundance of bacterivores was closely correlated with bacterial biomass (r = 0.917, P < 0.001), Gram-negative bacteria (r = 0.960, P < 0.001), and microbial carbon metabolism (r = 0.947, P < 0.001). However, the abundance of bacterivores exhibited a marginal relationship with actinomycetes (r = 0.572, P = 0.052). The abundance of fungivores was highly correlated with AMF abundance (AMF, r = 0.693, P = 0.012) but showed no significant relationship with fungal biomass (r = 0.558, P = 0.059).

Random forest modeling was performed to identify the important predictors of nematode assemblages, microbial communities, and carbon metabolic activities. The models were significant at the 0.01 level with R² = 0.56–0.92. Random forest modeling indicated that soil pH and SOC were the two greatest determinants of microbial communities (12.8–27.2% and 11.7–15.5%) and carbon metabolism (15.7 and 10.8%) (Table 1). However, total porosity (13.0–21.3%) exhibited a stronger effect on nematode assemblages than SOC (8.7–19.1%) and pH (6.8–19.8%), followed by tensile strength (7.1–15.2%) and mean weight diameter (1.6–9.6%) (Table 1). Similarly, the predation of bacterivores suggested a dramatic impact on microbial carbon metabolism (6.2–13.8% and 8.4%) (Table 1). Structural equation modeling was used to assess the direct and indirect effects of soil physicochemical properties and nematode predation on the microbial community structure and carbon metabolic profiles. We found that the effect of soil chemical variables on carbon metabolic activities was higher than that of physical features (Figure 6). Microbial carbon metabolism was significantly related to microbial biomass rather than community structure. However, soil physical properties showed a primary effect on nematode assemblages, while chemical characteristics made a pronounced contribution to microbial biomass and community structure (Figure 6). Notably, nematode predation may exhibit an indirect positive influence on microbial carbon metabolism by promoting microbial biomass.

The biomarker PLFA analysis highlighted the significant differences in microbial biomass and community composition under different fertilization treatments. We observed that the spatial distribution of soil microbial biomass in macroaggregates clearly reflected the uniform pattern. The biomasses of total microorganisms and microbial functional groups were simultaneously enriched by the higher SOC and pH because of fresh organic matter inputs (Figure 4). Fertilization treatments play a key role in diversifying the microhabitats of soil macroaggregates; the chemical conditions (e.g., pH, SOC, and total nitrogen) of which represent an important aspect of niche differentiation that selects different functional groups (Briar et al., 2011; Jiang et al., 2018). Moreover, amorphous minerals have high surface areas to adsorb SOC, which affects the turnover rates of SOC pools and the availability of energy to soil microorganisms (Torn et al., 1997). The results suggested that these groups were sensitive indicators of the substantial increases in soil pH, SOC, and nutrient availability under fertilization treatments.

Manure treatments drastically altered the microbial community composition in terms of the ratios of B/F and GP/GN. The observed changes in soil microbial community structure revealed that the abundance of functional groups was unevenly distributed under management regimes (Supplementary Figure 1). More specifically, the increasing B/F ratio was largely governed by the faster increase in bacterial abundance than fungal abundance. The high bacterial abundance contributes prominently to SOC turnover and the liberation of nutrients and improves the overall quality (Smith et al., 2014). However, fungi play a crucial role in both the physical and chemical fertility of soil through aggregation and mycorrhizal associations, respectively. For instance, the fungal hyphal length is greater in aggregated soil than in sieved soil, emphasizing the interdependence between soil structure and microorganisms. Fungal hyphae have been shown to enmesh microaggregates into larger macroaggregates and facilitate the formation of macroaggregates (De Gryze et al., 2005). The extensive mycelial network of AMF produces large quantities of a sticky proteinaceous hydrophobic substance (glomalin-related soil protein) and recalcitrant polymers (chitin and chitosan) (Rillig and Mummey, 2006). Our results suggested that the increase in AMF abundance could improve the mean weight diameter, which was an integrated measure of soil aggregate structure and stability. The increase in AMF abundance suggests that the symbiosis of AMF with plants is gradually established and concomitantly beneficial for plant growth (Montesinos-Navarro et al., 2012). Furthermore, the decrease in the GP/GN ratio suggests that Gram-negative bacteria play a relatively important role in manured soils. Proteobacteria, including Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria, are enriched in the macroaggregates that typically contain alkyl and aliphatic groups (Janik et al., 2007). It is well established that closely interconnected Gram-negative bacteria (e.g., Burkholderia, Rhizobiales, Rhodospirillales, and Pseudomonas) frequently occupy a wide range of ecological niches and exhibit cooperative partnerships together with plants (Coenye and Vandamme, 2003; Jiang et al., 2017). The complex coupling patterns of these specific genera reflect synergistic relationships to alleviate abiotic stresses and promote plant growth in many crop systems (Bouizgarne et al., 2015). Notably, manure treatments decreased the ratios of cy17:0/16:1ω7c and cy19:0/18:1ω7c, suggesting that pH and nutritional stress were less limiting for the growth of the microbial population. Two ratios of cy per precursor have been proposed as indicators of starvation and pH stress to the microbial community (Moore-Kucera and Dick, 2008). Lower values of these ratios are associated with increased SOC supplementation and rapid rates of microbial growth and multiplication (Sampedro et al., 2009).

Manure treatments enhanced the total amount of nematodes in soil macroaggregates compared to the M0 treatment, which was principally caused by explosive increases in bacterivores under the M2 and M3 treatments (Figure 5). We observed that soil porosity exhibited a major effect on the structures of soil nematode assemblages (Figure 6, Table 1). Changes in soil nematode assemblages corresponded largely to those in soil porosity of the macroaggregates. The porosity of soil macroaggregates with porosity >100 μm was significantly improved under the M2 and M3 treatments (Figure 2), likely supporting the observed high densities of bacterivores. Nematodes are typically sensitive to the changes in the proportions of large pore sizes that are associated with intensive management (Cardoso et al., 2012). The intra-aggregate microarchitecture of soil macroaggregates with an adequate pore size range offers a suitable space for nematodes moving and migrating to their food sources. The percentage of soil macoraggregates with pore size >100 μm was an important determinant of the activity and survival of nematodes with large body sizes. Nematodes are aquatic organisms that depend primarily on thin water films to live and move through existing pathways of habitable soil pores with diameters of 30–90 μm (Hassink et al., 1993). However, the pores with diameters of 200–1,000 μm are much larger than the diameters of nematodes, which may constrain colonization and habitation by nematodes in soil macroaggregates (Quénéhervé and Chotte, 1996). Furthermore, the abundance of bacterivorous nematodes theoretically provides information about bacterial contributions to the decomposition and turnover of soil organic matter (Vestergård, 2004). The high amounts of bacterial-feeding nematodes under high manure treatments implied that the bacterial-based energy channel contributed strongly to SOC decomposition. In marked contrast, the increased pH with lime seemed to inhibit the total nematode density, which was particularly caused by a sharp decline in plant parasites under the M3 treatment (Figure 5). Ammonia has been reported as a nematicide to effectively suppress plant-parasitic nematodes with the application of organic manure and lime (Ben-Yephet et al., 2006). Acid-base neutralization at neutral pH values can trigger a shift in the equilibrium from ammonium toward ammonia and enhance the nematicidal activities of controlling nematodes (Oka et al., 2007). Despite the moderate influence of pH, the abundance of plant parasites would plummet in response to neutral conditions as a result of the increased nematicidal activity (Oka et al., 2006; Jiang et al., 2013).

Our study took a further step toward describing the intimate interactions between nematodes and microorganisms in macroaggregates. The high value of the nematode structure index in macroaggregates suggests a complex community structure with many linkages in soil food web (Jiang et al., 2013). The positive bacterivore-bacteria relationship suggested that predation by bacterivores stimulates the increase in bacterial biomass. Recently, the incubation experiments have reported that the effects of nematode grazers on microbial biomass range from positive to negative (Trap et al., 2016). The slight to moderate grazing of bacterivores exerts positive feedback on bacterial biomass and activity (Fu et al., 2005). However, the population developments of each nematode species differed according to the prey bacteria. This study suggested that bacterivores showed strongly positive relationships with total bacteria and Gram-negative bacteria rather than actinomycetes. Bacterivores do not multiply in the presence of the filamentous bacteria Actinomyces sp. and the Gram-positive bacteria Arthrobacter sp. (Buchan et al., 2013). Large and filamentous bacterial cells can escape uptake by nematodes as a result of their small buccal cavity (Venette and Ferris, 1998). Gram-positive bacteria are likely to be less suitable food for protozoans and nematodes than Gram-negative bacteria. The lower edibility of Gram-positive bacteria may be related to a lower rate of digestion of their cell wall, which enables survival during passage through microbial grazers (Rønn et al., 2012).

Nematodes hold a central position in soil food webs, and the population dynamics of nematode communities are closely linked to SOC dynamics. The presence of bacterivorous nematodes substantially promotes microbial basal respiration contributing to SOC dynamics (Neher, 2010). Structural equation modeling quantitatively described that bacterivores may exert a stronger impact on microbial biomass than community structure and eventually increase microbial carbon metabolism in the macroaggregates (Figure 6). Our previous result reported that the predation of bacterivores exhibited indirectly positive associations with the sizes and turnover rates of SOC pools in the macroaggregates (Jiang et al., 2018). Microorganisms can dramatically increase SOC retention by directly incorporating living and senescing fungal and bacterial biomass into the stable SOC reservoir (Liang et al., 2017). Moreover, selective feeding by bacterivores on active bacteria suppresses soil metabolic quotient and thus stimulates SOC accumulation in the macroaggregates (Jiang et al., 2013). Therefore, this study supported that increased bacterivore abundance enhanced the grazing pressure on microorganisms and concurrently promoted microbial-derived SOC turnover. In addition, AMF generally prefer to grow in the macroaggregates with high porosity to facilitate SOC stabilization and protection by enhancing well-developed soil aggregation. The contribution of AMF hyphae to SOC sequestration contains extra-radical hyphae as well as glomalin-related soil proteins produced by the hyphae (González-Chávez et al., 2010). Glomalin is considered an important pool of recalcitrant SOC with a relatively long turnover time (Driver et al., 2005). Furthermore, AMF that form mutualistic associations with plants are vulnerable to grazing by fungivorous nematodes. Fungivores Aphelenchoides sp. populations grow rapidly on saprophytic (Agrocybe and Chaetomium) and mycorrhizal (Cenococcum, Hymenoscyphus, and Laccaria) fungi (Ruess et al., 2000). Furthermore, the predation by fungivorous nematodes on AMF Glomus spp. and Acaulospora spp. shows a positive feedback effect on AMF abundance (Hua et al., 2010). Collectively, soil macroaggregates could be identified as a well-organized soil food web with interrelated functional relationships between nematodes and microorganisms in regulating SOC turnover.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.