Manure application increased denitrifying gene abundance in a drip-irrigated cotton field

Site description and experimental design

A plot-based field experiment was conducted at the Cele Research Station (37°01′N, 80°43′E) of the Chinese Academy of Sciences in the 2015 and 2016 summer growing seasons. The region has a typically arid continental climate with an extremely low long-term average annual precipitation of only 42 mm, mainly distributed between May and July. The long-term average mean annual air temperature is 12.7 °C. The soil is classified as Aridisols in the USDA soil taxonomy system or Gypsisol in the FAO World Reference Base for Soil Resources. At the start of the study, the surface soil (0–20 cm) was a fine sand texture (sand 90%, silt 4%, clay 6%) with bulk density 1.46 Mg m⁻³, pH_H2O 8.0, electrical conductivity (EC) 144.4 μS cm⁻¹, total Kjeldahl N 0.31 g kg⁻¹, extractable NO₃⁻-N 25.7 mg kg⁻¹, 0.5M NaHCO₃-extractable P 14.6 mg kg⁻¹, 1.0M ammonium acetate K 153 mg kg⁻¹ and organic matter 6.9 g kg⁻¹. Analysis of soil characteristics was based on Carter (1993). Before this study, the experimental field was cropped to cotton for over 5 years and received both manure and urea applications each year, in accordance with local farmer practice.

The experimental design was previously described in Kuang et al. (2018) and only treatments under drip irrigation used in the current study. Briefly, the study was a randomized complete block design of four treatments with four replicate plots for a total of 16 plots. Each plot was 10 m long × 6 m wide and was separated from the other plots by a 1.1-m non-cropped buffer area. Treatments included (1) an unfertilized control (Control), and application of 240 kg of available N ha⁻¹ in the form of (2) granular urea (Urea, 46-0-0), (3) mixture of sheep and cattle manure compost (Manure), and (4) 50% urea with 50% manure by weight (U+M). Local producers commonly use such N application rates for high-yielding cotton fields. For the Urea treatment, 20% of added N was banded in the plant row before planting, and the rest was applied with irrigation water as a schedule of 5% at 9 weeks, and 15% of added N each at 11, 14, 15, 16 and 17 weeks after planting. The manure was all applied before planting by broadcast-incorporation at 10 cm depth. The manure had a moisture content of 25% and a dry weight-based total N, P, K content of 15.6, 2.0 and 16.8 g kg⁻¹, respectively. Analysis of manure was done on subsamples digested with a mixture of perchloric, sulfuric and hydrofluoric acid. Total P and K in the acid digestion were measured using the Mo-Sb colorimetric method and atomic absorption spectrometry (ICE3500; Thermo Fisher, Waltham, MA, USA), respectively. Total N was determined colorimetrically after Kjeldahl digestion. The manure had an available N concentration of 29 mg N kg⁻¹, determined by the alkaline hydrolyze method. In each year, cottonseed (c.v. Xinluzao 48; Huiyuan Tech, Shihezi, China) was planted at 75 kg ha⁻¹ in early to middle April under plastic-mulch and drip-irrigated, which is typical of recent cotton production in the region. Details on the system were described by Kuang et al. (2018). Before seeding, all plots received a broadcast-incorporated application of 120 kg P₂O₅ ha⁻¹ as calcium phosphate and 60 kg K₂O ha⁻¹ as K₂SO₄.

Soil sampling

Soil samples (0–20 cm depth) were collected with a hand auger (2.5 cm diameter) in September 2016 with cotton at the boll opening stage. For each plot, four soil cores were collected next to the drip tape and mixed thoroughly together for one composite sample per plot. The auger was cleaned using 95% alcohol and wiped with sterile paper before collecting the next soil sample. Each sample was passed through a two mm mesh screen and partitioned into three subsamples. One subsample was air-dried at room temperature for chemical analysis. The second subsample for analysis of denitrifying enzyme activity (DEA) and nitrifying enzyme activity (NEA) was stored at −20 °C and analyzed within 1 week. The third subsample for microbial molecular analysis was stored at −80 °C.

Soil chemical properties

Soil NH₄⁺ and NO₃⁻ was extracted using 0.01M CaCl₂ and measured with a continuous flow analyzer (SEAL Analytical, Norderstedt, Germany). Soil pH was measured at 1:2.5 soil:water ratio. Soil total C was measured using by wet oxidation method with potassium dichromate. Total N was analyzed by Kjeldahl acid-digestion method with a Kjeltec 1035 analyzer (Tecator AB; Sweden). Available Fe and Cu were extracted with diethylenetriamine pentaacetic acid (DPTA) (0.005M DPTA + 0.1M triethanolamine + 0.01M CaCl₂ set to pH 7.3) and analyzed using ICP-OES (730ES; VARIAN, Palo Alto, CA, USA). Soil dissolved organic carbon (DOC) was extracted using deionized water (1:5 soil:water ratio) and analyzed using a TOC analyzer (Aurora 1030W; OI Analytical, College Station, CA, USA). Soil C:N ratio was calculated on mass basis of total C and total N.

Determination of denitrifying and nitrifying enzyme activity

The frozen soil samples were pre-incubated to thaw at 25 °C for 2 days before analysis of DEA and NEA. Soil DEA was expressed as the rate of N₂O production (μg N h⁻¹ g⁻¹ soil) and determined using the anaerobic slurry technique (Beauchamp & Bergstrom, 1993). Briefly, 25 g thawed soil samples were placed into 125 ml plasma flasks. A total of 25 ml solution including 10 mM KNO₃, 10 mM glucose, 50 mM K₂HPO₄ and 0.1 g L⁻¹ chloramphenicol to inhibit new protein production was added to each flask. The flasks were evacuated and flushed with a 90:10 He-C₂H₂ gas mixture to create anaerobic conditions and suppress N₂O-reductase activity. Flasks were then shaken for 60 min and gas samples taken 0, 15, 30, 45 and 60 min after the onset of mixing using an orbital shake (180 rpm). Concentrations of N₂O in gas samples were immediately analyzed using gas chromatography equipped with an electron capture detector (Agilent 7890A; Agilent Technologies, Santa Clara, CA, USA).

Soil NEA was expressed as μg NO₃⁻-N h⁻¹ g⁻¹ dry soil and determined according to Hart et al. (1994). Briefly, a thawed soil sample (15 g dry soil equivalent) was placed into a 250 ml plasma flask with 100 ml solution of 1.5 mM (NH₄)₂SO₄ and one mM phosphate buffer (pH = 7.2). The flask was incubated at room temperature under constant agitation (180 rpm). Samples of the slurry were taken at 2, 4, 8, 12 and 24 h during incubation. Concentrations of NO₂⁻ and NO₃⁻ in the samples was then determined using the continuous flow analyzer. NEA rate was calculated from the slope of the regression model of NO₂⁻ plus NO₃⁻ concentrations and time.

Soil DNA extraction and real-time PCR

Soil DNA was extracted from 0.3 g of a soil sample using the Power Soil Total DNA Isolation Kit (MO-BIO Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer’s instructions. The quality and concentration of DNA were estimated using a Nanodrop 1000 Spectrophotometer (Thermo Fisher, Waltham, MA, USA) and gel electrophoresis (1.0% agarose). The DNA extracts were diluted at a ratio of 1:10 with double-distilled water (ddH₂O) to reduce the potential for PCR inhibition and then stored at −20 °C until use.

Quantitative PCR was used to quantify archaeal amoA and bacterial amoA, narG, nirK, nirS and nosZ gene in triplicate. All reactions were carried out using a CFX96™ (BIO-RAD, Laboratories Inc., Hercules, CA, USA). Each PCR reaction mixture contained one μl of 10-fold diluted soil DNA as a template, 10 μl SYBR^® Premix Ex Taq™ II (TaKaRa, Kusatsu, Japan), 0.8 μl of primer (10 μM), and 7.4 μl ddH₂O in a total volume of 20 μl. Primers and thermocycling conditions used in the qPCR reactions are given in Table 1. Plasmids containing respective sequences of the targeted genes were generated by cloning the targeted gene fragments from soil DNA into plasmid pMD™ 19-T Vector (TaKaRa, Kusatsu, Japan). Standard curves for each gene were created from 10-fold serial dilutions (10⁸–10¹) of the known quantities of linearized plasmid DNA harboring the target gene sequences. All qPCR reactions were conducted in triplicate. The qPCR efficiency and slope were 92% and −3.5 (R² = 0.990) for archaeal amoA, 105% and −3.2 (R² = 0.999) for bacterial amoA, 90% and −3.7 (R² = 0.999) for narG, 85% and −3.7 (R² = 0.997) for nirS, 96% and −3.4 (R² = 0.998) for nirK, and 80% and −3.5 (R² = 0.990) for nosZ, respectively. The generally low qPCR efficiency for nirS and nosZ genes agreed with previous studies reporting similar levels (74–90%, Ding et al., 2015; Harter et al., 2014).

High-throughput sequencing

The 16S rRNA gene of the V3-V4 hypervariable region was analyzed by MiSeq sequencing on the Illumina Miseq 2 × 300 bp platform at Shanghai Sangon Biotech Co., Ltd. with the universal primers 515F (GTGCCAGCMGCCGCGG) and 907R (CCGTCAATTCMTTTRAGTTT) that amplify both bacteria and archaea DNA (Li et al., 2014). Both forward and reverse primers were added with a barcode. The thermocycling program was set as: an initial denaturation at 94 °C for 3 min, five cycles at 94 °C for 30 s, 45 °C for 20 s, 65 °C for 30 s of extension, then 20 cycles of 94 °C for 20 s, 55 °C for 20 s, 72 °C for 30 s, with a final extension at 72 °C for 5 min. The reactions were set as: 15 μl 2 × Taq Master Mix (Thermo Scientific, Waltham, MA, USA), two μl of DNA template (about 20 ng), one μl of each appropriate primer (10 μM), 11 μl of ddH₂O. The PCR products were purified and quantified using the Agencourt AMPure XP (Beckman Coulter, Brea, CA, USA) reagent and Qubit™ ssDNA Assay Kit (Life Technologies, Carlsbad, CA, USA), respectively. Finally, the purified PCR products of each sample were equally combined based on their concentrations and produced a DNA pool which included 16S rRNA gene amplified fragments for sequencing.

Sequencing reads were allocated to each sample based on their unique barcodes. Raw sequences were firstly processed using cutadapt software (v 1.2.1) to trim the barcodes of their primers. Two short Illumina reads were then merged with PEAR (v 0.9.6) software (Zhang et al., 2014), and finally, PRINSEQ software (v 0.20.4, Schmieder & Edwards, 2011) was used for quality control of the merged reads. Only sequences >200 bp in length with an average quality score >40 were used for further analyses. Chimeras were filtered by comparing the sequences with those in the reference database using the UCHIME algorithm (v 4.2.40, Edgar et al., 2011). After the above screening, the remaining high-quality sequences were clustered into operational taxonomic units (OTUs) at a ≥97% similarity identity threshold. The singletons and low abundance OTUs were removed before further analyses. The Ribosomal Database Project classifier (Wang et al., 2007) was used to identify taxonomic information at the bootstrap cutoff of 80%. Based on the OTUs output, α-diversity, and β-diversity, and canonical correspondence analysis (CCA) were performed. Species richness and diversity indices, including coverage, Chao1, ACE, Simpson and Shannon were calculated using mothur (v 1.30.1, Schloss et al., 2009) to estimate the α-diversity of each sample.

Statistical analysis

Treatment effects on soil properties, NEA, DEA, α-diversity indices and bacterial gene abundances were conducted using a one-way ANOVA. Pearson correlation analysis was conducted to assess the relationships between the functional gene abundances, NEA, DEA and selected soil properties. ANOVA and Pearson correlation analysis was performed with SAS 9.3 (SAS Institute, Cary, NC, USA) and differences were considered significant at P < 0.05. Principal coordinates analysis (PCoA) was performed to determine the community β-diversity of bacteria using the Vegan package Version 1.17-7 (Oksanen, 2011) implemented with the R language, which was based on a bacterial weighted UniFrac metric matrix. CCA was performed with the Vegan package implemented with the R language to determine the relationships between soil physicochemical properties and bacterial communities. Untransformed data were used for the PCoA and CCA analyses. The relative abundances of bacterial community at the phylum level between treatments were compared using the Welch’s t-test with STAMP (Statistical Analysis of Metagenomic Profiles). Corrected P-values of the Welch’s t-test were calculated using the false discovery rate for multiple testing correction.

Soil chemical characteristics

Manure and U+M treatments increased soil total N content by half compared to the unfertilized control (Table 2). Soil NO₃⁻ concentrations with Manure and U+M treatments were 120 and 103 mg kg⁻¹, respectively, being 2.4–4.8 times greater than the Urea and Control treatments. In contrast, soil NH₄⁺ concentrations were not affected by the treatments. Soil total C and DOC were also greater in Manure and U+M compared to the Control and Urea treatments. As a result, treatments with manure addition (Manure and U+M) had 37–100% higher soil C:N ratios, compared to Urea and Control.

Denitrifying enzyme activity and nitrifying enzyme activity

Manure and U+M treatment significantly (P < 0.001) increased DEA levels compared to the Urea and Control treatments (Fig. 1A). In contrast, NEA levels did not respond significantly to any amendment treatment in spite of an increasing trend from control to urea and manure additions (Fig. 1B).

Bacterial community, and nitrifier and denitrifier gene abundances

The sequence coverage index ranged between 0.93 and 0.94, suggesting that the sequencing depth was sufficient to obtain the majority of genetic diversity of samples (Table 3). The average number of valid sequences were similar for treatments, being 27,131, 28,812, 26,413 and 30,405, for Control, Urea, Manure and U+M, respectively, with a mean read length of 376 bp. The number of OTUs was not affected by N addition and ranged between 4,072 and 4,295. The indexes of richness and diversity, Chao1, ACE, Simpson and Shannon were also not affected by the treatments.

The addition treatments resulted in an apparent clustering in β-diversity of the soil bacterial community (Fig. 2). Two groups with (Manure and U+M) and without (Control and Urea) manure application occurred along axis PCoA1, with a significant dissimilarity (P < 0.001). The PCoA explained 73% of the total variation in the composition of the bacterial community, with PCoA1 and PCoA2 explaining 61% and 12%, respectively. At the bacterial phylum level, the abundance of Planctomycetes, Bacteroidetes and Ignavibacteriae increased with manure additions, whereas that of Latescibacteria, Acidobacteria, Armatimonadetes, Actinobacteria and candidate division WPS-2 decreased (Table 4). There was no treatment effect on the abundance of Proteobacteria. At archaeal phylum level, the abundance of Thaumarchaeota and Euryarchaeota decreased with manure application.

Manure and U+M treatments doubled the gene copy number of AOA (Fig. 3A). Similarly, manure application also significantly (P < 0.001) increased the copy number of the AOB gene in manure than non-manure amended treatments. The copy number of AOA was generally one order of magnitude greater than that of AOB. Further, AOB copy number responded more to treatment additions than that of AOA. As a result, manure addition reduced the ratio of AOA/AOB, being 88.7, 27.6, 15.8 and 17.0 for Control, Urea, Manure and U+M, respectively.

Manure addition significantly (P < 0.001) increased the copy number of narG, nirK or nosZ genes, but did not affect that of nirS (Figs. 3B–3E). Copy number of narG was 27.5–39.0 times greater with manure (U+M and Manure) than non-manure (Control and Urea) addition treatments. Copy number of nirK was 3.4–3.7 times greater with manure than non-manure addition treatments. Similarly, the copy number of nosZ gene were 9.6-25.2 times greater in manure than non-manure addition treatments.

Relationships between DEA, NEA, soil properties and microbial abundance

Copy number of nitrifier (AOB and AOA), nitrate reducer (narG), nirK, and nosZ-type denitrifier genes, but not nirS-type denitrifier gene, were positively correlated with NO₃⁻, DOC, total N, total C and C:N ratio (Table 5). In contrast, soil NH₄⁺ and pH were not significantly correlated with the copy number of any of the functional genes. There were also significant positive correlations between DEA and abundance of AOB, AOA, narG, nirK and nosZ genes. NEA, however, was not correlated with the abundance of any functional gene, except for a positive correlation with nirK.

The bacterial β-diversity was highly associated with changes in soil environmental variables (Fig. 4). The first and second axes explained the variance in bacterial community structure to 38.6% and 13.9%, respectively. The bacterial community of manure addition treatments (Manure and U+M) was mainly associated with soil concentrations of NO₃⁻, DOC, total C, total N and C:N ratio. In contrast, the bacterial community of treatments without manure addition (Control and Urea) was mainly associated with soil pH.

Impact of N addition strategy on denitrification and nitrification

The increased activity of soil denitrifying enzymes with manure in the current study is inconsistent with our findings at the same field where we reported more N₂O emissions from manure compared with urea application under drip irrigation conditions (Kuang et al., 2018). The increased NO₃⁻ and carbon supply with manure application could likely have provided the primary substrate for denitrification and increased the N₂O/(N₂O + N₂) ratio (Francis et al., 2013). Chantigny et al. (2010) also suggested that manure can elevate soil respiration and deplete O₂ concentration to create temporary anaerobic conditions, thereby further increasing the proportion of N₂O production through denitrification. These studies highlight the importance of N addition source on soil N transformation processes and suggest that manure induced N₂O emissions are likely attributed to denitrification.

In contrast to DEA, NEA was not affected by manure application in the current study. Similarly, Shen et al. (2008) also reported that organic manure did not affect potential nitrification rates of an alkaline sandy loam soil in northern China. Several studies suggested that soil pH is the dominant factor for nitrification as it determines the availability of NH₄⁺, which is the primary substrate for ammonia oxidation, the initial and rate-limiting step of nitrification (Fan et al., 2011; Nicol et al., 2008). In our study, both pH and the availability of NH₄⁺ were not affected by N addition strategy, confirming the insensitivity of NEA to N sources.

In contrast with manure, urea did not significantly affect DEA and NEA compared to Control. Our results agree with those of Yin et al. (2015), who reported that manure but not inorganic fertilizer increased denitrification potential. In contrast, application of inorganic N fertilizers increased the activity of nitrification (Fang et al., 2018; Shi et al., 2016) and potential denitrification (Duan et al., 2017; Wang et al., 2018). The absence of the inorganic fertilizer effect in the current study was associated with the minor to no effect by urea application on soil properties such as pH, DOC and inorganic N (NO₃⁻ and NH₄⁺) compared with Control. It is likely the buildup of C and N substrates by urea application were not sufficient enough to affect the activities of functional genes.

Impact of N addition strategy on the abundance of functional genes and bacterial community structure

In the current study, the positive relationships of the abundances of narG, nirK and nosZ with DEA and further with soil DOC, total C and total N suggest that manure significantly increased gene abundance by providing C and N substrate. This result is in line with the previous findings that the denitrifiers abundance could be used as a predictor of DEA (Morales, Cosart & Holben, 2010). Our findings also agree with previous studies which reported that organic manure increased the abundance of nosZ-type denitrifier compared to inorganic fertilizers (Hallin et al., 2009; Tao et al., 2018). Also being consistent with previous studies (Zhou et al., 2011), the abundance of nirK but not nirS was increased by manure application in this study, suggesting that nirK was more susceptive to fertilizer regimes than nirS-type denitrifier. Hallin et al. (2009) also reported that denitrification rates were not correlated with the abundance of nirS genes in soils treated with different fertilizer regimes for 50 years. A possible reason for the lack of correlation could be that the denitrifier harboring the nirS gene might play a minor functional role for DEA (Attard et al., 2011). In the current study, the nitrate reducer (narG), and nirK- and nosZ-type denitrifiers, which encodes the main catalytic enzymes responsible for nitrate reduction, nitrite reduction, and N₂O reduction respectively, were more sensitive to manure application. The increase of the denitrifiers abundance with manure application thus increased the pool of denitrifying enzymes. Even though the limited soil sampling for microbial analysis hindered the possibility of directly linking results from the current study to the in-situ measurements of N₂O flux, the positive relationship between DEA and the abundance of denitrifiers suggest the manure-induced N₂O emissions in Kuang et al. (2018) was more likely determined by denitrification.

It is interesting to note that manure application increased AOA and AOB whereas did not affect NEA in this study, suggesting the abundance of ammonia-oxidizers are not necessarily associated with nitrification potential. Nicol et al. (2008) reported that the activity of ammonia-oxidizers was more associated with the relationships among transcription, translation, and enzyme function rather than the abundance of functional genes. It is also likely that the complicated subsequent hierarchical regulation of enzyme expression resulted in an uncouple effect between NEA and amoA gene abundance (Röling, 2007). Consistent with previous studies (Fan et al., 2011; Tian et al., 2014), N addition reduced the AOA/AOB ratio in this study, suggesting that AOA and AOB may occupy different soil niches due to the differences in physiological and metabolic pathways. Previously, AOA prefers low NH₃ substrate conditions for growth, whereas AOB prefers higher NH₃ levels (Di et al., 2010), thus potentially resulting in a lower AOA/AOB ratio following N addition.

Similar to previous studies (Ji et al., 2018; Kumar et al., 2018; Wang et al., 2019), manure application significantly changed β-diversity of soil bacterial community in the current study. The PCoA analysis revealed a dominant contribution of PCoA1 (61%) to total variation and a clear separation of manure vs. non-manure groups along the axis PCoA1. This suggests that the addition of manure was a key factor determining the variation in the bacterial community among treatment. Clearly, the increased N and C substrates with manure application have increased the growth of some specific microbial groups and suppress others and thus changed the composition of the soil microbial community. The absence of urea effect on the β-diversity of the bacterial community was attributed to the low organic matter content (6.6 g kg⁻¹), suggesting that the substrate deficiency of C limited microbial activities under the conditions in this study. In this study, manure or urea applications did not influence α-diversity of soil bacterial community, likely due to an absence effect on soil pH. Fierer & Jackson (2006) reported that soil pH is the primary driver determining the α-diversity and richness of the soil bacterial community.

In the current study, the changes of soil bacterial community structure in response to manure application were attributed to the increasing relative abundance of Planctomycetes, Bacteroidetes, lgnavibacteriae and decreasing abundance of Actinobacteria, Acidobacteria, Latescibacteria, Armatimonadetes and candidate division WPS-2. These results highlight the change of eutrophic and oligotrophic bacteria. For example, Fierer, Bradford & Jackson (2007) found that Bacteroidetes were typically copiotrophic bacteria and could thrive in soil with high available organic carbon. Planctomycetes are involved in the turnover of soil organic carbon, and nutrient availability, and the reproduction of this microbial group may increase intensively in response to the application of manure (Lupatini et al., 2016). The phyla which were negatively influenced by manure application were considered as slow-growing oligotrophs accustomed to nutrient-limited environments. For example, several studies had shown that Acidobacteria strains grew slowly with their growth being limited to substrate additions (Goldfarb et al., 2011).