Long-Term N Fertilization Decreased Diversity and Altered the Composition of Soil Bacterial and Archaeal Communities

: Soil microbial communities are essential in the cycling of nutrients that a ﬀ ect crop production. Our goal was to characterize the microbial community structure following 34 years of nitrogen (N) fertilization treatments in continuous maize production in highly fertile soils. Using 16S rRNA gene-based analysis of the V4 region via Illumina HiSeq2500 technology with downstream bioinformatics processing and analysis with QIIME 2.0, we aimed to characterize the prokaryotic communities under three increasing N fertilization rates. Factor analyses indicated that a high N level decreased the diversity of soil bacterial and archaeal communities and altered the relative abundance (RA) of the dominant ( > 1% RA) and minor ( < 1% RA) phyla. Among the 12 major phyla, we determined increases in Gemmatimonadetes, Proteobacteria, and Euryarchaeota, accompanied by reductions in Cyanobacteria, Chloroﬂexi, Firmicutes, and Planctomycetes with increasing N. Within the 29 minor phyla, N fertilization led to increases in Aquiﬁcae, WPS2, Parvarchaeota, AD3, FCPU426, Armatimonadetes, TM7, Chlamydiae, and OD1, along with reductions of Nitrospirae, WS3, Tenericutes, Lentisphaerae, OP3, Synergistetes, Thermotogae, and prokaryotes that could not be reliably assigned to a phylum (classiﬁed as Other).


Introduction
Nitrogen fertilizers have been a major contributor to the impressive crop yield increases realized since the 1950s [1], widely used to maintain soil fertility and crop productivity [2]. Soil microbial communities provide many benefits including assimilation of nutrients, plant disease resistance, and stabilization of the soil ecosystem [3,4]. In agronomic systems, these functions positively shift the soil quality and productive capacity over the long term [5]. Therefore, investigating the effects of nitrogen (N) fertilizers on soil microbial communities is of vital importance due to the critical part that microbes play in biochemical processes. Nitrogen additions directly affect plant residue quantity, quality, and rate of decomposition, thus shifting soil nutrient content [6] which can influence the size, structure, and function of microbial communities. Modern approaches, using phylogenetic surveys of bacteria employing universal primers, allow for characterization and comparison of the microbial diversity in different soil environments as well as the comparative analysis of changes in community structure due to management practices [7,8]. Although numerous studies have examined the influence of N fertilizer on the soil microbiome [9][10][11][12], most of these studies compared responses to organic and inorganic fertilizers; only a few studies examined the response of the soil microbiome to increasing N fertilizer rates in an agroecosystem [11,[13][14][15][16][17].
Changes in nutrients and type of fertilizer (mineral vs. organic) can cause a shift in bacterial communities in soils subjected to fertilization, favoring more active, copiotrophic microbial communities [14,18]. Intensive agricultural management practices, including fertilization levels, are generally considered to lead to simpler, less diverse, soil food webs [19] with varied impact on microbial biomass [20,21]. Therefore, nitrogen enrichment, as shown in many experiments [22], usually leads to declines in microbial community diversity and species richness. In addition to diversity, phylogenetic shifts of the soil microbial community have also been observed in response to N fertilization [14]. In a recent meta-analysis, Geisseler and Scow [20] researched these influential factors using a dataset based on 64 long-term field trials from across the world. Urea or ammonium salts were the most commonly used fertilizers and the durations of these trials ranged from 5 to 130 years, while the N fertilization rates ranged from 10 to 650 kg N/ha, averaging 136 kg N/ha. The authors concluded that long-term repeated mineral N applications can alter microbial community composition.
Well-maintained long-term field experiments can provide insights into the long-term effects of agricultural practices on soil properties [23]. Field trials of N fertilization rates on continuous corn in this study were established in 1981, which presents a unique opportunity to characterize the bulk soil microbiome and to test whether microbial shifts occur under increasing N fertilization levels. An improved understanding of the microbial community can help translate beneficial effects on the soil environment such as reducing nitrate runoff and greenhouse gas emission [24,25], into actionable practices. However, baseline characterization of the soil biological component and microbial shifts under N fertilization is currently lacking in Illinois. Therefore, the objective of this study was to characterize the prokaryotic community structure under long-term N rates applied to continuous corn plots using DNA sequencing and bioinformatics analyses.

Site Descriptions and Soil Sampling
The detailed set up for this field trial was described previously [25]. This study was established in 1981 at the Northwestern Illinois Agricultural Research and Demonstration Center in Monmouth, IL (40 • 90 N, 90 • 73 W) to study the effects of N fertilization on crop yields. The soils are Muscatune series (fine-silty, mixed, superactive, mesic, Aquic Argiudolls) with surface soil being usually dark, greyish brown and granular in structure. The field trial consisted of three N fertilization rates (0, 202, and 269 kg N/ha) in continuous corn (CCC) arranged in a randomized complete block design (RCBD) with three replications. These rates represent the current average (202 kg N/ha) annual application of N and what can be considered an overapplication of N (269 kg N/ha) for continuous corn production. Total experimental units include 3 × 3 = 9 plots. The main plots in this trial were 12 m × 9 m. We continued the agronomic practices that have been in use in this long-term trial [25]. Tillage consists of primary tillage with a chisel plow 20-25 cm deep in the fall after harvest, and secondary tillage with a field cultivator before planting in the spring. Corn is planted in April or early May each year, in 76-cm rows at a density of about 75,000 seeds/ha. The N fertilizer at 0, 202, and 269 kg N/ha rates for corn was injected at 10-15 cm depth before planting using urea ammonium nitrate solution (28% N). Weed control was achieved by using recommended rates and timing of appropriate herbicides along with hand weeding if necessary. Plots were maintained weed free during the study period.
Soil sampling of this study was conducted immediately after harvest and before fall tillage in 2015, following 34 years of N rate treatments in the same plots each year. Three composite subsamples (~500 g each) per plot were collected to 10 cm depth using an Eijkelkamp grass plot sampler. Each composited subsample was collected walking in a zig-zag pattern and contained about 10 pushes of the soil sampler probe. Samples were preserved with ice packs in the field, and frozen to −20°C upon arrival to our lab facilities. Results from laboratory analyses and determinations were averaged to get one value per plot per block, thus n = 3 for each N rate treatment.

Soil Characterization
After determining gravimetric water content for each soil sample, field moist soil was analyzed for available N, NO 3 − , and NH 4 + (mg/kg) using KCl extraction (

Soil DNA Extraction, qPCR, and Sequencing
Total genomic DNA from the soil samples was extracted employing the Power Soil DNA isolation kit (MoBio Inc., Carlsbad, CA, USA) from 0.25 g of the composited soil samples, with the soils carefully homogenized before subsampling. The extracted DNA quantity and quality were checked with a Nanodrop 1000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA) following the manufacturer's protocol. The extracted DNA samples were stored at −20 • C until analysis. An amplicon library with individual barcodes for each sample and compatible with the Illumina HiSeq was constructed. To build the amplicon library, 25 µL PCR reactions were performed using a BioRad T100 thermal cycler in 25 µL volumes with 1× buffer (GoTaq ® Flexi buffer; Promega Corp.), with the following composition: 2.5 mM MgCl 2 , 200 µM dNTPs, 0.4 µM each primer (forward and reverse), 1.0 µL template DNA (pooled amplicons), and 1.0 unit of GoTaq polymerase. PCR parameters were: initial denaturation at 95 • C for 10 min, followed by 34 cycles of amplification (45 secs at 95 • C; 45 s at 58 • C; 45 s at 72 • C), and a final extension at 72 • C for 10 min. PCR products were visualized on a 1.3% agarose gel containing GreenGlo™ Safe DNA dye (Denville Scientific, Inc. Metuchen, NJ, USA) under UV illumination. The 16S rRNA gene (V4 region) was amplified using the primer set 515F (GTGYCAGCMGCCGCGGTAA) and 806R (GGACTACNVGGGTATCTAAT) [26]. The primers were designed as 5 -PCR-specific + gene region + 3 -PCR-specific + 10 nt barcode and the Fluidigm platform utilized two primer sets simultaneously to create the final DNA amplicon.
The resulting amplicon libraries were quantitated with a Qubit Fluorometer and run on a bio-analyzer to evaluate the profile of fragment lengths. The barcoded libraries were pooled in equimolar concentrations and diluted to 10 nM. The diluted libraries were sequenced at the Roy Carver Biotechnology Center Functional Genomics lab at the University of Illinois at Urbana-Champaign (Urbana, IL, USA) using paired-end sequencing on the Illumina HiSeq2500 (Illumina, San Diego, CA, USA) in rapid mode yielding reads 250 nt in length.

Bioinformatic Analysis
The prokaryotic community composition and diversity were analyzed using the V4 region of 16S rRNA gene sequences. Quality control of the paired-end sequences included trimming of low-quality bases at the extreme of the sequences and filtering of sequences that had an average PHRED quality score <20, sequence length <200 nt or >800 nt, presence of ambiguous nucleotides, primer mismatch, and errors in the barcode. The primers on both forward and reverse reads were trimmed using the software Trimmomatic-0.36 [27]. The overlapping paired-end reads were joined using the fastq-join program [28]. The quality of reads was assessed using FASTQC. The chimera checking of the sequences was performed with USEARCH v5.2 [29] using GOLD as a reference database [30] and sequences identified as chimeric were filtered. The high-quality sequences were assigned to samples according to the barcodes. Sequences were clustered into operational taxonomic units (OTUs) using open reference OTU picking at 97% similarity threshold and clustered against the Greengenes 16S rRNA database [31].
The most abundant sequence of each OTU served as a representative sequence and was aligned using the Python Nearest Alignment Space Termination (PyNAST) algorithm with a minimum 80% identity [32]. Annotation of the taxonomic unit to representative OTU was based on classifications from the Ribosomal Database Project (RDP) and an OTU table was generated. Downstream processing of sequences was completed using the python programs built in the QIIME (Quantitative Insights into Microbial Ecology) pipeline v1.9.1 [33]. The rarefaction plot of the three treatments was generated following the QIIME tutorial protocols. The generated OTU table was rarefied to an even depth (135,639 sequences per sample) using the single_rarefaction.py script, thus all samples were retained for further analyses. The OTUs with a number of sequences <0.005% of the total number of sequences were discarded using filter_otus_from_otu_table.py script. Alpha diversity metrics including chao1 [34], Shannon's index, Abundance-Based Coverage Estimator (ACE) and Fisher's alpha [35] were extracted using alpha_diversity.py in the QIIME pipeline and compared among N fertilization rates as detailed in the statistical analysis subsection. To assess overall variation in community structure among treatments (i.e., beta-diversity), rarefied data from metabarcoding profiles were subjected to multivariate analysis using weighted UniFrac distances [36]. Ordination was performed using a principal coordinate analysis (PCoA) through_plot.py script in the QIIME pipeline to create data visualization of the relationships between samples from three increasing N fertilization rates. Non-parametric testing of beta diversity was performed employing jackknife resampling based on 5000 sequences per sample. A permutational multivariate analysis of variance using distance matrices to assess beta diversity was conducted in ADONIS (R built in to compare categories.py in the QIIME pipeline) [37].

Statistical Analysis
Linear mixed models were fit to each of the measured soil properties of pH, CEC, SOM, and available N forms using PROC GLIMMIX in SAS 9.4 (SAS Institute Inc., 2012); using the N rates and block factors as fixed and random terms, respectively. To analyze the prokaryotic community changes under increasing N fertilization rates, we deployed multivariate techniques followed by linear mixed models in SAS. The data set included 41 variables, which represented the relative abundances (RA, %) of each phylum extracted with the bioinformatics pipeline. Pearson's correlation coefficients were calculated using the CORR procedure in SAS to explore the correlation among soil prokaryotes (data not shown). Correlations between variable pairs were found to be >|0.3| (moderate to high range) in most cases which indicated the need to deploy a data reduction technique such as principal component analysis (PCA) to avoid problems of multicollinearity by compiling the information into a new smaller set of uncorrelated variables. PCA was performed using the PROC FACTOR procedure with the priors = 1 option [38]. The analysis resulted in 13 principal components (PCs) with eigenvalues >1 together explaining 88% of the variability of the data set. Only those PCs that in addition explained >5% of the variability were included in subsequent analyses. Next, linear mixed models were fit to the PCs extracted using PROC GLIMMIX to evaluate the effect of N fertilization rates on the RA of the main phyla now summarized within the PCs. Blocks were considered a random effect and N rates as fixed effects. When appropriate, least square means were separated using the lines option of the lsmeans statement, setting the probability of Type I error at 0.05(α). Figure 3 was built in Sigma Plot 12.5 (Systat Software, Inc.) plotting the standard error of treatment means as inferential error bars [39]. Table 1 shows the results for the general soil characteristics of pH, CEC, SOM, and available N forms (NO 3 − and NH 4 + ) under continuous corn production and in response to increasing N fertilization rates following 34 years of treatments. Results showed a gradual and statistically significant reduction in soil pH with increasing N fertilization rates, with a parallel increase in CEC. Soil organic matter (SOM) was not found to be affected by N additions, yet both available N forms showed statistically significant responses to fertilization. Nitrate levels (NO 3 − ) increased under the maximum N rate with no statistical differences detected between the 0 and the 202 kg N/ha fertilization rates. On the other hand, NH 4 + concentration showed a gradual increase with increasing N rates.

Sequencing Summary
After the elimination of low-quality sequences using Trimmomatic and chimeras using USEARCH, a total of 6,459,949 bacterial 16S rRNA sequences were obtained with a range of 107,802 to 350,409 and a mean of 239,257 sequences per sample. To measure the quality and sufficiency of sequences in the samples of each treatment, a rarefaction plot was generated ( Figure 1). In this study, all rarefaction curves tended to approach the saturation plateau indicating that only the rarest species remained to be sampled, and thus the sequences obtained were enough to capture the diversity information of the prokaryotic community within our samples.

Sequencing Summary
After the elimination of low-quality sequences using Trimmomatic and chimeras using USEARCH, a total of 6,459,949 bacterial 16S rRNA sequences were obtained with a range of 107,802 to 350,409 and a mean of 239,257 sequences per sample. To measure the quality and sufficiency of sequences in the samples of each treatment, a rarefaction plot was generated ( Figure 1). In this study, all rarefaction curves tended to approach the saturation plateau indicating that only the rarest species remained to be sampled, and thus the sequences obtained were enough to capture the diversity information of the prokaryotic community within our samples.

Soil Prokaryotic Community Diversity
After the upstream bioinformatics analysis, a rarefied and filtered OTU table was generated using QIIME (data not shown). Table 2 shows the observed OTUs, phylogenetic diversity, and other estimators of alpha diversity and richness of the soil prokaryotic community. Statistically significant reductions were determined in the number of observed OTUs and all richness (Chao1 richness estimate and ACE diversity index) and diversity (phylogenetic diversity (PD), Fisher's alpha, and Shannon's diversity) parameters at the highest fertilization level (269 kg N/ha) compared to the

Soil Prokaryotic Community Diversity
After the upstream bioinformatics analysis, a rarefied and filtered OTU table was generated using QIIME (data not shown). Table 2 shows the observed OTUs, phylogenetic diversity, and other estimators of alpha diversity and richness of the soil prokaryotic community. Statistically significant reductions were determined in the number of observed OTUs and all richness (Chao1 richness estimate and ACE diversity index) and diversity (phylogenetic diversity (PD), Fisher's alpha, and Shannon's diversity) parameters at the highest fertilization level (269 kg N/ha) compared to the parameters calculated from the control (0 kg N/ha) and the intermediate rate of 202 kg N/ha ( Table 2). To characterize the overall effect of N rates on the beta diversity of the total bacterial phyla, PCoA was applied to calculate the pairwise distances between the communities at 0, 202, and 269 kg N/ha N fertilization levels. The results of PCoA are shown in Figure 2 demonstrating a contrasting separation between the three N fertilization levels.  Table  2). To characterize the overall effect of N rates on the beta diversity of the total bacterial phyla, PCoA was applied to calculate the pairwise distances between the communities at 0, 202, and 269 kg N/ha N fertilization levels. The results of PCoA are shown in Figure 2 demonstrating a contrasting separation between the three N fertilization levels. The three axes of the PCoA respectively explained 26.81%, 5.66%, and 4.50% of the total variance. Furthermore, the results of the permutational multivariate analysis of variance using weighted UniFrac distance matrices showed that the increasing N fertilization levels had a significant effect on the dispersion of beta diversity within the soil archaeal and bacterial communities (R 2 = 0.34, p < 0.001).

Soil Prokaryotic Community Composition
Forty-one different phyla including three archaea phyla were identified from the samples. The community structure expressed as relative abundances (RA, %) of the 41 phyla following 34 years of N fertilization at 0, 202, and 269 kg N/ha is shown in Table 3 and 4. Table 3 shows those phyla with >1% RA (12 groups) together accounting for about 95% of the total OTUs. Ordered from highest to lowest abundance, Table 4 shows the remaining 29 phyla detected in our samples that had <1% RA and accounted for the remaining 5% of the RA. The three axes of the PCoA respectively explained 26.81%, 5.66%, and 4.50% of the total variance. Furthermore, the results of the permutational multivariate analysis of variance using weighted UniFrac distance matrices showed that the increasing N fertilization levels had a significant effect on the dispersion of beta diversity within the soil archaeal and bacterial communities (R 2 = 0.34, p < 0.001).

Soil Prokaryotic Community Composition
Forty-one different phyla including three archaea phyla were identified from the samples. The community structure expressed as relative abundances (RA, %) of the 41 phyla following 34 years of N fertilization at 0, 202, and 269 kg N/ha is shown in Tables 3 and 4. Table 3 shows those phyla with >1% RA (12 groups) together accounting for about 95% of the total OTUs. Ordered from highest to lowest abundance, Table 4 shows the remaining 29 phyla detected in our samples that had <1% RA and accounted for the remaining 5% of the RA.  Thus, Table 3 shows that members of the phyla Proteobacteria, Acidobacteria, and Bacteroidetes dominated all treatment levels, with these three groups accounting for over 58% of the 94% total OTU abundance, collectively represented by the 12 phyla listed therein. On the other hand, for the minor groups listed in Table 4, the highest representation (1%-0.5% RA) within the 0 kg N/ha was determined for WS3, Nitrospirae, and Armatimonadetes, together accounting for 2.4% of the remaining 5.3% total OTUs. Phyla Parvarchaeota, Armatimonadetes, and WS3 were also within the top groups with an intermediate rate of N added (202 kg N/ha) and Aquificae appeared within the top three groups at the highest N level (269 kg N/ha) following Parvarchaeota and Armatimonadetes.
To avoid issues of multicollinearity in our results we conducted a PCA on the relative abundances of the 41 main phyla. The PCA on 41 variables rendered a set of 13 uncorrelated variables (PC 1 to PC 13) with eigenvalues larger than 1 that explained about 88% of the variability contained in the main phyla database. Only the initial five PCs extracted by this technique contributed more than 5% of explained variance and are thus presented in Table 5.  The first PC1 had the largest eigenvalue (11.99) and explained about 29% of the community taxonomic variability while including both dominant and minor taxa (Table 5). Among the dominant taxa, PC1 included high positive loadings (>0.5) for Chloroflexi, Cyanobacteria, Firmicutes, and Planctomycetes, and high negative loadings (<−0.5) for Gemmatimonadetes, Proteobacteria, and Euryarchaeota. For the minor taxa, PC1 included high positive loadings (>0.5) for Nitrospirae, WS3, Tenericutes, Lentisphaerae, OP3, Synergistetes, Thermotogae, and prokaryotes that could not be reliably assigned to a phylum (grouped as "Other"), as well as high negative loadings (<−0.5) for Aquificae, WPS2, Parvarchaeota, AD3, FCPU426, Armatimonadetes, TM7, Chlamydiae, and OD1.

Microbial group Eigenvalues
As shown in Table 5, PC2 had an eigenvalue of 3.81 and explained an additional 10% of the variability in the dataset; PC2 eigenvector had high positive loadings for Chlamydiae, Chlorobi, Thermi, and TM7, and high negative loadings for AD3, Aquificae, and Elusimicrobia. The eigenvalue for PC3 was 3.37 and explained an additional 8% of the total variability. The PC3 eigenvector included high positive loadings for Acidobacteria, OP11, and Spirochaetes while including negative loadings for Actinobacteria and Choloroflexi ( Table 5). The eigenvalue for PC4 was 2.75 and explained an additional 6% of the total variability. The PC4 eigenvector included high positive loadings for Crenarchaeota, and Elusimicrobia while including a negative loading for Bacteroidetes. The eigenvalue for PC5 was 2.55 and explained an additional 7% of the total variability. Lastly, the PC5 eigenvector included a high positive loading for EM19 and negative loadings for Other and Verrucomicrobia. Each of these PCs were used as independent variables in follow up analysis of variance (ANOVA) to test the effect of N fertilization levels on the RA of the phyla represented by each PC. The probability values and degrees of freedom associated with the ANOVA for the effects of N rate are shown in the lower portion of Table 5 for each extracted PCs. The N fertilization effect was statistically significant for PC1 but not statistically significant for any of the remaining PCs.
Thus, PC1 mean value under no N fertilization (1.12) was statistically and significantly different (LSD 0.05 = 0.41) from the mean value measured under 202 kg N/ha (−0.9), which was also statistically different that the PC1 mean value following 34 years of N fertilization at a rate of 269 kg N/ha (−1.13). Thus, the phyla represented within PC1 showed a clear and gradual response to the increased N application. Scheme 1 depicts the changes in RA% of dominant and minor prokaryotic phyla with increasing N rates following 34 years of management as summarized by the measured changes in PC1 scores. Specifically, for the Proteobacteria phyla, Figure 3 shows the changes in the RA of the main classes under increasing N additions. Thus, Alphaproteobacteria and Gammaproteobacteria showed gradual and statistically significant increases in their RAs with each additional N rate (p < 0.0001 and p = 0.0028, respectively). While the RA of Deltaproteobacteria decreased with the rates of 202 kg N/ha and 269 kg N/ha (p = 0.0009), Betaproteobacteria and Epsilonproteobacteria were not affected by N fertilization. Scheme 1. An illustration of changes in the relative abundance (RA) of dominant (increased size font) and minor (regular font size) prokaryotic phyla brought about by increasing N fertilization levels after long-term N fertilization management. Phyla in green font show increases in their RA, while phyla in blue font show decreased RA with increasing N.

Discussion
A number of previous studies have revealed that N application level significantly influences microbial biomass, community composition, and function [14,15,40]. While many studies on 25 grassland sites across the world detected changes in soil bacterial communities across the N gradient over the long-term [41], research on Illinois agricultural field sites was missing in this regard. Using a long-term fertilization field experiment established in 1981, we monitored the structure of soil prokaryotic communities. We chose sequencing the V4 locus of 16S ribosomal RNA genes of prokaryotes in order to characterize the community shifts associated with long-term N-fertilization treatment. While the analysis of the response of microbial communities within this low level of taxonomic resolution (phylum) is done under the (false) assumption that all members of a phylum have similar biochemical and metabolic capabilities, it is however a convenient and commonly adopted approach to look for patterns within a manageable number of groups. Moreover, despite the enormous amount of phylogenetic and physiological diversity within the main bacterial phyla, functional diversity may be broadly mapped across major phyla of Bacteria and Archaea [42,43].
The highest rate of N fertilization caused a decline in all alpha diversity indexes after a 34-year treatment, which coincides with other research studies [43,44]. Direct effect of increased N input on prokaryotic community structure may be linked to the indirect effect of pH changes when soil nutrient availability and pH are concomitantly altered by fertilization [16]. Ding et al. [44] observed significant reduction of alpha diversity indexes in N P K fertilization treatment compared to unfertilized soil cropped with maize-wheat-soybean rotation. These contrasting changes are mainly associated with reduced soil pH [43], due to the long-term application of urea and ammonium fertilizers [45]. Likely, high N levels in mineral fertilizer have resulted in intensified soil acidification through nitrification of the ammonium applied [17]. Intermediate levels of N may not reflect pronounced shifts in soil pH, given that some of the acidity produced by nitrification is neutralized when plants take up more nitrate than cations [46]. Our results demonstrate a similar point, in that there is a significant difference between the intermediate and high N rate group on all alpha diversity and richness indicators, while no significant differences were found between the control and intermediate N rate group (Table 2). Shifts in the prokaryotic community structure could be attributed not only to pH changes but also to the change of the C/N ratio brought about by long-term N fertilization [47].
Results from our PCA showed that soil prokaryotic communities significantly differed among the three levels of N input, along the PC1. The variability among samples was explained by seven dominant taxa (with RA >1%) and 15 minor taxa (RA <1%). Therefore, we focused the analysis of composition shifts in that subset of phyla. Phyla showing negative correlations with PC1 are those that increased in RA under increased N input. Conversely, the phyla that have positive correlations with the PC1 diminished in their RA with increasing N levels. We found that the relative abundance of 11 bacterial and two archaeal groups increased with increasing N rates. The major bacterial group affected was Proteobacteria. This phylum comprises more than one third of the species that have been characterized up to date [42], and their members display enormous metabolic diversity. Increased nutrient inputs stimulate the growth of potential copiotrophic bacteria, while impairing oligotrophs [14]. It is of no surprise, then, that the overall group responded to increasing levels of N as most members of this group fall in the copiotrophic and generalist lifestyle, demanding high nutrient levels [48]. Moreover, the main steps in the biogeochemical N cycle such as N fixation, denitrification, and nitrification, are driven by Proteobacteria [42]. At the finer level of resolution of classes, N fertilization increased the relative abundances of the so-called copiotrophic Alphaproteobacteria and Gammaproteobacteria, while Deltaproteobacteria diminished and Betaproteobacteria and Epsilonproteobacteria were unaffected by N inputs (Figure 3). Our results are in agreement with those reported by Li et al. [43], in soil cultivated with maize with different N levels, as well as with similar patterns observed in soil under wheat-maize-soybean crop rotation and fertilized with inorganic N [44].
The phylum Gemmatimonadetes was the second largest group responding positively to N fertilization. This phylum is one of the most commonly detected phyla in 16S rRNA gene libraries from soil [49]. Other authors have previously reported increased relative abundance of Gemmatimonadetes in fertilized soils of distinct origin [50][51][52]. Under changing environments, it has been proposed that these microorganisms are able to fine-tune their C and N intake according to their metabolic needs, suggesting a generalist ecological strategy [52]. In that way, Gemmatimonadetes may be able to take advantage of the higher nutrient input provided with the highest N rate.
Aquificae and Armatimonadetes were among the bacterial groups which positively responded to long-term inorganic N input. Aquificae harbors some species that display autotrophic, chemolitotrophic, and microaerophilic growth, and some of them are able to use nitrate as an electron acceptor conducive to denitrification [42], while in previous studies, Armatimonadetes responded positively to the addition of inorganic NPK fertilizer [50]. These responses might help explain their higher abundance under high N levels in the present study.
Other bacterial groups that comprised a small proportion of total diversity, with RA<0.1% also responded positively to N input, namely Chlamydia, FCPU-426, WPS-2, AD3, TM7, and OD1 (Parcubacteria). Phylum Chlamydia comprises few species within a single order Chlamydiales, all of which are obligate intracellular parasites of eukaryotes, and are ubiquitous in nature within free living amoeba [42]. The rest of those groups are candidate phyla, with no cultivated representatives, that are collectively designated as the "microbial dark matter"; they remain to be characterized [53]. However, they have been detected previously in fertilized soils at comparable relative abundances (<1%) [49].
Both Euryarchaeota, which showed RA slightly higher than 1%, and Parvarchaeota (RA <1%), responded positively to N input. The Euryarchaeota is the largest archaeal phyla and the most metabolically diverse. This group includes strictly anaerobic methanogens; obligate aerobic extreme halophiles, and the acidophiles Thermoplasma and Ferroplasma, that thrive in low pH environments growing as chemolithotrophs oxidizing ferrous iron [42]. Members of Euryachaeota are capable of both N 2 fixation and denitrification, so it is not surprising that their relative abundance increases under higher N input. Parvarchaeota members are characterized by their small size and limited metabolic capabilities, as they lack biosynthetic pathways for all amino acid synthesis, and thus, they must scavenge them from the environment contributing to the N cycling. Studied representatives have been shown to closely interact with Thermoplasmatales (Euryarchaeota) which possibly act as a host providing amino acids and thus explain the co-occurrence in our study and others [54].
Groups that decreased under increasing N levels belong to four dominant bacterial phyla (Chloroflexi, Cyanobacteria, Firmicutes, and Planctomycetes) and eight minor taxa with RA <1%. Chloroflexi, Cyanobacteria, and Firmicutes are frequent inhabitants in the soil environment, contributing to the main steps of the N cycle (N fixation and denitrification), while members of the Planctomycetes phylum are responsible for the anaerobic ammonia oxidation. Their RA in our study where comparable to the reported values in Chinese black soils during the maize growing season [43]. Relative abundances of Cyanobacteria [55] and Planctomycetes [50] were reported to decrease with inorganic N fertilization compared to unfertilized control or organic fertilization. Both Chloroflexi and Planctomycetes were reported among the most abundant phyla in soils of the Argentine Pampas, with higher RA in pristine soil in comparison to fertilized agricultural soil [52].
Among the less abundant taxa (RA <1%) the phylum Nitrospirae was also negatively affected by increasing N inputs. This phylum comprises species that oxidize nitrite autotrophically, and have been classified as K-strategist, oligotrophic bacteria with a high affinity for nitrite scavenging. Thus, it can be argued that high levels of N fertilization may negatively affect the growth of this group by favoring the growth of other competing nitrite oxidizers, such as Nitrobacter with an r-strategist lifestyle [53]. The RA of candidate phylum Latescibacteria (WS3) was highest in unfertilized soils (0.98%). In contrast to our results, Latescibacteria sequences showed higher RA in fertilized agricultural soil samples from the Pampas region (4%) in comparison to pristine soil [9,52]. Both phyla have been shown to correlate positively to soil pH [16], thus the high N ratio might have indirectly reduced their RA by decreasing soil pH. Relative abundances of Tenericutes, Lentisphaerae, OP3, Synergistetes, and Thermotogae also decreased with increasing N fertilization, yet their RA was well below 1%. Little is known about these groups and their physiological and ecological roles, so discussion of the observed results is limited by the available information.
Usually, as it occurs in our study, data reporting prokaryotic community structure refer to relative, not absolute, abundances. As such, decreases in the relative abundance of certain groups could mean no change at all in the absolute number of these taxa, especially right after incorporation of crop residues and/or nutrient inputs that may greatly increase the abundance of fast-growing copiotrophic prokaryotes but not necessarily result in reduced abundance of other groups [47]. This present study was conducted after continuous corn was fertilized annually for 34 years, however, soil sampling was done immediately after harvest. Thus, we can safely attribute the detected changes in prokaryotic community structure among treatments to the long-lasting effects of fertilization on the soil environment.
Although we are aware that the prokaryotic community response may be better assessed by refining the analysis and examining bacterial and archaeal composition at a lower taxonomic level (at genus or OTU level), the broad taxonomic level used in this first exploratory analysis was efficient and sensitive enough to disclose the large impact of long-term N fertilization on overall prokaryotic diversity.
Our results suggest that main responding taxa may be those driving key steps in the N cycle. This hypothesis is partially supported by a recent research conducted to quantify genes encoding N cycle functions in the same field trial [13]. However, it should be confirmed by our future work towards prediction of functions based on 16S rRNA gene sequences by in silico analysis (with tools such as PiCrust or Tax4Fun).

Conclusions
Long-term N fertilization treatments have resulted in differences in the structure of soil prokaryotic communities modifying the relative abundance of archaeal and bacterial groups. A more complete understanding of the long-term interactions and feedback between fertilization practices, soil properties, and diversity of prokaryotic communities is essential for enabling the design of optimal agronomic practices that improve agricultural sustainability.