Complex soil food web enhances the association between N mineralization and soybean yield – a model study from long-term application of a conservation tillage system in a black soil of Northeast China

Long-term (10 years) application of conservation tillage following conversion from conventional tillage (CT) can achieve a new equilibrium in the soil environment, which is vital to reverse soil biodiversity declines and fulfil the goal of maintaining agroecosystem sustainability. However, in such a situation, how the soil community regulates nutrient cycling impacting crop yield is not well documented. Therefore, the relations between mineralized nitrogen (N) delivered by soil food web and soybean (Glycine max Merr.) yield were investigated after 14 years application of CT, reduced tillage (RT) and no tillage (NT) in a black soil (Typic Hapludoll) of Northeast China. We hypothesized that soil mineralizable N would increase with the complexity of the soil food web, and that the trophic groups involved in associating N mineralization with crop yield will vary with soil depth in the conservation tillage practice. During the soybean growing season, soil organisms, including bacteria, fungi, nematodes, mites and collembolans, were extracted and identified monthly from 0–5 and 5–15 cm soil depths to estimate the complexity of the food web indicated by the species richness and connectance indices, and to simulate the mineralized N using energetic food web modelling. The species richness and connectance of the food web at both soil depths were significantly affected by tillage practices, and their values decreased of the order of NT > RT > CT. A similar trend was also revealed for the simulated N mineralization, that is, the mineralized N released either from the functional feeding guilds or from the energy pathways of the food web were greater in RT and NT than in CT at both soil depths. Multiple linear regression analysis showed that soil organisms involved in coupling the mineralized N with soybean yield were different at different soil depths, in which fungal and root pathways at 0–5 cm and bacterial pathway at 5–15 cm were the driving factors for the supply of mineralized N to soybean in NT and RT soils. These results support our hypothesis and highlight the essential role of soil food web complexity in coupling N mineralization and crop yield after long-term application of conservation tillage. Additionally, the current modelling work provides basic hypotheses for future studies to test the impact of soil biodiversity or specific functional guilds on the fate of N in agro-ecosystems.

1 Introduction

Nitrogen (N) is the most important growth-limiting nutrient for crops (Fageria et al., 2010). In order to achieve the maximum yield, N fertilizer is applied to crops all over the world; even legumes that fix N through symbiotic N-fixing microorganisms require additional chemical N application for maximum yield (La Menza et al., 2020). However, globally, the N recovery rate by crops is only 60 % at most (Liu et al., 2010), which means that the rest of the fertilizer N is not available for the crop and is lost by leaching, runoff, ammonia volatilization or nitrous oxide emission, resulting in undesirable environmental issues. Hence, crop production, to a great, extent relies on N mineralization to meet the growth requirements (La Menza et al., 2020; Whalen et al., 2013).

The process of N mineralization mediated by soil organisms is closely related to the predation across multi-trophic groups because soil organisms require carbon (C), N and other nutrients from the prey to support their metabolic activities, ultimately converting the organic N compounds into the form of mineral N (de Ruiter et al., 1993; Whalen et al., 2013). The N immobilized in the biomass of the lower trophic groups can be released by the predation of the higher trophic groups. Furthermore, the predators usually have a higher C : N ratio than their prey, which results in more N obtained than their nutritional requirements, and the excess N is excreted into the soil ammonium (NH${}_{\mathrm{4}}^{+}$) pool (de Ruiter et al., 1993; Whalen et al., 2013). It is estimated that the N amount released from the predation of soil organisms accounts for 30 %–80 % of the annual N mineralization under field conditions (Carrillo et al., 2016; de Ruiter et al., 1993; Holtkamp et al., 2011), and the value of this contribution varies with the complexity of soil food webs (Carrillo et al., 2016; de Ruiter et al., 1993; Holtkamp et al., 2011). Several studies (Bender and van der Heijden, 2015; Thakur et al., 2014; Wagg et al., 2014) based on controlled (micro- or mesocosm) experiments demonstrated that the potential mineralizable N pool increases with the increase in complexity of the food web, which implies that a management practice that forms a complex soil food web is beneficial for improving N availability with less N fertilizer input.

Conservation tillage, including reduced tillage (RT) and no tillage (NT) with at least 30 % mulch cover of the soil surface, is becoming a popular practice around the world to counteract the disadvantage of conventional tillage (CT; soil inversion by mouldboard ploughing) on soil health (Lal, 2004). The benefits of conservation tillage on soil properties and processes, especially on crop productivity, are not immediately apparent, but can only be achieved after a period of time (5–10 years) when the soil environment reaches a new, stable equilibrium (Six et al., 2004). In such stable situations, soil biodiversity and its spatial heterogeneity are strongly enhanced, thereby constructing a more complex network among soil organisms relative to CT (de Vries et al., 2013; D'Hose et al., 2018; van Capelle et al., 2012). For example, bacteria and bacterivorous fauna dominate the whole plough layer of CT, while conservation tillage is typically characterized by fungi and fungivorous fauna near the surface and bacterial-based communities at deeper soil depths (D'Hose et al., 2018; van Capelle et al., 2012). Moreover, the increase in the richness and density of predaceous fauna reorganizes the topological structure of the food web through the modification of the prey–predator interactions (Bartley et al., 2019). Our understanding of how the entire food web assemblages mediate N mineralization to maintain crop yields after long-term conservation tillage is still limited.

Soybean (Glycine max Merr.) is a major crop produced in the black soil region of Northeast China and accounts for 50 % of the total national soybean production (Liu et al., 2019). Monoculture cropping, mouldboard ploughing, ridging, seeding into ridges and the removal of post-harvest residue is the typical practice in this region, which has caused serious land degradation that is threatening soil fertility and sustainability (Zhang et al., 2019). Consequently, a national project to promote the application of conservation tillage in Northeast China was launched in 2020 (Ministry of Agriculture of China and Ministry of Financing of China, MAC and MFC, 2020). Reliable information regarding the response of soil properties and grain yield to the conversion from CT to conservation tillage is needed to help the farmers adopt better agronomic measures.

The objective of this study was to investigate the relations between N mineralization by the soil food web and soybean yield under a long-term conservation tillage system. We hypothesized that (1) conservation tillage favours a greater release of mineralized N than CT as it forms a more complex soil food web, and (2) the trophic groups of soil organisms associated with N mineralization and crop yield will vary with soil depth in the conservation tillage system, given the strengthened heterogeneity of organisms along the soil profile.

To address these hypotheses, soil organisms, including bacteria, fungi, nematodes, mites and collembolans, were extracted monthly during the soybean growing season after a long-term (2001–2015) application of conventional tillage (CT), reduced tillage (RT) and no tillage (NT) in a black soil of Northeast China. The amount of mineralized N delivered by all trophic groups in the food web was calculated using the energetic food web modelling approach. This approach has been applied to a range of natural and agricultural systems and has proven very useful in simulating N mineralization and in understanding the ecological functions served by soil organisms (Barnes et al., 2014; Carrillo et al., 2016; Koltz et al., 2018; Pressler et al., 2017; Schwarz et al., 2017), although it cannot reflect the dynamics N flow in the same way as the isotope tracing technique.

2 Material and methods

2.1 Experimental design and soil sampling

This study was conducted at the Experimental Station (44^∘12^′ N, 125^∘33^′ E) of the Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, in Dehui County, Jilin Province, China. The station is located in a continental temperate monsoon zone. The mean annual temperature is 4.4 ^∘C; the lowest mean monthly temperature occurs in January (−21 ^∘C) and the highest in July (23 ^∘C). The mean annual precipitation is 520 mm, and >70 % occurs from June to August. The soil is classified as black soil (Typic Hapludoll; USDA Soil Taxonomy) with a clay loam texture (the average soil texture is 36.0 % clay, 24.5 % silt and 39.5 % sand).

The present experiment was conducted as part of an ongoing long-term tillage and crop rotation experiment. The long-term tillage experiment was established in the fall of 2001 and included conventional tillage (CT), reduced tillage (RT) and no tillage (NT) in a 2-year maize (Zea mays L.) and soybean (Glycine max Merr.) rotation system with residue return. Each treatment had four replicates, and the plot area was 5.2×20 m. Crops were sown at the end of April or early May and harvested in October every year and then fallowed for 6 months over the winter when the soil was frozen. The CT treatment consisted of fall mouldboard ploughing (20 cm), followed by a secondary seedbed preparation in spring by discing (7.5–10 cm), harrowing and ridge-building. In RT, the ridges (16 cm height and 75 cm width) were rebuilt with a cultivator in June of each year; a modified lister and scrubber were used to form and press the ridges. The soil of NT had no disturbance except for planting using a no-till planter (KINZE-3000NT; Kinze Manufacturing, Williamsburg, IA, USA). After harvest, the aboveground residues were returned to the soil surface in all treatments to prevent water and wind erosion in winter and early spring (Zhang et al., 2019). For RT and NT plots, maize residue was cut into about 30 cm pieces, leaving a 30–35 cm standing stubble; soybean residue was directly returned to the soil surface. Residues in CT plots were removed prior to the fall mouldboard ploughing, manually replaced on the soil surface after fall mouldboard ploughing and then mixed with the plough layer by discing and cultivation in the following spring.

Starter fertilizer was applied with the planter at a rate of 89 kg N ha⁻¹, 51 kg P ha⁻¹ and 51 kg K ha⁻¹ for maize and 40 kg N ha⁻¹, 49 kg P ha⁻¹ and 53 kg K ha⁻¹ for soybean. Additionally, 45 kg ha⁻¹ of N was top dressed at the V6 stage (six leaves) of maize. The application rates of N, P and K were the same in all tillage treatments, and the N application rate was reduced by about 30 % compared to the local conventional application rate (187 and 60 kg N ha⁻¹ for maize and soybean, respectively).

2.2 Soil sampling

Soils were sampled for the present experiment at the end of each month from April to September 2015 during the soybean growing season. The total precipitation during the growing season was 365 mm in 2015, which was located in the range of 330–605 mm across the past 10 years (2004–2014), and there were no typhoons in 2015 (data obtained from China Meteorological Data Service Centre; http://data.cma.cn/en; last access: 28 January 2021). All types of soil organisms, including microbes, nematodes and micro-arthropods, were determined monthly, except for nematodes, which were only determined in April, June and August due to a limited labour force. The nematode populations for non-sampled months were estimated by linear interpolation between adjacent sampling dates.

A total of seven soil cores (2.5 cm in diameter) in each plot were randomly collected from a depth of 15 cm and each core was separated into 0–5 and 5–15 cm sections. Soil cores were combined to form a single composite sample for each plot and depth. Samples were immediately taken to the lab and stored at 4 ^∘C. Soil bulk density for each plot was determined in the 0–5 and 5–15 cm depths, using a slide hammer probe with a 5 cm core diameter. The mean monthly bulk density is presented in Table S1 in the Supplement (hereafter, “S” refers to the Supplement). After the plants had reached physiological maturity, the soybean yield in 2015 was determined by hand-harvesting 3 m lengths of six interior rows from each plot. Grain yield samples were dried to a constant weight at 75 ^∘C in an oven and then corrected to 13.5 % grain moisture content.

2.3 Soil mineral N in the field condition

The content of soil mineral N, determined by summing NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$, in the field condition, was measured within 12 h after soil samples were collected each month. Mineral N was extracted by 1 M KCl (soil : KCl = 1 : 2 $w/v$) and determined by a continuous flow analyser (SAN$++$; Skalar Analytical B.V., the Netherlands).

2.4 Soil organism extraction

The microbial community was determined using the phospholipid fatty acid analysis (PLFA) as described by Bossio et al. (1998). Lipids were extracted from 8 g of freeze-dried soil with a Bligh and Dyer solution (chloroform : methanol : citrate buffer = 1 : 2 : 0.8 ($v:v:v$)). Polar lipids were separated from neutral lipids and glycolipids in a solid phase extraction column (Supelco Inc., Bellefonte, PA, USA) and transformed into fatty acid methyl esters with a mild alkaline methanolysis. Samples were then dissolved in hexane and analysed in an Agilent 6850 series gas chromatograph with MIDI peak identification software (Version 4.5; MIDI, Inc., Newark, DE, USA). Fatty acids were grouped as bacteria (14 : 0, i14 : 0, a14 : 0, 15 : 0, i15 : 0, a15 : 0, 15 : 1ω6c, 16 : 0, i16 : 0, a16 : 0, 16 : 1ω7c, 16 : 1ω9c, i17 : 0, a17 : 0, 17 : 1ω8c, 17 : 1ω9c, 18 : 1ω7c, 18 : 0 and 20 : 0), saprophytic fungi (18 : 1ω9c and 18 : 2ω6c) and arbuscular mycorrhizal fungi (AMF; 16 : 1ω5c; Bach et al., 2010; Dempsey et al., 2013). Microbial biomass was estimated using the following conversion factors of fatty acid concentrations (nmol): bacterial biomass, 363.6 nmol = 1 mg C, saprophytic fungal biomass, 11.8 nmol = 1 mg C, and AMF biomass, 1047 nmol = 1 mg C (Tsiafouli et al., 2015).

Nematodes were extracted from a 50 g soil sample (fresh weight) using a modified cotton wool filter method (Liang et al., 2009). At least 100 nematode specimens from each sample were selected randomly and identified to genus level (see Table S2 for the list of identified taxa and Table S4 for the abundance) using an Olympus BX51 microscope (Olympus Corporation, Tokyo, Japan), according to Bongers (1994). Nematodes were assigned into four trophic groups, namely bacterivores, fungivores, plant parasites and omnivores and predators (Ferris, 2010). Body length and maximum body diameter of nematodes were measured using an ocular micrometer to calculate the nematode fresh body mass (µg; Andrássy, 1956). Nematode biomass was estimated by assuming that the dry weight of a nematode is 20 % of the fresh weight, and the C in the body is 52 % of the dry weight (Ferris, 2010).

Microarthropods were extracted from 200 mL fresh soil using modified high-gradient Tullgren funnels (Crossley and Blair, 1991) for 120 h. Individuals were collected and stored in vials containing 95 % ethanol for identification. Mites and collembolans were identified to species or morphospecies level (see Table S3 for the list of identified taxa and Table S4 for the abundance), according to Christiansen and Bellinger (1980–1981), Balogh and Balogh (1992), Bellinger et al. (2019), Pomorski (1998) and Niedbala (2002). Soil microarthropods were allocated into four different functional groups, namely fungivorous (oribatid) mites, predaceous mites, fungivorous collembolans and omnivorous collembolans. Individual body length and width were measured to estimate the dry weight based on regression equations from the literature (Douce, 1976; Hódar, 1996). Mite and collembolan biomass were estimated by assuming the C in the body as being 50 % of the dry weight (Berg, 2001).

2.5 Modeling N mineralization by the food web

The first step in modelling the N mineralization by the food web was to construct a soil food web using the published feeding relationships (Fig. S1). All determined species of soil micro-flora and micro-fauna were grouped into six functional feeding guilds based on the trophic resources they exploit (Burns, 1989), namely bacteria, fungi, herbivorous feeders, bacterivorous feeders, fungivorous feeders and predaceous feeders. For omnivorous/predaceous species, we constructed every possible predator–prey interaction. Omnivorous/predaceous nematodes were assumed to feed on all other nematode groups (Yeates et al., 1993). Omnivorous collembolans, which mainly feed on bacteria, fungi, plant and microfauna (de Vries et al., 2013), were proportionally assigned to bacterivorous, fungivorous, herbivorous and predaceous collembolans according to the assumption that their diet consists of 25 % bacteria, 25 % fungi, 25 % plant and 25 % other microfauna. Taking into account the changes in abundance of soil organisms over time, the biomass during the soybean growing season was estimated by summing the monthly biomass. The biomass unit of each functional guild was converted from mg C g⁻¹ to mg C m⁻², based the soil bulk density (g cm⁻³) and thickness of the soil layer (15 cm). Furthermore, the complexity of soil food web was measured by the species richness (the number of taxa detected in the sample) and the connectance (expressed as the ratio of the number of actual links to the total number of possible feeding links) indices (Zhang et al., 2015).

The N mineralization was simulated using the energetic food web model based on a mass balance assumption that the energy flowing into a group in the food web is equal to the energy flowing out through metabolism and predation (Barnes et al., 2014; de Ruiter et al., 1993). The following equations were used to simulate the N mineralization delivered by soil organisms according to de Ruiter et al. (1993):

$$\begin{array}{}\text{(1)}& {\displaystyle }{F}_{ij}={\displaystyle \frac{w_{ij}{B}_i}{\sum _{k=\mathrm{1}}^n{w}_{kj}{B}_k}}\text{(2)}& {\displaystyle }F={\displaystyle \frac{d_j{B}_j+P_j}{e_{\mathrm{ass}}\times e_{\mathrm{prod}}}}\times F_{ij}\text{(3)}& {\displaystyle }{N}_{min}=e_{\mathrm{ass}}\times \left({\displaystyle \frac{\mathrm{1}}{{\mathrm{C}:\mathrm{N}}_i}}-{\displaystyle \frac{e_{\mathrm{prod}}}{{\mathrm{C}:\mathrm{N}}_j}}\right)\times F,\end{array}$$

where, in Eq. (1), F_ij is the feeding preference of predator (j) on prey (i), which was calculated based on the density-independent feeding preference of j on i (w_ij; dimensionless; listed in Table S5), n is the total number of potential prey types ($k=\mathrm{1},\mathrm{2},\mathrm{3}\mathrm{\dots }n$), and B is the biomass of prey (mg C m⁻²). In Eq. (2), F is the feeding rate of predator on prey (mg C m⁻² yr⁻¹), d_j is the natural death rate of j (yr⁻¹), B_j is the biomass of j (mg C m⁻²), P_j is the energy loss of j due to the predation (mg C m⁻² yr⁻¹), and e_ass and e_prod are the assimilation efficiency and production efficiency of j, respectively. In Eq. (3), N_min is the N mineralization mediated by the predation of j on i (mg N m⁻² yr⁻¹); C : N_i and C : N_j is the body C : N ratio of prey (i) and predator (j), respectively. The parameters of d, e_ass, e_prod and C : N of soil organisms were taken from the literature and are presented in Table S6. It is impossible to measure and confirm each parameter value under field conditions; therefore, these parameters were cited from the relevant studies that were also conducted on an agricultural system and updated according to the latest reports.

The simulation of N mineralization was started with the top predators, which are considered to have no energy loss from the predation, and then proceeded to the lower trophic groups. Based on the specific primary actors that drive energy flow from the basal resource to the food web, the soil food web was further categorized into the following three energy pathways: (1) the fungal pathway, in which the energy flux is driven by fungi and then flows to fungivores and their predators, (2) the bacterial pathway, in which the energy flux is driven by bacteria and then flows to bacterivores and their predators, and (3) the root pathway, in which energy flux is driven by herbivores and then flows to their predators. The N mineralization was first estimated at the functional guild level, by summing up the contribution of all species within a functional guild, and then estimated at the level of each energy pathway (bacterial, fungal and root pathways) and, finally, estimated for the whole food web.

2.6 Statistical analyses

In our initial model, the omnivorous collembolans were assigned into bacterivores, fungivores, herbivores and predators in equal portions to model the mineral N flux within the soil food web. To assess the impact of this assumption affecting mineral N flux, a sensitivity analysis was performed by re-assigning omnivorous collembolans into fungivores and herbivores (50 % each), according to Barnes et al. (2014). This acted as a null model with the least diet preference, and the difference from the initial model was expressed as a percentage.

Data were checked for normality and for homogeneity in the variances prior to statistical analysis. If necessary, the data were ln (x+1) transformed to meet the assumptions of an analysis of variance (ANOVA). A two-way ANOVA was performed to test the effect of tillage, soil depth and their interaction on the empirically observed soil mineral N, the biomass of each feeding guild, the complexity index and the simulated N mineralization of soil food webs. When their interaction was significant, multiple comparisons were performed based on a post hoc test to determine if tillage effects were significant in each soil depth. Tukey's honestly significant difference test was used for means comparisons, and a difference at the P<0.05 level was considered statistically significant.

Forward stepwise multiple linear regression (MLR) was used to identify which energy pathways closely link the release of mineralized N to soybean yield at each soil depth. In a stepwise regression, only one independent variable is considered at a time, and another variable is added to the model at each step until no significant (P value was set at 0.05) improvement in the percentage of explained variance is obtained. Prior to MLR, all parameters were min–max normalized to accurately preserve all relations of the data value and prevent potential bias from the domination of variables with large numeric ranges over those with small numeric ranges. Min–max normalization subtracted the minimum value of an attribute from each value of the attribute and then divided the difference by the range of the attribute. The normalized data lay in the range [0, 1] (Jayalakshmi and Santhakumaran, 2011). All statistical analyses were performed in R software (R 3.4.0; R Core Team, 2017) with the package called car for ANOVAs and the package called stats for MLR analyses.

3 Results

3.1 Soil mineral N and soybean yield

Tillage effect on the soil mineral N over the whole soybean growing season varied with soil depth (Table 1). At 0–5 cm, the amount of soil mineral N was higher (P=0.001) in RT than in CT, while an opposite trend was observed at 5–15 cm, with a lower (P<0.001) amount in RT and NT than in CT. For the entire soil layer (0–15 cm), NT significantly (P=0.027) decreased the amount of soil mineral N relative to CT. There was no statistical significance (P=0.065) for differences in soybean yield among tillage treatments (Table 1), although there was a general trend of NT > RT > CT.

Table 1 Effects of tillage systems on soybean yield and the cumulative mineral N concentrations (standard deviation – SD) during the soybean growing season.

^∗ Tillage effect was tested by one-way ANOVA. CT – conventional tillage; RT – reduced tillage; NT – no tillage. The same lowercase letter in the same row indicates no significant difference among tillage systems (P>0.05).

3.2 Metrics of soil food web

Tillage significantly influenced the complexity of the soil food web, as indicated by the indices of species richness and connectance (Fig. 1). Compared with CT, the whole food web richness (P=0.035) and connectance (P=0.001) were significantly increased in NT at both soil depths, while only connectance (P=0.045) was significantly increased in RT.

Compared to CT, NT and RT significantly increased the biomass of the whole food web by 33 %–56 % at 0–5 cm and by 28 %–42 % at 5–15 cm (Table 2). A similar trend was also found for the functional feeding guilds of bacteria, fungi, herbivores, bacterivores, fungivores and predators at both soil depths, with a higher biomass in RT and NT than in CT (P<0.05; Table 2).

Figure 1 Changes in soil food web complexity (indicated by richness, a, and connectance, b) under different tillage systems. CT – conventional tillage; RT – reduced tillage; NT – no tillage.

Table 2 Cumulative soil biotic biomass (expressed as mg C m⁻²) over the soybean growing season under different tillage practices (standard deviation – SD).

CT – conventional tillage; RT – reduced tillage; NT – no tillage; ns – no significant difference (P>0.05). Means for the different tillage systems at the same depth, and followed by the same lowercase letter, are not significantly different (P>0.05).

3.3 Mineralized N modelled by soil food web

To visualize the mineralized N within the food web, an N flux network calculated by the energetic food web model was constructed for different tillage systems at both soil depths (Fig. 2). When the omnivorous collembolans assigned in equal portions into four classes, namely bacterivorous, fungivorous, herbivorous and predaceous collembolans, were replaced with two classes, namely fungivorous and herbivorous collembolans in equal portions, there was a slight decrease in the total amount of mineralized N within the food web throughout all tillage systems (loss of 0.3 % at 0–5 cm and 2 % at 5–15 cm; Fig. 3). When the two functional feeding guilds were compared to four, a sharp decrease in the mineralized N from bacterivores, fungivores and herbivores to the top predators was observed for all tillage systems (decreasing 22 %–83 % at 0–5 cm and 2 %–24 % at 5–15 cm), although there was an increase in the mineralized N from roots to herbivores and from fungi to fungivores.

Figure 2 Modelled mineral N flux (mg N m⁻² yr⁻¹) through the fungal, bacterial and root pathways in the food webs of different tillage systems during the soybean growing season. Fungal, bacterial and root pathways are the energy fluxes driven on the base of fungal, bacterial and herbivorous communities. Vertical width of bars represents the amount of mineral N transferred between the trophic levels. CT – conventional tillage; RT – reduced tillage; NT – no tillage.

Figure 3 Differences (in percent) in the modelled mineral N delivered by the food web between model 1 that assigned the omnivorous collembolans into bacterivorous, fungivorous, herbivorous and predaceous collembolans in equal portions and model 2 that assigned the omnivorous collembolans into fungivorous and herbivorous collembolans in equal portions. CT – conventional tillage; RT – reduced tillage; NT – no tillage.

Across the entire soybean growing season, RT and NT significantly (P<0.001) increased the total amount of mineralized N within the food web by 33 %–41 % at 0–5 cm and 28 %–38 % at 5–15 cm relative to CT, and the maximum increase was observed in NT at both soil depths (Fig. 2 and Table S7). Furthermore, the mineralized N delivered by the energy pathways also varied with tillage systems (Fig. 2). Compared to CT, RT and NT significantly (P<0.001) increased the amount of mineralized N in the bacterial and fungal pathways at both soil depths, while only NT significantly (P=0.001) increased the amount of mineralized N in the root pathway at 0–5 cm (Table S7). The similar tillage effect pattern was also observed for the components within these energy pathways (Fig. 2 and Table S7). Compared with CT, NT significantly (P<0.05) increased the mineralized N released from each component in the fungal and bacterial pathways, while RT only significantly (P<0.05) increased the mineralized N from the basal resource to the primary decomposers and then to the intermediate microbial feeding fauna (fungivores and bacterivores) at both soil depths. For the components in the root pathway at both soil depths, a greater quantity of mineralized N from the basal resource to herbivores was released in RT and NT than in CT (P<0.001), and the mineralized N from herbivores to predators revealed no significant (P>0.05) difference among CT, RT and NT.

3.4 Relation between mineralized N in the food web and soybean yield

The multiple linear regression model (Table 3) showed that 83.6 % of the variation in soybean yield was explained by the mineralized N released from fungal and root pathways at 0–5 cm. Their relative contributions to the soybean yield decreased of the order of the fungal pathway (0.557) > root pathway (0.550), which means that when the min–max normalized N mineralization in the fungal pathway and root pathway increases by one, the min–max normalized soybean production would correspondingly increase by 0.557 and 0.550 respectively. At 5–15 cm, only the mineralized N delivered by the bacterial pathway significantly affected soybean yield and accounted for 37.3 % of the yield variance. Soybean yield would increase by 0.656 units when the bacterial pathway increased by one.

Table 3 Relations between modelled mineral N delivered by different energy channels and soybean yield based on multiple linear regression analysis. Data were min–max normalized and are dimensionless.

The ^∗ and ${}^{\ast \ast }$ indicate significance at 0.05 and 0.01, respectively.

4 Discussion

Soil N availability, which is generally linear with the crop yield, highly depends on the pool of soil mineral N and mineralizable N regulated by soil organisms (Fageria et al., 2010; Whalen et al., 2013). In this study, we monitored the variation in soil mineral N and modelled the amount of mineralizable N within the food web throughout the whole growing season of soybean under different tillage systems. The results showed that the variation pattern in the soybean yield among different tillage systems is counter to the empirically observed soil mineral N (Table 1), which was lower in RT and NT than in CT, at both the lower soil layer (5–15 cm) and at the entire layer (0–15 cm), but is consistent with the simulation of mineralizable N pool (Fig. 2 and Table S7) that decreased of the order of NT > RT > CT at both soil depths. Consequently, the mineralizable N pool has a greater contribution than the soil mineral N pool over the growing season for the soybean yield in RT and NT relative to CT; a detailed discussion is presented in Sect. 4.1 below.

We acknowledge that this study was based on only 1 year of data collection from a continental climate region, and therefore, due to the high variability of soil organisms in response to external disturbances, our work may not be directly applicable to other climate regions in the world that have also adopted conservation tillage system. Nevertheless, the current model work highlights the importance of soil food web complexity in coupling N mineralization and crop yield after long-term application of a conservation tillage system and can serve as hypotheses for future studies to test the impact of soil biodiversity or specific functional guilds on the fate of N in agro-ecosystems.

4.1 Performance of modelling N mineralization within the food web

A source of uncertainty in the simulation of mineralized N was the feeding preference assignment of omnivorous collembolans, which were allocated into equal portions of bacterivores, fungivores, herbivores and predaceous. The robustness of this assumption was tested using a sensitivity analysis by re-assigning the omnivorous collembolans into fungivores and herbivores (50 % each), which resulted in, at most, 3 % loss in the mineralized N of the whole food web. However, when considering the two functional feeding guilds, there was a dramatic decline in mineral N from bacterivores, fungivores and herbivores to the top predators (Fig. 3). These results indicate that the disparity between these two models highly depends on the feeding guilds, and accordingly, the assignment of the species into the functional guilds should be done with caution.

To our knowledge, there is no literature to date that has comprehensively identified the feeding habits of collembolans because they consume a wide spectrum of resources, including plant roots or litter, different types of soil microbes and metazoan soil fauna (Potapov et al., 2016). Additionally, collembolans can shift their diet from one food resource to another when choices are available (Chahartaghi et al., 2005; Endlweber et al., 2009). This inherently complex feeding nature of collembolans makes it difficult to correctly assign them to specific feeding guilds without using isotope tracer techniques. In this study, the diet of collembolans may change throughout the crop year, according to the availability of basal food resources of growing plants and crop residue, and organisms higher up in the food web in the different tillage systems. Therefore, except for those specific species that feed on fungi, other species classified as omnivorous collembolans (Table S3) are reasonably treated as generalists. In summary, our presented model is robust in calculating the N flux within the food webs under different tillage systems over the soybean growing season.

4.2 Tillage effects on the N mineralization within the food web

Consistent with our first hypothesis, the results showed that, as the structure of soil food web became more complex after the conversion from conventional tillage to conservation tillage, mineralized N released either from the functional feeding guilds or from the energy pathways of the food web was greater in RT and NT than in CT at both soil depths. Our result is in agreement with the reports of Bender et al. (2015), Carrillo et al. (2016) and de Vries et al. (2013) which state that farming practices favouring rich and abundant soil organisms can increase N availability. This may be due to the increase in the number of different kinds of species, leading soil organisms to release more N when they consume basal resources to create their own biomass (Holtkamp et al., 2011; Koltz et al., 2018). The higher biomass of the food web revealed in RT and NT than in CT (Table 2) further supports our results as more N would be released with the build-up of biomass. Additionally, the strengthened connectance between functional guilds (Fig. 1) also contributes to the increase in the amount of simulated N mineralization in RT and NT because the tight interlinkage within trophic levels in the food web stimulates the N release from predation (Bender et al., 2015; Carrillo et al., 2016; Wagg et al., 2014). Therefore, after 14 years continuous application of conservation tillage, a large variety of organisms and complex interlinks among them expand the potential mineralizable N pool.

4.3 Relations between N mineralization within the food web and soybean yield

Although RT and NT improved the amount of mineralized N within the food web over the soybean growing season, the multiple linear regression analysis showed that soil organisms involved in coupling the mineralized N with the soybean yield were different, along with the soil profile (Table 3). We found that the mineralized N released from fungal and root pathways was strongly related to the soybean yield at the surface of 0–5 cm, while at the 5–15 cm depth, only the mineralized N released from the bacterial pathway significantly contributed to the yield. These results strongly support our second hypothesis that the trophic groups of soil organisms responsible for N mineralization associated with crop yield vary with soil depth.

Fungal and bacterial pathways have been considered as being two very important parallel pathways in mediating the N mineralization rate, and their relative importance varies with the changes in soil environment resulting from changes in management practice (de Vries et al., 2013; Kou et al., 2020; Wardle et al., 2004). In this study, there was an obvious spatial difference in the distribution of fungal and bacterial pathways in that the fungal pathway at 0–5 cm and the bacterial pathway at 5–15 cm were the driving factors in relating N mineralization to the soybean yield. This vertical distribution pattern of fungal and bacterial pathways is not surprising, as many studies (D'Hose et al., 2018; Sun et al., 2016; van Capelle et al., 2012) have reported that fungal and bacterial communities, which are the primary decomposers of fungal and bacterial pathways, also exhibit this same spatial pattern within the plough layer under the conservation tillage system. Residues under conservation tillage were placed on the soil surface instead of being mixed with the soil, resulting in large soil pores (as indicated by lower soil bulk density in 0–5 cm; Table S1) and a longer distance for soil microbes to gain access to nutrients in the upper soil. These environmental conditions are recognized as being more suitable for the growth of fungal communities (Moore et al., 2005), thereby promoting energy transfer through the fungi-based pathway in the near-surface layer.

Fungal and bacterial pathways differ in the N process rate as their components have different metabolic strategies (de Vries et al., 2013; Wardle et al., 2004). In contrast to the slower turnover rate of the fungal pathway, which favours N retention in the soil by immobilizing N in the biomass and organism-processed compounds, the bacterial pathway supports a faster N turnover rate, releasing more mineral N from biosynthesis into the soil solution (de Vries et al., 2013; Wardle et al., 2004; Whalen et al., 2013). In the present study, we found that the N mineralization from the bottom bacteria to the intermediate bacterial feeders and then to the top predaceous feeders was greatly enhanced in NT soils, suggesting a tight interlinkage and effective energy transfer across trophic levels in the bacterial pathway. These features of the food web have been recognized as playing a prominent role in promoting N turnover among immobilized and mobile forms (de Vries et al., 2013; Pressler et al., 2017; Wagg et al., 2014). Therefore, the enhanced N mineralization of the bacterial channel is expected to stimulate N mineralization and release more mineral N that can be readily absorbed by plants. This may partially explain why the severe shortage of soil mineral N empirically observed at 5–15 cm in NT soils (Table 1) during the growing season did not result in a compromise of soybean yield relative to CT.

The root pathway has been considered to have a very minor effect on N mineralization (Holtkamp et al., 2011; Pressler et al., 2017). In this study, the amount of mineralized N in the root pathway was indeed the least among different energy pathways across tillage systems (Fig. 2). However, to our surprise, the multiple linear regression analysis showed that there was a positive association between the mineralized N in the root pathway at 0–5 cm and soybean yield. This may be primarily due to the significant increase in mineralized N delivered by herbivores in the root pathway under RT and NT (Fig. 2 and Table S7). Verschoor (2002) reported that the N mineralization of herbivores accounted for 10 % of total N mineralization in a grassland system and attributed these beneficial effects of herbivores to the activity of soil microbes that were stimulated by the increase in root exudates after infection by herbivores. In our study, most groups classified into herbivores are the facultative feeders. For example, herbivorous collembolans can switch their diet from plant roots to decaying litter (Endlweber et al., 2009). Therefore, we propose that the positive role of herbivores at 0–5 cm in RT and NT soil may partly be due to their manipulation of surface residues by fragmenting and mixing. And then, the surface area of litter in contact with soil microbes would be increased, which is beneficial for N mineralization (Soong and Nielsen, 2016).

5 Conclusion

Combining the experimental data and the energetic food web modelling approach, our results suggest that, after long-term (14 years) application, conservation tillage has a larger potential mineralizable N pool as the soil food web becomes more complex relative to conventional tillage. Furthermore, soil organisms involved in associating mineralized N with soybean yield are different along the soil profile in which the fungal and root pathways at 0–5 cm and the bacterial pathway at 5–15 cm are the key driving factors for the supply of mineralized N to plants. Given that our finding is based on simulations and assumptions of steady-state soil biological communities resulting from a long duration of conservation tillage in the continental climate region, more studies, using isotope-tracing technique across different management practices, duration periods and climate regions, are needed to gain insight into how the soil food web processes energy and nutrients to maintain agroecosystem service and sustainability.