Ammonia-oxidizing bacterial communities are affected by nitrogen fertilization and grass species in native C₄ grassland soils

Experimental design and sample collection

This study was conducted at a small-plot trial at the University of Tennessee East Tennessee Research and Education Center (ETREC) in Knoxville, TN (35.53°N, 83.06°W). Soils at this site are classified as sandy loam (fine-loamy, mixed, semiactive, thermic Typic Hapludults) with 59.9% sand, 23.5% silt, and 16.6% clay (Bandopadhyay et al., 2020). The experimental design was a randomized complete block with split-plot treatment arrangements. The main plots had two C₄ native grass species (switchgrass (SG) (Panicum virgatum) and big bluestem (BB) (Andropogon gerardii)) and the split-plots had different urea fertilization rates (0, 67, and 202 kg N ha⁻¹ (0N, 67N, and 202N, respectively)). Each combination of grass species and N rate had three replications. Grasses were planted in 2013, treatments were applied in 2014. Full experimental details, including the application of other fertilizers, lime, and herbicides were described in detail in our previous study (Hu et al., 2021b).

After five years of the experimental treatments, soil samples were collected as previously described (Hu et al., 2021b). Briefly, samples were collected at three grass growing seasons in 2019: grass green up (G, late April, one week before the first N fertilization), initial grass harvest (H1, late June, within one week after harvest and before the second N fertilization), and second grass harvest (H2, Mid-August, within one week before harvest). Six soil cores (2.5 cm diameter, 10 cm long) were collected from each sub-plot and composited. A 10 g sub-sample was immediately stored at −80 °C for DNA and RNA extraction, and the remainder of the soil sample was used for soil physicochemical properties and nitrification potential measurements.

Soil physicochemical analysis and nitrification potential

Soil pH, soil water content (SWC), NH₄⁺-N, NO₃⁻-N, and nitrification potential were analyzed using freshly collected soil. Total organic carbon (TOC), total nitrogen (TN), dissolved organic carbon (DOC), and dissolved organic nitrogen (DON) were analyzed using air-dried soils using established methods described in our previous study (Hu et al., 2021b). Soil nitrification potential was measured using a chlorate block method (Belser & Mays, 1980). Briefly, soil samples (2 g) were mixed with 14 mL nitrification potential medium (4 mL 0.2M K₂HPO₄, 0.5 mL 0.2M KH₂PO₄, 2.5 mL 0.2M (NH₄)₂SO₄ and 10 mL 1M NaClO₃ in 983 mL deionized water) in 50 ml tubes and shaken for 4-6 h (175 rpm, 25 °C). The nitrification potential medium is a solution containing adequate NH₄⁺ and PO₄³⁻ to maximize nitrification rates along with chlorate to block nitrite oxidation to nitrate. one mL slurry subsamples were removed at 3 h intervals and centrifuged to remove soil particles. NO₂⁻ concentrations were determined by adding 50 µL of NO₂⁻ color reagent (4 g of sulfanilamide and 0.1 g N-(1-naphthyl)-ethylenediamine hydrochloride and 17 mL 85% phosphoric acid in 87.6 mL deionized water) to 200 µL of centrifuged samples in a 96-well plate. After 15 min, the absorbance was measured in a microplate reader spectrophotometrically at λ = 543 nm. We used absorbance values for the NO₂⁻ standards to compute the NO₂⁻ concentration for each sample, subtracting the non-soil controls. Nitrification potential was calculated as µg NO₂⁻-N g⁻¹ dry weight soil hour⁻¹ (gdw⁻¹ h⁻¹).

Soil N₂O emissions

Soil N₂O emissions were measured using the static chamber method (Collier et al., 2014) on the same days that soil samples were taken. A circular stainless-steel base (15 cm in diameter, 16 cm in height) was inserted into the soil to a depth of approximately eight cm at the center of each sub-plot. The chamber was closed with a lid and gases were sampled through a valve in the lid with plastic syringes (30 mL) at four-time intervals (0, 20, 40, and 60 min). The samples were transferred to pre-evacuated glass vials (12 mL). Soil temperature was measured by inserting a digital thermometer (Thermo Scientific, USA) one inch into the soil adjacent to the chamber.

Gas samples were analyzed in a gas chromatograph (Shimadzu GC-2014, Japan) with an electron capture detector for N₂O determination. N₂O flux was calculated by linear interpolation of the four sampling times. The N₂O-N emission flux was calculated as follows: ${\displaystyle F=\rho \times h\times \mathrm{d}c/\mathrm{d}t\times 273/\left(273+T\right)\times k}$ where F is emission flux (µg N m⁻² h⁻¹); ρ is the N₂O-N density in standard state (1.25 kg m⁻³); h is the net height of static chamber (0.08 m); dc/dt is the rate of gas concentration change (µL L⁻¹ min⁻¹); T is the average soil temperature in the static chamber during gas sampling (°C); k is the time conversion factor (60 min h⁻¹).

Soil nucleic acid extraction and cDNA synthesis

Nucleic acid extraction and cDNA synthesis methods were described in our previous study (Hu et al., 2021b). Briefly, soil genomic DNA and RNA were extracted using the DNeasy PowerSoil Kit and RNeasy PowerSoil Total RNA Kit (Qiagen, Hilden, Germany), respectively, according to the manufacturer’s protocols. DNA and RNA concentrations were quantified using the NanoDrop One spectrophotometry (NanoDrop Technologies, Wilmington, DE). RNA was checked for DNA contamination by attempting to PCR amplify amoA genes; samples that produced a negative electrophoresis result (i.e., no DNA contamination) were used for reverse transcription. cDNA was produced with RNA as the template using SuperScript IV Reverse Transcriptase and random hexamer primers (Invitrogen, Paisley, UK) according to the manufacturer’s protocol. RNA was stored at −80 °C; DNA and cDNA samples were stored at −20 °C.

Quantitative analysis of amoA genes and gene transcripts

Quantitative PCR (qPCR) of AOB amoA genes and transcripts were performed on a CFX96 Optical Real-Time Detection System (Bio-Rad, Laboratories Inc., Hercules, CA, USA) using the primer pair amoA-1F (5′- GGGGTTTCTACTGGTGGT -3′) and amoA-2R (5′- CCCCTCKGSAAAGCCTTCTTC -3′) (Rotthauwe, Witzel & Liesack, 1997) as we have previously described (Hu et al., 2021a).

amoA gene amplicon sequencing

A two-step PCR was used for amoA amplicon sequencing library preparation. We used methods similar to our previous study on nifH genes (Hu et al., 2021b) for the amoA genes. First, amoA genes were amplified from diluted DNA extracts (10 ng µl⁻¹) with primers amoA-1F/amoA-2R containing Illumina-compatible adapters (Klindworth et al., 2013). The amplifications were performed in a 25-µl mixture containing 12.5 µl KAPA HiFi HotStart ReadyMix (2 ×) (Kapa Biosystems), 1 µl of each primer (final concentration of 0.4 µM), 2.5 µl of template DNA, and 8 µl of PCR grade nuclease-free water. The PCR protocol was: 95 °C for 3 min, then 35 cycles of 98 °C for 20 s, 58 °C for 15 s and 72 °C for 15 s, followed by 72 °C for 5 min. PCR products were purified using SparQ PureMag Beads (Quantabio). For the second step, index PCR was performed in a 50-µL reaction containing 25 µl KAPA HiFi HotStart ReadyMix (2 ×) (Kapa Biosystems), 5 µl of each forward/ reverse NexteraxT index/primer, 5 µl of PCR product from the first step as template DNA, and 10 µl of PCR grade nuclease-free water. The PCR protocol was: 95 °C for 3 min, then 8 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s, followed by 72 °C for 5 min. PCR products were purified using SparQ PureMag Beads, quantified by NanoDrop, pooled, and run on an Agilent Bioanalyzer to ensure quality (Hu et al., 2021b). 2  × 275-bp paired-ends sequencing was run on an Illumina Miseq platform (Illumina, CA, U.S.) in the University of Tennessee Genomics Core laboratory.

Bioinformatics analysis

amoA sequences were processed using Mothur v.1.39.5 (Schloss et al., 2009). Forward and reverse reads were merged using make.contigs. Primer sequences were trimmed using trim.seqs. Reads were removed using screen.seqs if the overlap region of the merged reads was less than 10 bp, if the reads were shorter than 400 bp, or if the reads contained any ambiguous bases. Similar to the pipeline used in our previous study (Hu et al., 2021b), chimera.uchime was used to remove chimeras based on a manually created database containing 9682 AOB amoA gene reference sequences FunGene (http://fungene.cme.msu.edu/) (Fish et al., 2013). For all downstream analyses, all libraries were subsampled to a uniform sequencing depth (4,307 reads). Operational taxonomic units (OTUs) were classified based on 97% nucleotide sequence similarity (Wang et al., 2014; Zheng et al., 2014) using vsearch in QIIME2-2019.7 (Bolyen et al., 2019). Singletons (OTUs appearing only once across all samples) were removed, and representative sequences for each OTU were selected. Abundant OTU representative sequences (total reads ≥ 100) were blasted against the FunGene AOB amoA-like database using BLASTn. The raw sequencing data can be accessed from NCBI Sequence Read Archive (SRA) Database with BioProject accession number PRJNA733230.

Statistical analysis

The statistical analyses used in this research were similar to our previous study (Hu et al., 2021b). Diversity analyses of AOB community were performed with QIIME2-2019.7. Chao1, observed OTUs, Shannon index, and Pielou’s evenness were calculated to determine alpha-diversity (Gourmelon et al., 2016; He et al., 2013). A weighted-UniFrac method (Lozupone et al., 2007) was used for analyzing beta-diversity, similar to our previous study (Hu et al., 2021b). Permutational multivariate analysis of variance (PERMANOVA) was used to examine beta-diversity by treatment factors (Hu et al., 2021b). A neighbor-joining phylogenetic tree was constructed by QIIME using representative sequences of 28 abundant OTUs (OTUs with total reads ≥ 100) and was visualized using Interactive Tree of Life (iTOL, v5) (Letunic & Bork, 2011) as described in our previous study (Hu et al., 2021b). Principal coordinates analysis (PCoA) and redundancy analyses (RDA) based on weighted-UniFrac distances was performed in R 3.6.1 (The R Foundation for Statistical Computing) with packages vegan (2.5–5) and phyloseq.

A mixed model ANOVA within the GLIMMIX procedure in SAS 9.4 (SAS Inst., Cary, NC, USA) was performed to test treatment effects on soil properties, abundance and expression of amoA genes, and alpha-diversity of AOB communities, similar to our previous approach (Hu et al., 2021b). Gene and transcript abundances were log transformed to achieve normal distributions. Fixed effects included sampling season, N fertilization rate, and grass species, as well as their interactions. Random effects included the block and block by grass species. A least significant difference (LSD) method was used to compare the means of each two groups. Pearson correlation analyses were performed in IBM SPSS Statistics v26 (Hu et al., 2021b). IBM SPSS Amos 27.0 software was used to perform structural equation modeling (SEM) to characterize interactions among measured variables. We followed the procedures of developing and modifying a structural equation model in Li et al. (2019) and Hu et al. (2021b). Because community is a theoretical concept and not directly measurable, the “AOB community” component of the SEM model was a latent variable inferred from the measured relative abundances of the top 28 most abundant amoA OTU. These were OTUs with ≥ 100 reads, which comprised 97.9% of the total reads and were were classified into four species: Nitrosospira sp. 56-18, Nitrosospira sp. APG3, Nitrosospira multiformis, and Nitrosospira sp. Nsp14. For grass species, we defined switchgrass = 1, big bluestem = 2 due to the higher biomass of big bluestem than switchgrass. A well fit model should have a Chi-square/Df value in the range of 1.0 to 3.0 with probability level > 0.05, RMSEA value < 0.05, and fit indices (TLI, IFI, and CFI) both > 0.95 (Blunch, 2012; Hooper, Coughlan & Mullen, 2008; Kline, 2015).

Soil physicochemical characteristics, N₂O emission, and nitrification potential

The soil physicochemical characteristics under different treatment combinations were described in detail in our previous study (Hu et al., 2021b). Briefly, SWC, DOC, and DON showed significant differences by sampling season (p < 0.001 for SWC and DOC; p = 0.006 for DON) (Table 1; Table S1). SWC was highest (0.30 ± 0.01 g H₂O g⁻¹ dry weight soil (gdw⁻¹)) at initial grass harvest (H1) and lowest (0.15 ± 0.01 g H₂O gdw⁻¹) at second grass harvest (H2). DOC and DON were highest (223.00 ± 5.47 µg gdw⁻¹ and 22.44 ± 1.98 µg gdw⁻¹, respectively) at H2 and lowest (163.85 ± 4.07 µg gdw⁻¹ and 16.16 ± 0.65 µg gdw⁻¹, respectively) at H1 (Table S1) (Hu et al., 2021b).

Nitrogen fertilization rate and C₄ grass species affected several soil properties (Table 1). With increased N fertilization rate, NO₃⁻-N increased from 0.39 ± 0.16 to 2.77 ± 0.39 µg gdw⁻¹, soil pH decreased from 6.36 ± 0.05 to 6.12 ± 0.05, and DOC decreased from 210.18 ± 8.60 to 186.07 ± 5.23 µg gdw⁻¹. Soil pH was higher under switchgrass (6.32 ± 0.05) compared to big bluestem (6.19 ± 0.03). Fertilized plots had higher SWC and lower C:N for switchgrasss than 0N plots. In contrast, fertilization of big bluestem resulted in lower SWC and no change in C:N (Table S1).

N₂O emission rate was affected by the interaction of sampling season, N fertilization rate, and grass species (Table 1; Table S1): at H1, N₂O-N increased with N-rate, from 10.01 ± 2.00 to 23.77 ± 2.08 g N₂O-N ha⁻¹d⁻¹ in switchgrass plots and from 8.74 ± 3.92 to 74.98 ± 29.33 g N₂O-N ha⁻¹d⁻¹ in big bluestem plots; but at G and H2, no effect of N fertilization rate and grass species on N₂O-N was observed.

There was a significant interaction effect between sampling season and N fertilization rate on nitrification potential (Table 1; Table S1): at G, nitrification potential was not significantly different under different N-rates, but at H1 and H2, nitrification potential was significantly higher under 202N. Nitrification potential was significantly higher in switchgrass plots (0.14 ± 0.02 µg NO₂-N gdw⁻¹h⁻¹) than in big bluestem plots (0.10 ± 0.02 µg NO₂-N gdw⁻¹h⁻¹) (Fig. 1; Table S1).

Abundances of amoA genes and transcripts

amoA gene abundances significantly varied with N fertilization rate, and grass species (Table 1). In general, amoA gene abundances were lowest at H1 (2.51  × 10⁵ copies g⁻¹ dry weight soil (gdw⁻¹)) and highest at H2 (1.14 × 10⁷ copies gdw⁻¹) (Fig. 2A). N fertilization increased amoA gene abundance (5.03  × 10⁶ and 7.00  × 10⁶ copies gdw⁻¹ for 67N and 202N, respectively) compared to no N addition (1.70  × 10⁶ copies gdw⁻¹) (Fig. 2B). In addition, amoA genes were more abundant in switchgrass plots (5.34  × 10⁶ copies gdw⁻¹) than in big bluestem plots (3.81  × 10⁶ copies gdw⁻¹) (p = 0.0066) (Fig. 2C).

amoA gene transcript abundances were significantly affected by the interaction effect of sampling season and N fertilization rate (Table 1): at G, amoA transcript abundance was not different between N-rates; at H1, amoA transcript abundance was highest under 67N; at H2, amoA transcript abundances increased with increasing fertilizer. Grass species did not affect amoA transcript abundance (Table 1 and Fig. 2C).

Ammonia-oxidizing bacterial structure and community

Using a 97% sequence similarity cutoff, a total of 215,943 high-quality bacterial amoA gene sequences were clustered into 538 OTUs, with 28 to 134 OTUs per sample. N fertilization rate and grass species significantly affected the alpha-diversity of AOB community (Table S2). AOB species richness (observed OTUs, Chao1) was lower at H1, but alpha-diversity estimates were similar across sampling times (Fig. 3A). Observed OTUs, Shannon index, and Pielou’s evenness increased with N fertilization rate (Fig. 3B). Higher richness was observed in switchgrass plots than in big bluestem plots (Fig. 3C).

AOB community structure was significantly affected by N fertilization rate and grass species (Table 2; Fig. 4). N fertilization rate had stronger effect (PERMANOVA R² = 0.504, p = 0.001) than grass species (R² = 0.042, p = 0.021) (Table 2). PCoA showed that AOB community composition in 202N was different from 0N and 67N treatments (PERMANOVA p < 0.001). Under 0N and 67N, the communities were much more variable under big bluestem compared to switchgrass (Fig. 4).

Phylogenetic analysis of ammonia-oxidizing bacterial community

There were 28 abundant amoA OTUs (OTUs with ≥ 100 reads) which comprised 97.9% of the total reads. A neighbor-joining phylogenetic tree showed that most of the 28 abundant amoA OTUs classified as Nitrosospira, with 8 OTUs (accounting for 11.8% of total reads) classified as Nitrosospira multiformis (Fig. 5). Many amoA OTUs were treatment-specific; for example, 8 OTUs (OTU4, OTU8, OTU9, OTU10, OTU12, OTU16, OTU19, and OTU22) classified as Nitrosospira multiformis and 2 OTUs (OTU11, OTU23) classified as Nitrosospira sp. 56-18 were more abundant in 202N compared to 67N or 0N. OTU19 was only detected in big bluestem plots with 202N, whereas OTU3 and OTU6, belonging to Nitrosospira sp. Nsp14, were more abundant under 0N and 67N in big bluestem plots. OTU2 (classified as Nitrosospira sp. Nsp14) had higher relative abundance under 0N and 67N in both switchgrass and big bluestem plots (Fig. 5).

Relationships between soil properties and AOB populations

The abundances of amoA genes and transcripts were not correlated to each other (Table S3). The amoA gene abundance was negatively correlated to soil pH (R = −0.303, p < 0.05), SWC (R = −0.792, p < 0.01), NH₄⁺-N (R = −0.471, p < 0.01) but positively correlated to DOC (R = 0.477, p < 0.01), DON (R = 0.445, p < 0.01), and nitrification potential (R = 0.530, p < 0.01). In contrast, the amoA gene transcript abundance was positively correlated to SWC (R = 0.561, p < 0.01), and available inorganic N (NH₄⁺-N: R = 0.448, p < 0.01; NO₃⁻-N: R = 0.339, p < 0.05), and negatively correlated to DOC (R = −0.640, p < 0.01) (Table S3). Alpha-diversity indices (except Chao1 index) were negatively correlated to soil pH but positively correlated to NO₃⁻-N and nitrification potential (Table S3). The Chao1 index, observed OTUs, and Shannon index were positively correlated to amoA gene abundance but not transcript abundance (Table S3).

Structural equation modeling (SEM) was conducted to identify the impacts of soil physicochemical parameters on AOB populations (Fig. 6; Table S4). SWC had a direct negative effect on amoA gene abundance (standardized path coefficient = −0.711, p < 0.001) and a direct positive effect on amoA transcript abundance (0.433, p = 0.010). amoA transcript abundance was directly and positively affected by NH₄⁺-N concentration (0.331, p = 0.002) but directly and negatively affected by DOC (−0.493, p < 0.001). In addition, amoA gene abundance had a direct positive effect on its transcript abundance (0.541, p < 0.001). The Chao1 index was directly and negatively affected by NO₃⁻-N (−0.366, p < 0.001), whereas the Pielou’s evenness of AOB community was directly and positively affected by NH₄⁺-N (0.268, p < 0.001). N₂O emission rates were directly and positively affected by SWC (0.723, p < 0.001) and NH₄⁺-N (0.167, p = 0.022). The total effects of N fertilization rate on amoA gene and transcript abundances, alpha-diversity of AOB community, nitrification potential and N₂O emission rate were positive (Table S5). For the most abundant species within the AOB community, the total effects of N fertilization rate were positive for Nitrosospira sp. 56-18, Nitrosospira multiformis, and Nitrosospira sp. APG3 but negative for Nitrosospira sp. Nsp14. Compared to big bluestem, switchgrass had higher amoA gene and transcript abundances, alpha-diversity of AOB community, nitrification potential, and soil NO₃⁻-N content but grass species had no effect on AOB community composition (composed of the four Nitrosospira species mentioned above).

A redundancy analysis (RDA) showed that soil physicochemical properties, amoA gene and transcript abundances, sampling season, N fertilization rate, grass species, and alpha-diversity indices measured in this study explained a total of 62.5% of the variability of AOB community structures, with the first two axes explaining 50.2% of this variability (Fig. 7). Soil pH, NO₃⁻-N, N₂O emission, nitrification potential, amoA transcript abundance, N fertilization rate, grass species, and alpha-diversity significantly correlated to the variability of AOB community structures (p < 0.05) (Fig. 7).