Effects of different proportions of organic substitution for mineral fertilizers on soil methanogenic and methanotrophic communities in paddy fields

Site description

The long-term experiment, initiated in 1984, was conducted at the Institute of Soil Fertilizer and Resource Environment, Jiangxi Academy of Agricultural Sciences, Nanchang, Jiangxi Province (28°33′92″N, 115°56′25″E; 25 m). The average annual temperature and rainfall were 17.5 °C and 1,600 mm, respectively. The soil at the test site was a medium-retention yellow mud field developed from a Quaternary sub-red clay matrix, classified as Chromic Alisols (Cutanic, Dystric, and Loamic) according to the World Reference Base for Soil Resources (WRB) in 2022 (IUSS Working Group WRB, 2022). The cropping system consisted of double-season rice (early rice–late rice). The basic chemical properties of the soil from 0 to 20 cm depth before the experiment in 1984 were as follows: pH 6.5, organic matter 25.6 g/kg, total nitrogen 1.36 g/kg, total phosphorus 0.49 g/kg, alkaline hydrolyzable nitrogen 81.6 mg/kg, available phosphorus 20.8 mg/kg, and available potassium 35.0 mg/kg.

Experimental design

The experiment was designed as random blocks, with each plot covering 33.3 m² (11.1 × 3 m²) and an interval of 0.50 m between plots. The plots were separated by a 0.50-m-wide cement ridge, and each plot was independently drained and irrigated. The experiment included five treatments, each with nine replicates: no fertilizer (CK), 100% mineral nitrogen, phosphorus, and potassium (NPK) fertilizers, 70% mineral fertilizer and 30% organic fertilizer (M30), 50% mineral fertilizer and 50% organic fertilizer (M50), and 30% mineral fertilizer and 70% organic fertilizer (M70). For early rice, the 100% mineral fertilizer was applied with pure N at 150 kg/ha, P₂O₅ at 60 kg/ha, and K₂O at 150 kg/ha. For late rice, it was applied with pure N at 180 kg/ha, P₂O₅ at 60 kg/ha, and K₂O at 150 kg/ha (Table S1). The organic fertilizer for early rice was Chinese milk vetch (A. sinicus L.), containing N at 3.03 g/kg, P₂O₅ at 0.8 g/kg, and K₂O at 2.3 g/kg. The organic fertilizer for late rice was rotted pig manure, containing N at 4.5 g/kg, P₂O₅ at 1.9 g/kg, and K₂O at 6.0 g/kg. The treatments (except CK) were designed to have equal nutrient quantities based on nitrogen quantity, with mineral fertilizers filling any phosphorus and potassium shortfalls. Phosphorus and organic fertilizers were applied as basal fertilizers. Notably, 50% of nitrogen fertilizer was used as basal fertilizers (Early Rice: April 20; Late Rice: July 18), 25% in the tillering stage (Early: May 26; Late rice: August 20), and 25% in the spike stage (Early Rice: June 10; Late Rice: September 9). Further, 50% of potassium was used in the tillering stage and 50% in the spike stage.

Sample collection and analysis

Early rice was transplanted in late April and harvested in mid-July, whereas late rice was transplanted in late July and harvested in late October. Both early and late rice were collected in the maturity stage, and their yields were measured. Soil samples were collected on October 15, 2021. A 3-cm auger was used to collect soil from 0 to 20 cm depth at four corners and the center of each plot, and these samples were mixed to create one composite sample per plot. The collected soil was processed to remove debris such as stones and plants and subsequently transported to the laboratory within 12 h. Upon arrival, fresh soil samples were sieved through a 2-mm stainless steel sieve and divided into three parts: one for chemical property determination (see supporting information), another for the microcosm experiment, and the last part for microbial sequencing.

Microcosm experiment

A 36-day microcosm experiment was conducted at Jiangxi Academy of Agricultural Sciences to investigate the impact of different treatments on methane emissions in paddy fields. The 200 g of freshly sieved soil was added to a 500-mL culture flask for each treatment, and then 200 mL of distilled water was added to simulate the flooded conditions of a rice field. The flask was then tightly sealed to maintain a closed environment. Nine replicates were set up for each treatment, and the incubator was maintained at a constant temperature of 25 °C throughout the experiment. A 25-mL gas sample was collected from each flask daily using a syringe, and methane levels were measured immediately using a gas chromatograph (Agilent 7890B; Agilent Technologies, La Jolla, CA, USA). The supplementary materials provide detailed information for this section.

DNA extraction

Genomic DNA extraction from 0.5 g of fresh soil was performed using the FastDNA SPIN Kit for soil (MP Biomedicals, Irvine, CA, USA) following the manufacturer’s protocols.

Quantitative real-time polymerase chain reaction

The primers MLf/MLr and A189F/mb661R were used to quantify the mcrA gene and pmoA gene, respectively (Bourne, McDonald & Murrell, 2001; Luton et al., 2002). The primer sequence for mcrA was GGTGGTGTMGGATTCACACARTAYGCWACAGC/TTCATTGCRTAGTTWGGRTAGTT, and the primer sequence for pmoA was GGNGACTGGGACTTCTGG/CCGGMGCAACGTCYTTACC. The real-time polymerase chain reaction (PCR) mixture consisted of 16.5 µL of ChamQ SYBR Colour qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China), 0.8 µL of the forward primer (5 µM), 0.8 µL of the reverse primer (5 µM), and 2 µL of template DNA (10 ng). The qPCR amplification program for the mcrA and pmoA genes included an initial hold at 95 °C for 5 min, followed by 40 cycles of denaturation at 95 °C for 5 s, annealing at 55 °C for 30 s, and extension at 72 °C for 40 s.

Illumina MiSeq sequencing

The aforementioned two primer sets were used to amplify mcrA and pmoA genes for sequencing methanogenic and methanotrophic communities. The PCR reaction system comprised 25 µL, including 4 µL of 5 × FastPfu buffer, 0.4 µL of TransStart FastPfu Polymerase (Transgen Biotech, Beijing, China), 2 µL of dNTPs (2.5 mM), 0.8 µL of forward primer (5 µM), 0.8 µL of reverse primer (5 µM), and 1 µL of template DNA. Then, the amplicons were extracted using a 2% agarose gel, purified following the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) instructions, and quantified using the QuantiFluor-ST system. The purified amplicons were used for paired-end sequencing (2 × 250) on the Illumina MiSeq platform. The reaction conditions for Illumina MiSeq sequencing of methanogens were initial denaturation at 94 °C for 3 min; 35 cycles of denaturation at 94 °C for 25 s, annealing at 55 °C for 45 s, and extension at 72 °C for 60 s, followed by a final extension at 72 °C for 5 min. The reaction conditions for Illumina MiSeq sequencing of methanotrophs were as follows: initial denaturation at 95 °C for 5 min; 35 cycles of denaturation at 92 °C for 1 min, annealing at 55 °C for 1.5 min, and extension at 72 °C for 60 s, followed by a final extension at 72 °C for 5 min.

Bioinformatics analysis

The data analysis was performed using QIIME 1 (version 1.9) software. Sequences with poor quality (average quality score below 30) were excluded from the analysis. The optimized sequences were clustered into Operational Taxonomic Units (OTUs) at a 93% similarity level using Usearch for both methanogen and methanotroph datasets (Lüke & Frenzel, 2011; Zhou et al., 2020a). The chimera and singleton OTUs were eliminated from the dataset. The representative sequence for each OTU was identified using the FunGene database (Fish et al., 2013). Randomly selected subsets comprising 52,767 and 30,462 sequences (the lowest sequence read depth) per sample were used to compare the compositions and diversity of methanogenic and methanotrophic communities, respectively, to mitigate potential calculation errors arising from varying sequencing depths across samples.

Statistical analysis

The one-way analysis of variance module in SPSS 20 software was used to analyze differences in soil chemical properties, total methane emissions, mcrA and pmoA abundance, and α diversity index across various treatments (Tukey HSD, P < 0.05). The Pearson correlation analysis assessed the correlation between methane emissions, microbial α diversity index, and soil chemical properties. Non-metric multidimensional scale analysis was performed to elucidate differences in methane-associated microbial communities across various treatments.

Rice yield and soil chemical properties

The long-term fertilization significantly increased both grain and straw yields compared with CK (Fig. 1). The organic substitution for mineral fertilizers led to higher grain yield than mineral fertilizer alone (Fig. 1A). In contrast, little difference in straw yield was observed among treatments except for CK (Fig. 1B). Fertilization resulted in significant changes in soil chemical properties (Table S2). Compared with CK (pH 5.60), mineral fertilization (NPK, pH 5.06) and low proportions of organic substitution (M30, pH 5.35) decreased soil pH levels. Conversely, a high proportion of organic substitution (M70, pH 5.93) caused a significant increase in soil pH. Compared with CK, mineral fertilization and organic substitution treatments increased soil fertility and soil organic (SOC) content except for total potassium (TK). In addition, organic substitution treatments, especially a high proportion of organic addition, appeared to be superior to mineral fertilization alone in increasing soil fertility and SOC content (Table S2).

CH₄ emission dynamics observed through microcosm experiments

The dynamic change of CH₄ emission flux during the microcosm experiment showed that CH₄ emissions were not prominent in the early stages of cultivation, with the peak occurring on the 18th and 23rd days (Fig. 2A). Significant differences were observed in the cumulative CH₄ emissions following flooding culture across various treatments, with cumulative emissions ranking as M50 > M30 > M70 > NPK > CK (Fig. 2B). Cumulative methane emissions peaked in the M50 group, reaching 10.21 µg/kg. Compared with mineral fertilization, treatment with a medium proportion of organic substitution significantly increased cumulative CH₄ emissions, whereas treatments with low and high proportions of organic substitutions had minimal effects (Fig. 2B).

Abundance of methanogens and methanotrophs

Compared with the CK group, both mineral fertilization and organic substitution treatments significantly increased the gene copy numbers (i.e., abundance) of mcrA and pmoA genes. For methanogens, the gene copy number of mcrA was the highest in the M70 group (5,520 copies/g), followed by the M50 group (4,123 copies/g) and the M30 group (2,965 copies/g). For methanotrophs, treatments with a high proportion of organic substitution exhibited an even higher abundance than those with low and medium proportions. The gene copy number of pmoA was the highest in the M70 group (26,672 copies/g), followed by the M50 group (6,418 copies/g) and the M30 group (3,119 copies/g) (Fig. 3). The Pearson correlation analysis revealed a significant positive correlation between the abundance of mcrA and the cumulative methane emission flux. Furthermore, the abundances of mcrA and pmoA showed significant positive correlations with soil pH, SOC, total nitrogen (TN), total potassium (TP), available nitrogen (AN), available phosphorus (AP), and available potassium (AK), whereas they displayed negative correlations with TK (Table S3).

Alpha diversity analysis of methanogens and methanotrophs

A total of 3,699,882 mcrA sequences were obtained across all samples, ranging from 52,767 to 119,523 sequences per sample. Additionally, 2,631,826 pmoA sequences were obtained, with per-sample counts ranging from 30,462 to 83,877 sequences. Further, 3,912 and 1,799 methanogenic and methanotrophic OTUs were detected from 45 samples, respectively. The organic substitution treatments (M30, M50, and M70) increased the OTU richness of methanogens compared with CK (Fig. 4A). Both mineral fertilization and organic substitution treatments increased the Shannon index of methanogens but had little effect on their evenness compared with CK (Figs. 4B and 4C). Mineral fertilization and organic substitution treatments showed little effect on the OTU richness of methanotrophs compared with CK (Fig. 4D). Mineral fertilization (NPK) decreased the Shannon index and evenness of methanotrophs. In contrast, treatment with organic substitution maintained similar levels of Shannon index and evenness of methanotrophs as CK (Figs. 4E and 4F). The OTU richness of methanogens was significantly positively correlated with SOC, TN, TP, AN, and AP, and negatively correlated with TK. The methanogenic evenness was significantly negatively correlated with SOC, TN, and AN. For methanotrophs, the OTU richness was significantly positively correlated with TK and negatively correlated with AN and AK. The Shannon index of methanotrophs was significantly positively correlated with pH and negatively correlated with AK. The methanotrophic evenness was significantly positively correlated with pH, SOC, TN, TP, and AP (Table S4).

Community composition of methanogens and methanotrophs

Fertilization significantly altered soil methanogenic and methanotrophic community compositions compared with CK. Treatment with organic substitution, especially at high proportions, exhibited a significant effect on these communities than mineral fertilization alone (Fig. 5 and Table S5). Soil pH and SOC were primary factors affecting both methanogenic and methanotrophic community compositions (Fig. 5). For methanogens, fertilization decreased the relative abundance of Methanocella paludicola and increased the relative abundance of Methanobacterium oryzae and Methanolinea mesophile. Treatment with low and medium proportions of organic substitution decreased the relative abundance of Methanospirillum psychrodurum and Methanoperedens nitroreducens. Treatment with organic substitution increased the relative abundance of Methanobacterium aggregans and Methanobacterium bryantii, whereas it decreased the relative abundance of Methanosarcina thermophila (Fig. S1).

The relative abundance of Type I and Type II methanotrophs tended to increase and decrease along organic substitution gradients, respectively (Fig. S2). The relative abundance of Methylomonas methanica increased and Methylocystis bryophila decreased under NPK fertilization and organic substitution treatments compared with CK. NPK fertilization resulted in the highest relative abundance of Methylocystis echinoides and Methylocystis parvus, and the lowest relative abundance of NC10, Methylomarinovum caldicuralii, and Methylococcus capsulatus compared with other treatments. The relative abundance of Methylogaea oryzae increased with a higher proportion of organic fertilizer substitution treatment (Fig. S3).