In March 2018, a composite soil sample (0–20cm) was collected from an experimental site at the Universidad Nacional del Sur Campus (38°41.64′ S, 62°14.46′ W) in Bahía Blanca (Argentina), to fill 5.5dm³ pots, each with 5.6kg of soil. The sample belongs to a well-drained Petrocalcic Paleustoll with 15years of A. sativa L. cropping history. Aldana and Carioni (2001) reported the following properties for this soil: sandy loam texture and bulk density 1.28gcm⁻³, while a soil analysis in 2018 indicated a pH_(H2O) 7.6 (1:2.5 soil:water), organic matter 2.46%, extractable (Bray) P 7.1mgkg⁻¹, N-NH₄⁺ 13.3mgkg⁻¹, N-NO₃⁻ 9.1mgkg⁻¹ and of total N 0.129%.

The greenhouse assay was conducted with A. sativa L. var. Cristal INTA. The experimental factors were fertilization level (with and without inorganic N fertilizer) and CC termination method (M: mowing, G: glyphosate, and U: untreated control) arranged as a 2 × 3 factorial in a completely randomized design, resulting in a total of six treatment–combinations, with four replicates (pots) per treatment. For glyphosate termination, the herbicide was applied as a commercial formulation (3lha⁻¹, Eskoba Full II, Red Surcos, 662gl⁻¹, monopotassium salt) using a knapsack sprayer. The total growth period from CC planting (April 3) to termination date was 92days (Z3.1 stage) with automatized irrigation twice a day during 5min through sprinklers. The N fertilizer (urea, 46% N) was surface applied at planting and tillering stage (day 76) at a total rate of 100kgNha⁻¹. Untreated control pots were destructively sampled at the time of termination to collect rhizospheric soil of living roots. For the experimental units terminated with glyphosate and mowing, the soil was sampled 12days later, to collect rhizospheric soil of plants with notable desiccation symptoms. In all cases, loosely adhering soil was removed by gentle shaking of the root system to discard bulk soil, and the tightly adhering soil (rhizospheric soil) was removed by brushing the root system with sterile brushes, collected on sterile trays, and stored in sterile plastic bags (Yanai et al., 2003). The soil was stored at −80°C for DNA extraction and molecular analysis and at 4°C for assessment of PNA.

Nitrification rate was measured according to the potential nitrification assay described by Hart et al. (1994), with modifications in the slurry preparation (10-fold scaling, i.e., 1.5g of fresh soil in 10ml of reaction buffer) to adapt volumes and soil quantities to a microplate assay, similar to Hoffmann et al. (2007). PNA reaction buffer (0.3mM KH₂PO₄, 0.7mM K₂HPO₄, 0.05mM (NH₄)₂SO₄, and 10mM KClO₃) was freshly prepared on the same day in which soil slurries were incubated. The slurries were prepared in 125ml sterile flasks loosely capped with aluminum foil to reduce water evaporation while leaving enough space for aeration. Flasks were incubated in an orbital shaker (25°C, 180rpm in the dark) for 1h to achieve a homogeneous suspension, and, afterward, 1ml samples were taken at 0, 2, 4, 20, and 22h (Hart et al., 1994). Samples were centrifuged at 15,000rpm (RCF: 21,379×g, 4°C) in a Hermle Z32-HK microcentrifuge (Wehingen, Germany) and the supernatant was transferred to sterile 1.5ml tubes. The tubes were immediately frozen at −20°C.

Quantitative analysis of nitrite concentration was conducted in 96-wells, polypropylene, flat-bottom microplates (Eppendorf AG, Hamburg, Germany, reference number: 0030 602.102) according to the diazotization method (Griess reaction) with sulfanilic acid (0.8% w/v in acetic acid; reference number: B1550261, Laboratorios Britania S.A., Buenos Aires, Argentina) and N-(1-Naphthyl)ethylenediamine 2HCl 0.1% w/v in acetic acid (Carbosynth^TM, Compton, United Kingdom). Absorbance was measured at 540nm and 25°C in a FLUOstar Optima microplate reader (BMG Labtech, Offenburg, Germany). The rate of production of nitrite was calculated by linear regression of solution concentration over time. After correction for soil moisture content of the soil and conversion to a soil dry weight (dw) basis, the potential nitrification activity (μgN-NO₂⁻ g⁻¹ dw soil h⁻¹) was calculated following Drury et al. (2006).

The collection of root exudates was done from mowed or glyphosate treated plants (unfertilized) with a solution of CaCl₂, as described by Egle et al. (2003). Briefly, all mineral particles attached to the root system in each pot were removed with tap water and the roots were submerged in 100ml of collection solution (CaCl₂ 0.05mM, pH 5.5) for 1h. The liquid was discarded and the roots were submerged again in 50ml of the collection solution for 16h. The resulting exudate obtained in each case was centrifuged 10min at 13,000×g, filter-sterilized, lyophilized, and stored at −20°C. The lyophilized exudate from each pot was suspended in 10ml of PNA reaction buffer and immediately used to assess the effect on ammonia oxidation. Two pools of root exudates (one from glyphosate-terminated plants and one from mowed plants) were obtained by mixing equal quantities of exudates from each pot of the same treatment group.

To assess the effect of root exudates on nitrifying communities, soil samples (approximately 1.5g) from each of the four replicates (pots) collected at time zero (before termination) were divided into three portions (each about 0.4g) to prepare soil slurries with either root exudates collected from glyphosate-treated plants, root exudates collected from mowed plants, or soil slurries with the PNA reaction buffer (control). The reaction buffer and the conditions of incubation were as described previously except for the sampling times (0, 2, 20, 22, and 24h). Soil slurries were prepared by combining the soil samples with 9ml of reaction buffer and 1ml of root exudate, or soil samples with 10ml of PNA reaction buffer for the control flasks.

DNA was extracted from 250mg of rhizospheric soil samples with the commercial kit PowerSoil DNA Isolation kit (Qiagen, Hilden, Germany) according to manufacturer instructions. DNA was quantified using QuantiFluor dsDNA kit in a Quantus fluorometer (Promega, Promega, Madison, WI, United States). The quality of the DNA was assessed from 1% agarose gel and absorbance ratios (260:230 and 260:280nm ratio) in a DS-11 FX spectrophotometer (DeNovix Inc., Wilmington, DE, United States).

To study the treatment effects on the estimated abundance of AOB and AOA, a quantitative real-time PCR (qPCR) of amoA gene of each group was conducted (amoA_AOB and amoA_AOA, respectively). Primers used for AOA and AOB as well as their amplicon lengths are indicated in Supplementary Table S1. To avoid introducing errors (mainly related to an unknown number of operons per cell), we analyzed the copy numbers (Ouyang et al., 2016; Morrison et al., 2017); no attempt was made to convert these copies into cell numbers. The amplification program for AOB and AOA as well as the PCR master mix volumes and reaction setup were previously reported by Zabaloy et al. (2016) and Zabaloy et al. (2017), respectively. All amplifications, baseline corrections, melting curve analysis, and standard curve assessments were conducted in ABI 7500 Real−Time System and its associated software (7500 Software v2.0.3, Applied Biosystems, Foster City, CA). The equation of the standard curve for amoA_AOB (Ct=39.084–3.881 log₁₀ copy number) represented an assay efficiency of 81% with an R² value of 0.993, while the equation for amoA_AOA (Ct=33.55–3.76 log₁₀ copy number) represented an efficiency of 84% with an R² value of 0.996.

Three replicates of each treatment were analyzed by amplicon sequencing of amoA gene. Gene libraries were prepared with the primers amoA1F/amoA2R (amoA_AOB) and CrenamoA23f/CrenamoA616r (amoA_AOA; Supplementary Table S1), using the Fluidigm™ protocol at the DNA Service Laboratory, Roy J. Carver Biotechnology Center at University of Illinois (Urbana-Champaign, United States). The library was quantitated by qPCR and sequenced on an Illumina MiSeq platform (Illumina, San Diego, CA, United States) using one MiSeq flowcell for 251cycles, from each end of the fragments. FASTQ files were generated and demultiplexed by the sequencing service with the bcl2fastq v2.20 Conversion Software (Illumina, San Diego, CA, United States).

The FASTQ files were processed in QIIME2 (Bolyen et al., 2019) using the recommended pipelines for paired-end reads (AOB) or single-end reads (AOA) based on DADA2 denoising algorithm. Briefly, primers were removed using p-trim-left argument within dada2 denoised-paired script (AOB) or within dada2 denoised-single script (AOA). For AOA, we used dada2 denoised-single script (i.e., only forward reads) considering that the length of the reads (300bp) was not enough for merging forward and reverse reads of the amplicon and, thus, a pipeline for paired reads would be inappropriate in this case. The trimming argument within dada2 denoised-paired script was set on 273bp for forward reads and 240 for reverse reads (median Q values >27 reported by Interactive Quality Plot tool of QIIME2), while a default value of 2 was set for the expected number of errors. For dada2 denoised-single, script we used the same expected number of errors, while no trimming was applied as the quality values were acceptable (median Q values >33). Single-ends analysis yields sequences that are inherently shorter than in paired-end analysis, thus a trimming step would shorten the length of the amplicon even further for taxonomic analysis.

The resulting amplicon sequence variants (ASVs) were tabulated and the ASV table was rarefied to a value equal to the sum of reads of the sample with the lower number of reads. The rarefied data of ASV table was used in vegan package 2.5–7 (Oksanen et al., 2020) of R Statistical Software v4.1.0 (R Core Team, 2021) for calculation of the following alpha diversity metrics: Shannon diversity index (H′), observed richness (S′=number of ASVs), and Pielou’s evenness index (J’). The ASVs obtained after all processing steps (removal of primers, trimming and filtering, denoising, and chimera removal) were clustered into operative taxonomic units (OTUs) using qiime vsearch cluster-features-de-novo at the appropriate species-level identity threshold for amoA gene (90%; Ouyang et al., 2016), similarly to other studies (Ouyang et al., 2016; Lourenço et al., 2018). For calculation of UniFrac distances, the representative sequences of the resulting OTUs were then aligned using Muscle algorithm in MEGAX (Kumar et al., 2018) and the aligned sequences were loaded in phangorn v2.7.0 package (Schliep, 2011) of R Statistical Software v.4.1.0 to obtain the distance matrix and a neighbor-joining tree. The likelihood of this tree was computed (pml function) and then optimized (optim.pml) using the most appropriate model of nucleotide evolution according to the Bayesian Information Criterion (BIC; Schliep, 2018). The phylogenetic tree and the OTU table were used in GUniFrac package v1.2 (Chen, 2021) of R Statistical Software v.4.1.0 to calculate generalized UniFrac distances.

Phylogenetic trees containing the OTU sequences and database sequences of amoA gene of AOB and AOA were constructed using maximum likelihood method with 1,000 bootstraps in MEGAX (Kumar et al., 2018). Tamura 3-parameter model and Hasegawa-Kishino-Yano model were used for AOA and AOB, respectively, according to the lowest Bayesian Information Criterion of model selection. We used high-quality amoA sequences from FunGene database (Fish et al., 2013) also included by Lourenço et al. (2018), with a score above 350, a size greater than 200 amino acids in length, a hidden Markov model (HMM) coverage of more than 85%, and a defined organism. Additional sequences from members of Nitrosospira clusters 0, 2, and 4, also included in reported AOB phylogenetic trees of rhizospheric soil (Glaser et al., 2010), were included as well. For AOA, we used several amoA sequences from marine environment, soil, sediments, and the rhizosphere of Zea mays L., all previously used in the construction of phylogenetic trees of AOA (Park et al., 2006; He et al., 2007; Tourna et al., 2011; Yao et al., 2011; Ai et al., 2013). Interactive Tree of Life (iTOL) tool was used for tree visualization (Letunic and Bork, 2016).

All statistical analyses were conducted in R Statistical Software v.4.1.0 (R Core Team, 2021). Quantitative PCR data (copy number μg⁻¹ DNA) and PNA values were analyzed using two-way ANOVA and least-square means procedure for post hoc multiple comparisons of means (α=0.05) after a log₁₀ transformation in emmeans package v.1.6.1 (Lenth, 2021). The interaction term was considered significant at p<0.2 (Littell et al., 2002). The alpha-diversity metrics were analyzed using two-way ANOVA (α=0.05) with the rarefied ASV table. The PNA in response to the addition of root exudates was analyzed using a linear mixed model in nlme package v3.1–152 (Pinheiro et al., 2021) considering the pot (replicate) from which soil was originally extracted as a random effect, and the treatment (glyphosate root exudate, mowing root exudate, or control buffer) as a fixed effect. Correlation analyses between the copy number of amoA genes and the PNA were performed using Pearson method (α=0.05) with cor.test function in base R.

Amplicon sequencing data was analyzed using both the standard approach as well as the compositional approach (Gloor and Reid, 2016; Gloor et al., 2017) using several packages of R Statistical Software v.4.1.0. In the first case, the UniFrac distances were used as input in vegan package v2.5–7 (Oksanen et al., 2020) for multivariate analysis of beta diversity through non-metric multidimensional scaling (NMDS) using metaMDS function. In the second case, centered log-ratio (clr) transformed data (Aitchison, 1986) was first calculated after replacing zero values with zCompositions package v.1.3.4 (Palarea-Albaladejo and Martin-Fernandez, 2015) and then used as input for principal component analysis (PCA) with PCA function (singular value decomposition) using FactoMineR package v.2.4 (Lê et al., 2008). The covariance biplots of compositional data were obtained using the packages FactoMineR, factoextra v.1.0.7 (Kassambara and Mundt, 2020), and ggplot2 v.3.3.3 (Wickham, 2016). The principal components (PCs) with eigenvalues larger than the average of all the eigenvalues (Kaiser-Guttman stopping rule) were retained (Jackson, 1993). The coordinates of the individuals in these components were used as input in a two-way ANOVA to assess the significance of fertilization and termination method on each component, with Tukey’s test (α=0.05) for post hoc multiple comparison of means with the emmeans package (Lenth, 2021).

In a complementary multivariate analysis, the generalized UniFrac distance and the Aitchison distance (i.e., Euclidean distance of clr-transformed data) were used as input in PERMANOVA (α=0.05; Anderson, 2001) with the adonis function (1,000 permutations) of vegan package v2.5–7 (Oksanen et al., 2020) to study the effects of the termination method, fertilization, and their interaction on community structure. When necessary, pairwise PERMANOVA were conducted with FDR p-value adjustment method in RVAideMemoire package v 0.9–80 (Hervé, 2021). The homogeneity of multivariate dispersions among treatments was visualized in vegan package v2.5–7 using betadisper function and a permutation test (α=0.05; 999 permutations).

A statistically significant main effect of the termination method (p<0.001) was detected for amoA_AOA. Mowing resulted in a higher amoA_AOA abundance in the rhizosphere compared to control plants (1.49-fold change, p<0.01), and to those terminated with glyphosate (1.76-fold change, p< 0.001). No significant differences were observed for amoA_AOA between glyphosate termination and the untreated control (p>0.05; Table 1, Supplementary Figure S1). For the N fertilizer, a marginally significant main effect was observed (p= 0.06) with a higher copy number of amoA_AOA detected in the rhizosphere of unfertilized plants relative to the fertilized plants (1.14-fold change, Table 1).

A differential response of amoA_AOB to the addition of N fertilizer was observed in the rhizosphere of glyphosate-treated plants relative to the response observed in mowed or in untreated plants (p= 0.053 for the interaction term). While no N-induced stimulation of AOB was observed in glyphosate-treated plants, a significantly greater amoA copy number was observed after N fertilization in the rhizosphere of mowed and control plants relative to the unfertilized plants (p< 0.05; lower case letters, Table 1; Supplementary Figure S1). A 2.5 and 3.6-fold change was observed after fertilization in mowed and control plants, respectively (Table 1). When comparing among termination methods, a significantly lower abundance of amoA_AOB was observed for glyphosate termination relative to both mowing and control within the fertilized condition (Table 1, upper case letters), a result that was not observed within the unfertilized condition.

A statistically significant main effect of N fertilization was observed for PNA (p<0.001). A greater ammonia oxidation activity was observed in N fertilized treatments, regardless of the termination method applied (Table 2). In contrast, termination method had no influence on this functional assay (p>0.05). Further, these results were consistent with the finding from the addition of root exudates from glyphosate-treated plants or mowed plants to soil slurries. The addition of exudates did not show a significant effect in the nitrifying response measured by PNA (p> 0.05; Supplementary Figure S2).

When analyzing abundance-activity correlations, a strong positive correlation was observed between PNA and amoA_AOB copies (r=0.78, p< 0.001), while no significant correlation was observed between amoA_AOA copy number and PNA showing a negative and weak correlation (r=−0.19, p> 0.05).

The total number of amoA_AOA and amoA_AOB sequences in each sample after each step (filtering, denoising, and chimera removal) is shown in Supplementary Tables S2, S3, respectively. As shown in these tables, each step was conducted without removing a large number of reads, reflecting the high quality of the datasets. For AOA, a total of 216 ASVs were detected with a final length of 283bp. A total of 113 ASVs were obtained for AOB with a final length of 452nt. The amoA_AOA and amoA_AOB amplicon sequencing datasets have been deposited in Sequence Read Archive (SRA) repository under the accession PRJNA701453. Rarefaction curves (Supplementary Figure S3) indicated that sequencing depth was enough to cover the full diversity of AOA and AOB, even at the rarefaction value (4,543 for AOA and 2,971 for AOB).

After clustering ASVs at the appropriate amoA identity threshold (90%), we obtained OTU tables for AOA (15 OTUs) and AOB (6 OTUs) with a considerably lower sparsity than the zero-enriched ASV tables, rendering them appropriate for further beta diversity analyses with multivariate methods and compositional data. The BLASTn search using the nucleotide database reported significant hits with amoA sequences of uncultured AOA (Supplementary Table S4) and uncultured AOB (Supplementary Table S5).

For AOA, no significant differences were observed for any of the alpha-diversity metrics (p> 0.05) between N-fertilized and unfertilized plants or among termination methods (Supplementary Table S6). Statistically significant greater values of alpha diversity metrics were observed for AOB in the rhizosphere of N fertilized plants compared to the unfertilized (p= 0.003 for Shannon diversity index H′ and p= 0.013 for the observed richness S′), except for Pielou’s evenness index (J’; p= 0.19). Conversely, no significant differences among termination methods were observed for any of the parameters studied (p> 0.05, Supplementary Table S7).

The multivariate statistical analysis through PERMANOVA to investigate the beta diversity of AOA indicated no interaction between fertilization and the termination method (Table 3). A marginally significant effect of fertilization was detected on AOA communities when the dataset was analyzed under the standard approach with the phylogenetic distance (UniFrac), while no significant effect of fertilization was observed with Aitchison distance (p> 0.05, Table 3). Contrary to AOB, a statistically significant effect of the termination method was observed on the community structure of AOA (p< 0.05, Table 3) for both Aitchison distance and the UniFrac distance. Pairwise PERMANOVA indicated significant differences between communities of glyphosate-treated and mowed plants (FDR adjusted p= 0.04) or between mowed plants and the untreated control (FDR adjusted p= 0.04). No significant differences were detected between glyphosate and the untreated control (FDR adjusted p> 0.05), as visualized in NMDS analyses based on Aitchison distance (Figure 1) and generalized UniFrac distance (Supplementary Figure S4).

Principal component analysis of AOA (Figure 2) showed that the first two PCs explained 53.5% of the total variance of the dataset. A separation across PC1 was observed between communities of mowed plants and the other treatments (Figure 2), as previously observed in NMDS analyses (Figure 1). The two-way ANOVA with PCs as variables indicated no significant interaction and a statistically significant effect of termination method on PC1 (p< 0.05, Table 4). Pairwise comparisons revealed no differences between glyphosate termination and the untreated control, while mowing termination was significantly different from both treatments (p< 0.05, Tukey’s HSD test). While an interaction between fertilization and termination method was observed for PC2 (Table 4), no significant differences were detected between fertilized and unfertilized samples within each termination method, nor among termination methods within fertilized or unfertilized samples (p> 0.05, Tukey’s HSD test). The contributions of OTU5 and OTU6 to the ordination on PC1 were notably greater than the contributions of the other OTUs (53.5 and 35.8%, respectively). A positive correlation was observed between OTU5 and PC1, while a negative correlation with this component was detected for OTU6 (Table 4). As shown in the biplot, OTU5 was enriched in communities from glyphosate-treated or control plants and depleted in the rhizosphere of mowed plants, the opposite trend was observed for OTU6 (Figure 2). Further phylogenetic analyses (maximum likelihood tree, Figure 3) indicated that the sequence of responsive OTU5 is related to the amoA sequence of an uncultured AOA isolated from a suboxic soil (red label, Figure 3). No clustering with the other OTUs in this study or with amoA sequences from maize rhizosphere, bulk soil, or sediments was observed for OTU6 (green label, Figure 3).

For AOB, the statistical analysis of Aitchison and UniFrac distances through PERMANOVA indicated no interaction between fertilization and the termination method (p> 0.05, Table 3). N fertilization showed a marginally significant effect on the community structure of AOB for the Aitchison distance (Figure 4A) and a statistically significant effect for the phylogenetic distance (generalized UniFrac; Table 3 and Figure 4B). No significant effect of the termination method was detected on the community structure of AOB (PERMANOVA, p> 0.05, Table 3). The phylogenetic analysis of amoA_AOB amplicons (maximum likelihood tree) showed most OTUs within the Nitrosospira cluster 3 (green and light green labels, Figure 5) and only one OTU (OTU4) closer to Nitrosospira cluster 0 (light blue label, Figure 5). No Nitrosomonas spp. were detected in the rhizosphere of A. sativa plants grown in this soil.