The responses of CO₂ emission to nitrogen application and earthworm addition in the soybean cropland

Study site

This research was performed at the farm of Jinming Campus of Henan University, Kaifeng City, Henan Province, China (34°49′N, 114°18′E). The crop was soybean and the variety was Kaidou 1104. A permanent 25 × 10 m² rainout shelter with steel frames and covered with a clear polyethylene roof was built in late summer 2021 to control precipitation inputs each year, to avoid the death of soybean and earthworms caused by extreme rainfall. The appropriate rainfall amounts were selected to simulate natural precipitation in the local area with the long-term rainfall means. The region belongs to the temperate continental monsoon climate with 14 °C of annual mean temperature and 650 mm (80% occurring between July and August) of annual average precipitation. The soil texture is sandy loam.

Experimental design

The experiment used a randomized block design involving two factors of nitrogen and earthworm, including six treatments: control (C), earthworm addition (E), high-frequency nitrogen application (NH), low-frequency nitrogen application (NL), earthworm and high-frequency nitrogen application (E*NH), and earthworm and low-frequency nitrogen application (E*NL). Each treatment was replicated five times with an area of 1 m × 1 m per plot. The total gram of earthworms (Metaphire guillelmi) was controlled at 8.0–8.9 g/m² (about 2–4 earthworms/m²) (Li et al., 2022). The endophytic earthworms (Metaphire guillelmi) used in the experiment were purchased from breeding company. Each plot was surrounded by glass to prevent earthworms from escaping. The glass fits tightly into the plot and has enough height (40 cm) to prevent earthworms from escaping. Nitrogen (urea) was added by water dissolving and root topdressing. The total amount of nitrogen added was the same as the conventional local field nitrogen application. The experiment included two frequencies of N application (two times vs. 12 times): High frequency nitrogen was applied once every seven days from July 15th to September 30th, and low frequency nitrogen was applied once every 30 days from July 15th to September 30th. From seedling stage (VE) to drum stage (R6), high-frequency nitrogen was uniformly added 12 times with 29 N kg·hm⁻² each time and low-frequency nitrogen was uniformly added twice with 174 N kg·hm⁻² each time. The first addition of earthworms was made at the soybean seedling stage in July 15th (Houida et al., 2022), and we continuously monitored the biomass of earthworms by electroshocking, and a second addition was made at the soybean podding stage in August 21st of 2022 to maintain a constant biomass (Wang, 2022).

Measurement of soil respiration

Soil respiration and temperature were measured every seven days during the soybean growing season using a Li-8100 portable soil CO₂ flux system (Li-Cor, Inc., Lincoln, NE, USA) and a thermocouple probe (Li-8100-201; Li-Cor, Inc., Lincoln, NE, USA) connected to the Li-8100 from June 2022 to October 2022. Soil volumetric water content at 0–10 cm soil depth was determined adjacent to each collar using a soil detector (TR-6D). All measurements were performed between 9 a.m. and 11:30 a.m. To avoid the respiration of aboveground parts of plants and litter decomposition, all living plants and litter inside the collars were removed by hand 2 days before soil respiration was measured. If it rains heavily, our measurement would be postponed for 2 days.

Soil sampling and analysis and plant index measurements

In September 2022. Three soil cores were randomly collected from each plot at a depth of 0–10 cm using a soil corer (inner diameter 5 cm) and mixed into one, sieved through a 2 mm mesh to separate gravel and roots, and divided into three parts. One subsample was stored at 4 °C for the analysis of the chemical properties of soil i.e., the available ammonium and nitrate using a colourimetric method (Smart Chem 200 Discrete Auto Analyser; Systea, Anagni, Italy). Another subsample was air-dried and ground for analysis of pH, total N (TN) and total carbon (TC). The soil pH was measured with a soil pH meter (TR-6D). The soil TN and TC concentrations were measured by a Vario ELIII Elementar (Elementar Analysensysteme GmbH, Langenselbold, Germany) elemental analyzer. The third part was stored at −20 °C for the analysis of the microbial community diversity composition spectrum (Wang et al., 2022).

After removing the whole plant from the soil, rinsed it slowly with running water to separate the above and below-ground parts of the plant from the cotyledon nodes. The washed roots were dried in an oven at 65 °C to a constant weight, and the root biomass (RB) was weighed. Aboveground biomass (AGB) was weighed after two weeks of natural air drying. Grain yield (GRY) was measured by removing the mature pods from the plants, placing them in paper bags, and leaving them in a ventilated place for drying to constant weight. The plant height (PLH) was measured by selecting three plants from each plot and measuring them with a tape measure. The number of pods per plant (NPP) and number of grains per plant (NGP) were also artificially measured by choosing three plants from each plot (Ji et al., 2017). A total of 100 grain weight (W100) was chosen at random from the grain yield of each plot and weighed using a precision scale (Ji et al., 2017).

DNA extraction, PCR amplification, and Illumina sequencing

Data were collected as previously described in Li et al. (2022). Specifically using E.Z.N.A Soil DNA Kit (Omega Bio-Tek, Norcross, GE, USA), soil DNA was extracted from each sample in accordance with the manufacturer’s protocol. The purity and concentration of the extracted DNA were determined by a NanoDrop-2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) (Li et al., 2022). Purified soil DNA was fully pooled together after quantitative determination and then for downstream manipulations (Wang et al., 2022). The V3–V4 of bacterial 16S rRNA genes were amplified with the following universal primer set: upstream primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and downstream primers 806R (5′-GGACTA CHVGGGTWTCTAAT3′). For fungi, the primers ITS5 (5′-GGAAGTAAA AGTCGTAACAAGG -3′) and ITS2 (5′-GCTGCG TTC TTCATCGATGC-3′) (Usyk et al., 2017) were used to amplify the ITS_V1 region of the rDNA gene. PCR reactions were performed in 25 μL reaction mixtures containing 5 μL 5 × reaction buffer, 5 μL 5 × GC buffer, 2 μL 2.5 mM dNTPs, 1 μL Forwardprimer (10 uM), 1 μL Reverseprimer (10 uM), 2 μL DNA Template, 8.75 μL ddH₂O, 0.25 μL Q5 DNA Polymerase. The reaction conditions were programmed of an initial denaturing step at 98 °C for 2 min, denaturation 98 °C 15 s, annealing 55 °C 30 s, extension 72 °C 30 s, final extension 72 °C 5 min and 10 °C hold 25–30 cycles. Samples were sequenced in an Illumina MiSeq High-Throughput Sequencing (HTS) platform (Illumina, San Diego, CA, USA) at Personal Biotechnology Co. Ltd Shanghai, China to determine soil microbial community composition.

Statistical analyses

Data were collected as previously described in Wang et al. (2022). Specifically, the sequenced data was performed using QIIME 2 2019.4 with slight modification. Raw sequence data were demultiplexed using the demux plugin followed by primers cutting with cutadapt plugin. Sequences were then merged, filtered and dereplicated using functions of fastq_mergepairs, fastq_filter, and derep_fulllength in Vsearch. All the unique sequences were then clustered at 98% (via cluster_size) followed by chimera removing. At last, the non-chimera sequences were re-clustered at 97% to generate OTU representative sequences and OTU table. Representative sequences were aligned with mafft and used to construct a phylogeny with fasttree. Alpha-diversity metrics (Observed_species, Simpson) were estimated using the diversity plugin with samples were rarefied. Meanwhile, principal coordinates analysis (PCoA) was selected to illustrate the clustering of different samples. In this study, Pco1 and Pco2 were used to represent the β diversity of microbial communities. OTUs was given a taxonomy using the Silva v132 99% OTU reference sequences and the classify-sklearn nave Bayes taxonomy classifier in the feature-classifier plugin (Liu et al., 2020b).

Two-way ANOVAs (p < 0.05) was used to analyze the significant differences between nitrogen application frequency and earthworms on CO₂ emission and soil properties. One-way ANOVA with Duncan testing (p < 0.05) was used to evaluate the significant differences in soil properties and CO₂ emission among the six treatments. Linear regression analysis was used to study the relationship between soil respiration, soil property and soil microbial community under six treatments. Spearman’s correlation analyses were performed to assess the relationships between soil properties, respiration, plant biomass and microbial community. Soil chemical properties data were analyzed with SPSS software (version 26; IBM, Chicago, IL, USA). We conducted a classification random forest analysis to identify the major statistically significant microbial predictors of the composition (relative abundance: number of sequences of major phyla/class/order level) of bacteria and fungi acting on soil respiration. The analysis was conducted using the rfPermute package of the R (4.2.2; R Core Team, 2023) statistical software. The significant predictors from random forest analysis were further selected for structural equation modeling (SEM) analysis. SEM analysis was applied to determine the direct and indirect contributions of soil properties and the bacterial community to CO₂ emission. SEM analysis was performed using AMOS 22.0 software (SPSS, Chicago, IL, USA). The model fitness was evaluated by χ2 (p > 0.05), comparative fit index, and root mean square error of approximation.

Nitrogen and earthworm application effects on soil respiration

Soil respiration varied with Soybean growth period, showing obvious seasonal variation (Fig. 1A). It was the lowest during the early vegetative stage (June 23rd–July 14th), reached maximum at the reproductive stage (July 21st–Sep. 7th), and started declining during the maturity period (Sep. 7th–Sep. 17th). Total soil respiration varied during the study from 1.13 to 4.63 µmol m⁻² s⁻¹, with an average of 2.55 ± 0.12 µmol m⁻² s⁻¹ (Fig. 1B). Compared with C, soil respiration in the NL increased significantly by 25%. E*NH, E*NL, NH, and E increased by 21%, 21%, 12%, and 11%, respectively. Soil respiration was significantly affected by nitrogen application, but not by earthworm addition or interaction between nitrogen and earthworm addition (Table 1).

Nitrogen and earthworm application effects on soil properties and plant index

Nitrogen application significantly influenced soil pH, TN, NO₃⁻-N, and NH₄⁺-N. No effects of earthworm on soil properties were found. A significant interactive effect between earthworms and nitrogen application was found on grain yield (Table 1). The E*NL treatment had the highest soil TN, which was also significantly higher than the C and E treatments (Fig. 2). In comparison to the C treatments (Fig. 2), the soil NO₃⁻-N was considerably greater in nitrogen application treatments (NH, NL, E*NH, and E*NL). Soil NH₄⁺-N was significantly higher in NH and E*NH than in other treatments (Fig. 2). Soil pH of E*NL, NL, NH, and E*NH decreased by 3.38%, 3.03%, 2.17% and 1.64% respectively and E increased 1.26% compared with the C treatments. The highest grain yield was observed in E*NH, which was also significantly greater than NH (Fig. 2). However, neither the application of nitrogen nor its interaction with earthworms had a significant influence on soil temperature, TC, plant biomass, plant height, NPP, NGP and W100.

Soil properties, plant index and microbial community regulated CO₂ emission

Soil respiration showed a positive correlation with Nitrate N (R² = 0.10, p < 0.05, Fig. S1), aboveground biomass (R² = 0.16, p < 0.05, Fig. S1), plant height (R² = 0.15, p < 0.05, Fig. S1), grain yield (R² = 0.12, p < 0.05, Fig. S1), and negatively with soil pH (R² = 0.18, p < 0.05, Fig. S1), but it was not correlated with root biomass (Fig. S1).

The random forest modeling showed the major bacterial and fungal phyla, classes, and orders for predicting soil respiration. These taxa include several bacterial and fungal such as Corynebacteriales, Myxococcales, Sordariomycetes, Verrucomicrobia, Thermoleophilia, Agaricomycetes, Gammaproteobacteria, Blastocladiomycetes, and Proteobacteria (Fig. 3). In addition, soil respiration was significantly positively correlated with the relative abundance of Gammaproteobacteria (R² = 0.149, p < 0.05, Fig. S2) and Corynebacteriales (R² = 0.154, p < 0.05, Fig. S2), as well as bacteria Pco1 (R² = 0.109, p < 0.05, Fig. S2), and negatively correlated with the relative abundance of Myxococcales (R² = 0.161, p < 0.05, Fig. S2), Verrucomicrobia (R² = 0.122, p < 0.05, Fig. S2) and Thermoleophilia (R² = 0.115, p < 0.05, Fig. S2). The Sordariomycetes, Agaricomycetes, Blastocladiomycetes and Proteobacteria were not correlated with soil respiration (p > 0.05). Notably, the relative abundance of Thermoleophilia, Myxococcales and Verrucomicrobia was low in NL, while the relative abundance of Corynebacteriales and Gammaproteobacteria was high. In general, NL treatment raised the relative abundance of copiotrophs while decreasing the relative abundance of oligotrophs.

The Spearman’s correlation coefficients between the microbial characteristics and the soil respiration as well as soil properties were estimated (Table 2). Soil respiration rate was negatively correlated to bacterial Simpson. In particular, there was a correlation between the respiration and the relative abundance of major bacterial and fungal phyla. The strong correlation between soil respiration and the abundances of Blastocladiomycetes and Agaricomycetes in the fungal compositions, as well as Gammaproteobacteria, Verrucomicrobia, Thermoleophilia, and Myxococcales in the bacterial compositions was found (Table 2).

The optimum structural equation model (SEM) revealed that the interactive networks of abiotic and biotic factors regulated soil respiration (χ² = 65.907, p = 0.515, CFI = 1.000, RMSEA < 0.001). Results of SEM explained 58% of the variations in soil respiration (Fig. 4). Changes in pH and NO₃⁻-N caused by nitrogen application altered bacterial richness and the relative abundance of Myxococcales, which affected soil respiration.

Effects of nitrogen application frequency on soil respiration and microbial community

Nitrogen application increases crop yield and quality (He et al., 2022). But fertilization practices can significantly impact on soil CO₂ emissions (Wang et al., 2021). In the study, fertilizer treatments considerably improved soil CO₂ emission by 12–25% in comparison to the control treatment, which is consistent with previous study findings (Yan et al., 2021; Lamptey et al., 2019). The frequency of nitrogen application (two and twelve times) had different pulse effects and influenced response variables such as soil nitrogen, bacterial diversity and root foraging requirements, which demonstrated the importance of nitrogen application frequency on nitrogen deposition in ecosystems (Cao et al., 2019). The NL resulted in lower pH (Ning et al., 2015; Ning et al., 2022) and NH₄⁺ (Zhou et al., 2016), but greater acidity than that the NH. The reason is that the H⁺ generated by nitrification reacted with the salt ions adsorbed on the colloidal surface of the soil, and the H⁺ was adsorbed on the surface of the soil, resulting in a decrease in soil pH and soil acidification (Sun et al., 2006). Moreover, the single nitrogen application of NL was higher than that of NH, which was difficult for plants to absorb in a short time. More nitrogen is used for microbial nitrification, so more NH₄⁺ is available for nitrification, and ultimately more H⁺ is produced (Zhou et al., 2016). Lower soil pH due to nitrogen application led to increases in soil respiration through changing the bacterial community (He et al., 2022; Whitaker et al., 2014; Liu et al., 2020a). In addition, high-frequency applications resulted in stronger levels of soil N content. Soil NO₃⁻-N is higher in the NH than the NL treatments, the reason is that higher NH₄⁺ exacerbates the toxicity of alkali ions inhibiting growth and NO₃⁻-N uptake of root (Hao et al., 2022; Ning et al., 2022). Our study revealed that NO₃⁻-N and soil pH were important environmental factors explaining changes in bacterial community composition, which is consistent with previous studies (Chen et al., 2021a).

Most bacteria are sensitive to pH variations and prefer being neutral (Anil, Menka & Preeti, 2019). Moreover, the bacterial communities control the major of ecological activities in agricultural ecosystems (Van der Heijen, 2008). In our study, acidification of alkaline soil made an ideal habitat for the soil microbial community. NO₃⁻-N is a substrate that provides nutrients to microorganisms and is associated with increased richness of soil bacterial communities, which is consistent with previous studies (Chen et al., 2021a). Consequently, the reduced pH and NH₄⁺ toxicity resulted in the highest bacterial richness in the NL, which directly enhanced soil respiration. In this study, NL treatments increased the bacterial richness, whereas decreased the fungi proportion. Higher bacterial richness was consistent with enhanced soil respiration, suggesting that increased metabolic rates drive this reaction (Hagerty et al., 2014).

Effects of earthworm on soil respiration and microbial community

Earthworm activity can increase the contents of soil active organic carbon, soil inorganic nitrogen, microbial biomass carbon, and microbial biomass nitrogen (Yu et al., 2007; Li et al., 2022). Moreover, earthworm addition may also increase soil CO₂ emissions (Lubbers et al., 2013). In this study, earthworm addition enhanced soil CO₂ emissions, which is consistent with previous research results (Lubbers et al., 2013; Song et al., 2020; Yang et al., 2019). Earthworm addition did not significantly increase soil respiration in our study, which is attributed to the fact that low density addition and type of earthworm were endophytic. Endophytic earthworms are mainly active in deeper soil layers, which have the characteristics of slow sexual maturity, limited population and weak respiration (Zhang et al., 2011).

The increase in pH and NH₄⁺-N caused by the addition of earthworms was consistent with previous studies (Wang et al., 2013; Ferlian et al., 2020; Na et al., 2023). Earthworm activity increases soil pH, which is related to the large amount of mucus secreted by earthworm epidermis and the amino acids, sugars and inorganic salts in earthworm gut excreta (Ferlian et al., 2020; Na et al., 2023). In the study, earthworm addition affected soil respiration by altering bacterial community composition. It has been shown that soluble carbon and other compounds released in the digestive tract of earthworms contribute to bacterial proliferation (Barbosa et al., 2017). Changes in pH and available nitrogen after the application of earthworm casts altered microbial community composition (Zhao et al., 2016), which is consistent with the changes in NO₃⁻-N caused by earthworm addition altering the relative abundance of Myxococcales and affecting soil respiration in our study.

Effects of nitrogen application frequency and earthworm on soil respiration

In the study, the addition of earthworms and the interaction between earthworms and nitrogen had no significant effect on soil respiration, which may be because the reduction of respiration rate caused by the absorption of sufficient nitrogen by plants and the increase of respiration rate caused by earthworm addition offset each other (Yang et al., 2019). Another possible explanation is that earthworm addition is not adapted to the habitat and cannot survive for a long time, therefore it has a short-term effect on soil respiration (Frouz et al., 2014). The interaction between earthworms and high-frequency N application showed an additive effect, whereas the interaction between earthworms and low-frequency N application reduced the effect of low-frequency N application on soil respiration. This suggests that the presence of earthworms buffers the pulse of low-frequency nitrogen application.

Microbial community regulated CO₂ emission

Important microbial classification and functional properties have been reported to potentially predict changes in soil respiration. Based on the study of the random forest model, we were able to pinpoint the primary microbial taxa that predict soil respiration. We also found that the dominant bacteria in the microbial community explained 35.82% of the variation in soil respiration (Fig. 3). Dominant bacteria of Myxococcales have been shown to be key factors in regulating soil respiration in our study, which is consistent with a previous study (Liu et al., 2021). Proteobacteria have been found to prefer soils with abundant carbon availability, hence this phylum might encourage increases in respiration. In addition, soil respiration was substantially connected with Gammaproteobacteria, Sordariomycetes, and Agaricomycetes (Ren et al., 2018), and Verrucomicrobia (Han & Wang, 2023), which is consistent with predicted species by random forest model in our study. The Sordariomycetes, Agaricomycetes, Blastocladiomycetes and Proteobacteria did not correlate with soil respiration in our study, which may be attributed to competition among microorganisms and the importance of dominant species (Ren et al., 2018). According to generally believed accounts, copiotrophs and oligotrophs have ability to use C for respiration. In general, it has been suggested that oligotrophs have slower respiration rates than copiotrophs (Liu et al., 2018; Chen et al., 2021b; Liu et al., 2020b). Myxococcales and Gammaproteobacteria are considered potential copiotrophs, but Verrucomicrobia, Thermoleophilia, and Corynebacteriales are classified as oligotrophs (Chen et al., 2021b). Moreover, the relative abundance of oligotrophs (Thermoleophilia and Corynebacteriales) is lowest in the low-frequency nitrogen application treatments, which can help explain the highest soil respiration in the NL plots. Importantly, our revealed that soil properties like pH and NO₃⁻-N affected soil respiration rates indirectly by changing the soil microbial community rather than directly (Liu et al., 2021; Wang et al., 2022). Our results highlight the importance of keystone taxa and microbial community composition as predictors in the development of the soil carbon model.