Bacterial Community in Soils Following Century-Long Application of Organic or Inorganic Fertilizers under Continuous Winter Wheat Cultivation

Fertilization is one of the most common agricultural practices to achieve high yield. Although microbes play a critical role in nutrient cycling and organic matter decomposition, knowledge of the long-term responses of the soil bacterial community to organic and inorganic fertilizers is still limited. This study was conducted to evaluate the effects of century-long organic (manure), inorganic (NPK), and no fertilization (control) treatments on soil bacterial community structure under continuous winter wheat (Triticum aestivum L.) cultivation. Fertilization treatments altered the richness, diversity and composition of the soil bacterial community. Compared with the control, manure significantly increased the operational taxonomic units (OTUs), Chao 1 and Shannon indices, and taxonomic groups, while NPK significantly decreased these parameters. Fertilization treatments did not alter the types of dominant phyla but did significantly affect their relative abundances. Acidobacteria and Proteobacteria were the most dominant phyla in all treatments. Manure led to enrichment of most phyla, with a diazotrophic group, Cyanobacteria, being an exception; NPK reduced most phyla, but enriched Chloroflexi; control led to promotion of Cyanobacteria. Soil pH and NO3− were two dominant parameters influencing the bacterial community structure. Soil pH positively correlated with the relative abundances of Proteobacteria and Gemmatimonadetes but negatively correlated with those of Acidobacteria and Chloroflexi; NO3− negatively correlated with the relative abundance of Cyanobacteria, which was 14–52 times higher in control than the fertilized soils. Cyanobacteria, especially M. paludosus and L. appalachiana, could be the key players in maintaining wheat productivity in the century-long unfertilized control.


Introduction
Understanding the soil microbial community structure presents understated challenges, largely due to the complex nature of the system and the enormous abundance and diversity of the community residing within [1]. Bacteria comprise over 70% of the total soil microorganisms [1][2][3]. With as many as one billion bacterial cells and an estimated 8.3 million species in each gram of soil [4], the challenges in revealing and understanding the richness, diversity and composition of soil microbial community are thought-provoking. Nevertheless, the importance of the soil microbial community in nutrient acquisition, cycling, and availability is amply evidenced [5][6][7][8].
Numerous studies have demonstrated that healthy soil harbors abundant and diverse microbes, while environmental perturbation and management practices could lead to changes in microbial community structure that impact soil health and productivity [5][6][7][8][9]. Although long-term manure application resulted in enrichment of soil organic matter and microbial abundance [5,[10][11][12][13], it did not always translate into higher crop yield than soils supplemented with inorganic fertilizers [5,14]. Long-term inorganic fertilizer addition, on the other hand, may lead to soil acidification and decreases in microbial abundance and diversity [5,11,13,15]. Based on an analysis of 22 soils from a liming experiment spanning > 60 years, Rousk et al. [10] found that both relative abundance and diversity of bacteria were positively correlated with soil pH in the range from 4.0 to 8.3. In contrast, Lauber et al. [16] found that bacterial richness peaked in the near-neutral pH range based on analysis of 132,088 16S rRNA sequences from 88 soils. The observed inconsistency could, in part, be due to gradual system-level changes that make their detection challenging [7,8,17]. Although the dominant bacterial phyla in soil are reported to be Proteobacteria, Acidobacteria, Actinobacteria, Verrucomicrobia, and Bacteroidetes [18], their relative abundances varied considerably in pristine forest, grassland, winter wheat fields, and rice paddies [13,[19][20][21][22]. Of soil properties, carbon and N availabilities and pH are recognized as key factors governing changes in the soil bacterial community composition. However, the relative importance of these factors in driving changes in specific bacterial groups that impact ecosystem function and soil productivity is still not clear.
A century-long continuous winter wheat (Triticum aestivum L.) experiment provided an opportunity for extending our understanding of the effect of long-term organic and inorganic fertilizations on the bacterial community structure. We hypothesized that different fertilization practices would differentially affect soil properties and bacterial communities over century-long wheat cultivation. Insights gained from long-term studies are crucial to support sustainable wheat production. The specific objectives were to (1) evaluate richness, diversity, and composition of the soil bacterial community in response to century-long organic, inorganic, and no fertilization regimes under continuous winter wheat cultivation; and (2) identify drivers in these soils that govern the richness, diversity and composition of the bacterial community.

Site Description and Soil Sampling
This study was conducted in a century-long continuous winter wheat (Triticum aestivum L.) experimental field established in 1892 in central Oklahoma, USA (36 • 7 11" N, 97 • 5 19" W). The soil is a Kirkland silt loam (fine, mixed, thermic Udertic Paleustolls) with 37.5% sand, 40% silt, and 22.5% clay. Annual precipitation in this region (CD5 Central Oklahoma, OK, USA) ranged from 483 to 1448 mm (average 874 mm) and annual air temperature ranged from 14.4 to 17.7 • C (average 15.8 • C) based on data from 1895 to 2019 (https://climate.ok.gov). The field was conventionally tilled before planting every year. Comprehensive descriptions of this experimental site, soil properties, and winter wheat yield from 1892-2014 have been reported by Omara et al. [14], Girma et al. [23], and Aula et al. [24]. In this study, soil treatments were: no fertilizer added (Control), organic fertilizer applied in the form of cattle manure (Manure), and inorganic nitrogen, phosphorus, and potassium fertilizers (NPK). The manure treatment was initiated in 1899 and NPK treatment was added in 1929. Cattle manure from a feedlot was applied every four years at 269 kg N ha −1 and inorganic NPK fertilizers were applied annually in October before planting at 67 kg N ha −1 , 14.6 kg P ha −1 , and 28 kg K ha −1 in the forms of NH 4 NO 3 , Ca(H 2 PO 4 ) 2 , and KCl, respectively.
When the field experiment was initiated in 1892, the application of statistics to research was not yet firmly established. To compensate for the limited degrees of freedom from replication restrictions, we used an adapted sampling design to obtain three replicated samples from each fertilization treatment and repeated the study in two consecutive years (2012 and 2013). In each year, each treatment plot (30.5 × 6.1 m) was divided into three subplots (10.2 × 6.1 m) after wheat harvest. Five soil cores (0-15 cm) were collected from each subplot and composited as one soil sample. A total of 18 composite soil samples were collected for the three treatments in two years. Following sample collection, soils were placed in a cooler with dry ice and transported to the laboratory. Each field-moist sample was then passed through a 2 mm sieve, mixed thoroughly, and divided into two portions. One portion was air-dried and stored at room temperature for chemical analyses, and the other was freeze-dried and stored at −80 • C for subsequent biological and genetic analyses.
Soil pH was determined using a soil: water ratio of 1:2.5 (v/v) and a combination glass electrode (Mettler-Toledo International Inc., Columbus, OH, USA). Soil organic carbon (SOC) and total nitrogen (TN) were determined using a Carlo-Erba NA 1500 nitrogen/carbon/sulfur analyzer [25]. Soil NH 4 + and NO 3 − were extracted using 2 M KCl and analyzed with a BioTek Epoch 2 microplate reader (BioTek Instruments, Inc., Winooski, VT, USA) at OD660 and OD550, respectively. Microbial biomass carbon (MBC) was determined using the chloroform-fumigation incubation method with a K c factor of 0.45 [26,27].

Soil DNA Extraction, PCR Amplification and Sequencing
Soil DNA was extracted from 0.5 g soil using an Ultra Clean Soil DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA) according to the manufacturer's instructions. The quality and quantity of DNA was evaluated utilizing agarose gel electrophoresis and a NanoDrop-1000 spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, USA). All extracted DNA had an A 260 /A 280 > 1.8 and was of sufficient purity to be used in a polymerase chain reaction (PCR) following a 10-fold dilution with nuclease-free water.

Bioinformatics and Statistical Analyses
Prior to downstream analyses, de-noising and chimera checking were performed on the sequences using USEARCH [32] and UCHIME [33] to remove low quality sequences and error reads. Sequences with less than 250 bp, sequences with quality scores less than 25, and all chimeric sequences were removed [32][33][34][35]. The remaining 16S rRNA sequences were clustered into operational taxonomic units (OTUs) at 97% similarity using UPARSE [36]. Clusters with <2 members (singleton clusters) were removed.
The filtered high-quality sequences were analyzed to reveal the richness, diversity, and composition of soil bacterial communities under different fertilization treatments. Microbial α-diversity analyses were conducted using MOTHUR 1.44.1 and the Schloss standard procedure (http://www.mothur.org/) with modifications [37,38]. The 16S rRNA RDP reference database (Trainset9_032012.rdp) [39] was used to support analyses in MOTHUR. The evaluated α-diversity indices included two richness indices, Chao 1 [40] and abundance-based coverage estimator (ACE) [41], and two diversity indices, Shannon [42] and reverse Simpson (Invsimpson) [43]. Taxonomic assignment was conducted using the USEARCH global alignment algorithm against a database of high-quality 16S rRNA sequences derived from NCBI [32]. Three-way Venn diagrams (Venny 2.0, http://bioinfogp.cnb.csic.es/tools/venny/index.html) were used to reveal taxonomic groups that were shared by two or three treatments and taxonomic groups that were specific to one of the treatments.
One-way multivariate analysis of variance (MANOVA) was conducted in R 3.6.3 (https://www. r-project.org/; [44]) to determine whether soil properties, number of OTUs, α-diversity indices, and the relative abundance of bacterial phyla differ among treatments. Mean separation was carried out according to the Tukey's honestly significant difference post hoc test (HSD, p ≤ 0.05). Pearson's correlations between soil properties and the relative abundance of each bacterial phylum and their significances (p) were assessed with the "corr" function at p ≤ 0.05.
The variations in bacterial community structure among different treatments (β-diversity) were evaluated by non-metric multi-dimensional scaling (NMDS) analysis using the "vegan" package in R. Bray-Curtis distance matrix and the "ordiellipse" function were used to generate the NMDS plot. Permutational multivariate analysis of variance (PERMANOVA) was performed to test effects of different fertilizers on β-diversity using the "adonis" function (999 permutations) at a p value of 0.05. The multi-response permutation procedure (MRPP) was used to compare the within-groups dissimilarities (WD) and between-groups dissimilarities (BD) at p ≤ 0.05. Pearson's correlation coefficients (r) between the NMDS scores and tested soil properties were calculated with the "cor" function. Significant treatment differences (p ≤ 0.05) were tested with the post hoc permutation tests with 999 permutations. Vectors of Pearson's correlations between the NMDS scores and tested soil properties were added to the NMDS plot using the "envfit" function.

Bacterial Richness and Diversity
A total of 42,087 high-quality 16S rRNA gene sequences were used for the microbial community structure analyses. Bacterial rarefaction curves of the observed OTUs numbers approached plateaus in all treatments (Figure 1), indicating that the number of 16S rRNA sequences obtained in each treatment was sufficient to reveal detectable genotypes in the community. Bacterial richness and diversity were evaluated based on the number of OTUs observed in each soil, richness indices, diversity indices, and taxonomic groups specific to each treatment (Tables 1 and 2, Figures 1 and 2). The sequences were clustered into 1463, 1580, and 834 OTUs at 97% similarity in control, manure, and NPK treatments, respectively (Table 1). Bacterial richness and diversity indices, including Chao 1, Shannon, and Invsimpson, were significantly higher in manure and control treatments, and lower in NPK treatment ( Table 1). The number of bacterial groups' at all taxonomic levels from phylum to species were consistently and significantly lower in NPK than control or manure treatments. Manure had 1-1.3 times the bacterial groups compared with the control, while NPK had 72-84% of the number of bacterial groups detected in the control ( Table 2).    Control, century-long continuous wheat cultivation without fertilization; Manure, century-long continuous wheat cultivation received organic fertilizer cattle manure every four years at 269 kg N ha -1 ; NPK, century-long continuous wheat cultivation received inorganic fertilizer NPK annually at 67 kg N ha -1 , 14.6 kg P ha -1 , and 28 kg K ha -1 .
The number of bacterial groups that were shared by two or three treatments and the number of bacterial groups that were specific to one treatment are illustrated in the three-way Venn diagrams ( Figure 2). Higher grouping resolutions yielded higher percentages of bacterial groups that were specific to one treatment. At the phylum level, bacterial groups specific to control and manure treatments accounted for 7.7% and 10.3% of total phyla, respectively; bacterial groups specific to the NPK treatment were not detected at this taxonomic level. At the species level, bacterial groups specific to control, manure, and NPK treatments comprised 16.1%, 31.4%, and 10.4% of the total species, respectively. About 64% of the detected phyla were shared among the three treatments; while only 19% of the species were shared among the three soils. Of the soils tested, the number of bacterial groups specific to one treatment was generally highest in manure, followed by control, and lowest in NPK. Of the 399 species detected, 125 were specific to manure, 64 were specific to the control, while 42 species were unique to the NPK treatment. Compared to manure, NPK had a greater degree of impact on soil bacterial richness and diversity. Phyla or classes specific to NPK were not detected. At Figure 1. Rarefaction curves of 16S rRNA sequences obtained from soils tested. Curves were constructed based on the number of OTUs observed at 97% similarity per 100 sequences analyzed. Control, century-long continuous wheat cultivation without fertilization; Manure, century-long continuous wheat cultivation received organic fertilizer cattle manure every four years at 269 kg N ha −1 ; NPK, century-long continuous wheat cultivation received inorganic fertilizer NPK annually at 67 kg N ha −1 , 14.6 kg P ha −1 , and 28 kg K ha −1 .
The number of bacterial groups that were shared by two or three treatments and the number of bacterial groups that were specific to one treatment are illustrated in the three-way Venn diagrams ( Figure 2). Higher grouping resolutions yielded higher percentages of bacterial groups that were specific to one treatment. At the phylum level, bacterial groups specific to control and manure treatments accounted for 7.7% and 10.3% of total phyla, respectively; bacterial groups specific to the NPK treatment were not detected at this taxonomic level. At the species level, bacterial groups specific to control, manure, and NPK treatments comprised 16.1%, 31.4%, and 10.4% of the total species, respectively. About 64% of the detected phyla were shared among the three treatments; while only 19% of the species were shared among the three soils. Of the soils tested, the number of bacterial groups specific to one treatment was generally highest in manure, followed by control, and lowest in NPK. Of the 399 species detected, 125 were specific to manure, 64 were specific to the control, while 42 species were unique to the NPK treatment. Compared to manure, NPK had a greater degree of impact on soil bacterial richness and diversity. Phyla or classes specific to NPK were not detected. At the species level, about 10.4% were specific to NPK, which was significantly lower than the 16.1% and 31.4% specific to control and manure, respectively.

Bacterial Community Composition
The non-metric multi-dimensional scaling (NMDS) analysis based on the relative abundance of bacterial groups indicated significant (p ≤ 0.01) differences in the composition of soil bacterial communities under different fertilization regimes (Figure 3). The result was further supported by the significantly (p ≤ 0.01) larger between-groups distance (dissimilarity) than within-group distance based on the MRPP analysis ( Figure 3). To further evaluate the effects of different fertilization regimes on the composition of the soil bacterial community, we compared the relative abundance of bacterial phyla among the treatments (Figure 4). Of the 14 phyla detected in one or more soils in both years, Acidobacteria and Proteobacteria were most dominant, averaging 45.8% and 22.0%, respectively. Additional phyla included Actinobacteria (8.6%), Chloroflexi (6.2%), Firmicutes (4.9%), Bacteroidetes (4.4%), Verrucomicrobia (2.7%), Nitrospirae (2.4%), Gemmatimonadetes (1.3%), Cyanobacteria (0.94%), and TM7 (0.40%). Although the relative abundance of Firmicutes was not significantly different among soil treatments, the relative abundance of some bacterial groups varied considerably. For example, the relative abundance of Actinobacteria was highest in control, followed by manure, and lowest in NPK. The relative abundance of Chloroflexi, on the other hand, was highest in NPK, followed by control, and lowest in manure. The abundance of Cyanobacteria, a diazotrophic group, was significantly higher in the unfertilized control (2.59%), compared with manure (0.19%) or NPK (0.05%) treatments. Focusing on Cyanobacteria, Microcoleus paludosus and Leptolyngbya appalachiana were the two most dominant species, with the control population reaching 39 and 55 times of those in manure and NPK treatments, respectively ( Figure 5). Compared to the control, manure application

Bacterial Community Composition
The non-metric multi-dimensional scaling (NMDS) analysis based on the relative abundance of bacterial groups indicated significant (p ≤ 0.01) differences in the composition of soil bacterial communities under different fertilization regimes (Figure 3). The result was further supported by the significantly (p ≤ 0.01) larger between-groups distance (dissimilarity) than within-group distance based on the MRPP analysis ( Figure 3). To further evaluate the effects of different fertilization regimes on the composition of the soil bacterial community, we compared the relative abundance of bacterial phyla among the treatments (Figure 4). Of the 14 phyla detected in one or more soils in both years, Acidobacteria and Proteobacteria were most dominant, averaging 45.8% and 22.0%, respectively. Additional phyla included Actinobacteria (8.6%), Chloroflexi (6.2%), Firmicutes (4.9%), Bacteroidetes (4.4%), Verrucomicrobia (2.7%), Nitrospirae (2.4%), Gemmatimonadetes (1.3%), Cyanobacteria (0.94%), and TM7 (0.40%). Although the relative abundance of Firmicutes was not significantly different among soil treatments, the relative abundance of some bacterial groups varied considerably. For example, the relative abundance of Actinobacteria was highest in control, followed by manure, and lowest in NPK. The relative abundance of Chloroflexi, on the other hand, was highest in NPK, followed by control, and lowest in manure. The abundance of Cyanobacteria, a diazotrophic group, was significantly higher in the unfertilized control (2.59%), compared with manure (0.19%) or NPK (0.05%) treatments. Focusing on Cyanobacteria, Microcoleus paludosus and Leptolyngbya appalachiana were the two most dominant species, with the control population reaching 39 and 55 times of those in manure and NPK treatments, respectively ( Figure 5). Compared to the control, manure application generally led to enrichment of Proteobacteria, Nitrospirae, Bacteroidetes, and Gemmatimonadetes, and declines in Acidobacteria, Chloroflexi, Verrucomicrobia, and Cyanobacteria; NPK application resulted in higher relative abundances of Acidobacteria, Chloroflexi, Nitrospirae, Bacteroidetes, and TM7, but lower relative abundances of Actinobacteria, Verrucomicrobia, Cyanobacteria, and Gemmatimonadetes.
generally led to enrichment of Proteobacteria, Nitrospirae, Bacteroidetes, and Gemmatimonadetes, and declines in Acidobacteria, Chloroflexi, Verrucomicrobia, and Cyanobacteria; NPK application resulted in higher relative abundances of Acidobacteria, Chloroflexi, Nitrospirae, Bacteroidetes, and TM7, but lower relative abundances of Actinobacteria, Verrucomicrobia, Cyanobacteria, and Gemmatimonadetes. Figure 3. Non-metric multidimensional scaling (NMDS) ordination of soil bacterial community structure in different fertilization regimes. Control, century-long continuous wheat cultivation without fertilization; Manure, century-long continuous wheat cultivation received organic fertilizer cattle manure every four years at 269 kg N ha -1 ; NPK, century-long continuous wheat cultivation received inorganic fertilizer NPK annually at 67 kg N ha -1 , 14.6 kg P ha -1 , and 28 kg K ha -1 . The degree of separation between groups ranges from −1 to 1, where 0 means not different, 1 means completely different. Ellipses represent the 95% confidence interval for the standard error of the distances of samples in each microbial community. With increasing distance between two ellipses the community structure becomes more dissimilar; * p ≤ 0.05 means that there was a significant difference in the structure between two microbial communities based on the permutational multivariate analysis of variance (PERMANOVA). WD stands for within-groups distance and BD stands for between-groups distance; * p ≤ 0.05 means that there was a significant difference between WD and BD based on the MRPP analysis. Vectors show Spearman's correlations of variation in bacterial community composition and soil properties tested. Soil pH and NO3 − had significantly strong correlations to the first axis (NMDS1) or second axis (NMDS2) of NMDS ordination. . Non-metric multidimensional scaling (NMDS) ordination of soil bacterial community structure in different fertilization regimes. Control, century-long continuous wheat cultivation without fertilization; Manure, century-long continuous wheat cultivation received organic fertilizer cattle manure every four years at 269 kg N ha −1 ; NPK, century-long continuous wheat cultivation received inorganic fertilizer NPK annually at 67 kg N ha −1 , 14.6 kg P ha −1 , and 28 kg K ha −1 . The degree of separation between groups ranges from −1 to 1, where 0 means not different, 1 means completely different. Ellipses represent the 95% confidence interval for the standard error of the distances of samples in each microbial community. With increasing distance between two ellipses the community structure becomes more dissimilar; ** p ≤ 0.01 means that there was a significant difference in the structure between two microbial communities based on the permutational multivariate analysis of variance (PERMANOVA). WD stands for within-groups distance and BD stands for between-groups distance; * p ≤ 0.05 means that there was a significant difference between WD and BD based on the MRPP analysis. Vectors show Spearman's correlations of variation in bacterial community composition and soil properties tested. Soil pH and NO 3 − had significantly strong correlations to the first axis (NMDS1) or second axis (NMDS2) of NMDS ordination.  Deinococcus-Thermus, Fibrobacteres and Thermodesulfobacteria were less dominant bacterial phyla detected in control and manure treatments, but not in NPK. Deinococcus-Thermus populations were near the detection limit in the control but averaged 0.02% of the community in manure-treated soils.   Deinococcus-Thermus, Fibrobacteres and Thermodesulfobacteria were less dominant bacterial phyla detected in control and manure treatments, but not in NPK. Deinococcus-Thermus populations were near the detection limit in the control but averaged 0.02% of the community in manure-treated soils. Deinococcus-Thermus, Fibrobacteres and Thermodesulfobacteria were less dominant bacterial phyla detected in control and manure treatments, but not in NPK. Deinococcus-Thermus populations were near the detection limit in the control but averaged 0.02% of the community in manure-treated soils. Application of manure also promoted Fibrobacteres, averaging 0.22% of the community compared to 0.04% in the control. The impact of manure treatment on Thermodesulfobacteria was limited, and these bacteria were not detected in NPK-treated soils.

Soil Parameters and Their Relationships with Bacterial Community Structure
Century-long research plots with different fertilization regimes had significantly different soil pH, ranging from 4.6 (NPK) to 6.6 (manure) ( Table 3). The trends observed for soil organic C (SOC), total N (TN), and microbial biomass C (MBC) were similar, being significantly higher in the two fertilized soils compared to the unfertilized control. Differences between manure and NPK treatments were not significant. Focusing on the plant-available forms of N, NH 4 + was significantly higher in manure than NPK, but NO 3 − was significantly higher in NPK compared to the manure treatment. NMDS analysis and Pearson's correlation further showed the interrelationships between soil properties and the bacterial community composition (Table 4, Figure 2). Of the evaluated soil properties, soil bacterial community structure was significantly correlated to pH (r = 0.786, p ≤ 0.05) and NO 3 − (r = 0.933, p ≤ 0.01). Of the 20 dominant phyla detected, pH had significant positive correlations to Proteobacteria (r = 0.703, p ≤ 0.05) and Gemmatimonadetes (r = 0.700, p ≤ 0.05), and significant negative correlations to Acidobacteria (r = −0.744, p ≤ 0.05) and Chloroflexi (r = −0.766, p ≤ 0.05). Similarly, NO 3 − had a significant positive correlation to Acidobacteria (r = 0.750, p ≤ 0.05), but significant negative correlations to Actinobacteria (r = −0.773, p ≤ 0.05), Verrucomicrobia (r = −0.924, p ≤ 0.0001), and Cyanobacteria (r = −0.911, p ≤ 0.0001) ( Table 5). Note: Pearson's correlations between the NMDS scores and tested soil properties were assessed with the "cor" function and "envfit" function and post hoc permutation tests (999 permutations). Correlations with p ≤ 0.05 * and p ≤ 0.01 ** were considered as significant. NMDS1 is the first axis of the NMDS plot, NMDS2 is the second axis of the NMDS plot. SOC, soil organic carbon; TN, total nitrogen; MBC, soil microbial biomass carbon.

Soil Parameters and Their Relationships with Bacterial Community Structure
All soils used in this study were of the same soil series under long-term continuous winter wheat cultivation. The detected differences in the richness, diversity, and structure of bacterial community would likely reflect impact induced by different long-term soil fertilization regimes, evidenced by the consistent trends in results from samples in both years. Clear separation of the tested treatments in NMDS plot suggested significant effect of soil fertilization regimes on the bacterial community structure (Figure 3). The richness and diversity of the bacterial community in manure treatment were higher or similar to those of control, and both were significantly higher than NPK treatment, which was evidenced by the number of OTUs observed, richness indices, diversity indices, and the number of bacterial groups at each taxonomic level (Tables 1 and 2; Figures 1 and 2). Compared to manure, NPK had a greater impact on soil bacterial community structure, which was evidenced by the lower percentage of bacterial groups shared between NPK and control than the percentage of bacterial groups shared between manure and control at all taxonomic levels (Figure 2), as well as the larger distance (dissimilarity) in microbial community structure between NPK and control than the distance between manure and control in the NMDS plot ( Figure 3).
Promotion of bacterial richness and diversity by long-term application of manure could be due to multifaceted changes in soil, including enrichment of organic C and reduction of acidity. Previous studies using these soils showed significantly higher microbial biomass and enzyme activities in manure treatment than control or NPK treatment [5,7]. Application of manure led to a 16-20% increase in microbial activity [22,45], a marked increase in the abundance of gram-negative bacteria [21], and proliferation of selected microbial groups and changes in community composition [46,47].
A loss of richness and diversity by long-term application of NPK was clearly demonstrated in this study. It is well established that long-term application of inorganic fertilizer often leads to soil acidification [13,15,48]. Addition of available N, P, and K could shift the microbial community in favor of those that are less competitive in acquiring these nutrients but have higher tolerance to acidity. Therefore, the NPK treatment may have exerted selective pressures that led to reduced diversity and development of a less evenly distributed bacterial community. The relatively low number of bacterial groups and diversity indices observed among the treatments (Tables 1 and 2; Figure 2) were consistent with shifts in the microbial communities.
Century-long wheat production without fertilization led to marked depletion of soil nutrition and productivity [5]. The average grain yield in the control between 1890 and 2015 was about 52% and 46% of manure and NPK treatments, respectively [14,23]. Nevertheless, control continuously produced an average of >1100 kg ha −1 of wheat grain every year for over a century [14]. Aside from limited nutrients added by rainfall, microbes with the capability to fix atmospheric N 2 may have played an important role in maintaining soil productivity.

Soil Bacterial Groups under Century-Long Fertilization Regimes and Their Ecological Implications
Nutrient availability, especially C and N, may markedly impact the predominant microbial groups found in soils [49,50]. High nutrient availability promoted fast-growing copiotrophs while nutrient-limited soils favored slow-growing oligotrophs [49,51]. In the present study, Acidobacteria and Proteobacteria were the two most dominant phyla in soils under all three fertilization regimes, comprising 43-52% and 18-27% of the bacterial community within each soil, respectively (Table 5). This finding is consistent with a number of reports showing that Acidobacteria and Proteobacteria were the two most dominant bacterial groups, with their relative abundance varying considerably across different soil environments [3,13,18,52]. The abundance ratio of Proteobacteria:Acidobacteria can serve as a broad indicator of trophic status across a range of terrestrial soils [53]. In this study, the ratio was highest in manure (0.64), followed by control (0.47), and lowest in NPK (0.35) treatment ( Figure 4).
Of the three fertilization treatments, manure generally led to enrichment of all phyla except Cyanobacteria, NPK resulted in shrinking of most phyla but thriving of Chloroflexi, the control led to promotion of Cyanobacteria but detraction of Nitrospirae (Table 5). The lowest relative abundance of Chloroflexi was found in the manure treatment, while the highest abundance was in NPK (Figure 4). The results suggested that Chloroflexi might be less competitive in a diverse community but thrive in disturbed, high-stress environments, consistent with previous reports [19,20,54]. Suleiman et al. [20] reported that Chloroflexi, the most dominant rare (with ≤ 1% abundance) phylum detected, constituted a higher abundance in deforested grassland (0.6%) than pristine forest (0.4%). In agroecosystems, its abundance was higher, especially in soils under long-term cultivation. Sheng et al. [19] reported that Chloroflexi was significant higher in a century-long rice paddy (15.1%) than soils associated with ≤30 years of vegetable cultivation (1.1-4.3%). In the wheat field of this study, Chloroflexi accounted for 1.7-10.4% of the bacterial community, and was the third dominant phyla. Assuming pristine forest is a non-disturbed system, deforested grassland and agricultural fields would be systems with higher disturbance and stress. The higher occurrence of Chloroflexi in stressed environments has been reported previously, including permafrost at subzero temperatures [54], CO 2 -rich hypoxic soil [55], and PCB-contaminated soils [56]. Thus, the thriving of Chloroflexi in NPK could be related to acid stress resilience.
In the century-long unfertilized control plot, microbes with the capability to fix atmospheric N 2 may have played an important role in maintaining soil nutrient availability and wheat productivity. Cyanobacteria were most abundant in the control, with relative abundance values 13.6-51.8 times those in the two fertilized soils (Table 2). Another phylum that was promoted in the unfertilized control was Actinobacteria. Interestingly, Actinobacteria also contains strains capable of nitrogen fixation [57]. The dominance of Actinobacteria in soil environments is often reported [58,59], which was also the third most dominant groups in the tested soils. Actinobacteria contribute to biological buffering of soils and organic matter decomposition, as well as producing many bioactive metabolites (antibacterials, antifungals, and growth promoting substances for plants and animals) [60,61].

Soil Parameters Drive Soil Bacterial Community Structure
Soil properties play key roles in shaping microbial community structure [50]. In this study, soil pH and NO 3 − content were the two most-dominant factors influencing soil bacterial community structure under different fertilization regimes (Tables 4 and 5; Figure 3). The results were in agreement with previous studies [12,13,50] indicating nutrient availability (especially C and N) and soil pH are main factors for shifts in soil microbial communities under different fertilization regimes. Soil acidity, as a crucial factor governing bacterial richness and diversity, could serve as a selective pressure to promote or suppress the growth of specific bacterial groups. In this study, soil pH was lowest in NPK, followed by control, and highest in manure (Table 3), following the trends in bacterial richness and diversity in the three treatments (Tables 1 and 2; Figures 1 and 2). Furthermore, a significant correlation between pH and NMDS scores suggested an association with bacterial community structure. Fierer and Jackson [62] found that bacterial richness and diversity differed by ecosystem type, and these differences could largely be explained by soil pH (richness r 2 = 0.70, p ≤ 0.0001; diversity r 2 = 0.58, p ≤ 0.0001). Rousk et al. [10] also reported positive relationships between bacterial richness and soil pH (r 2 = 0.75, p ≤ 0.001), with the number of OTUs doubling between pH 4.0 and 8.3. Fierer and Jackson [62] and Zhang et al. [48] reported that most bacteria require a pH of 6.0-8.0, which is suitable for most proteins to function [63]. Of the bacterial phyla detected, pH had significant positive correlations with Proteobacteria and Gemmatimonadetes but significant negative correlations with Acidobacteria and Chloroflexi ( Table 5), suggesting that the variations in bacterial richness and diversity in three treatments could be primarily due to the promotion or suppression of microorganisms within these four phyla. With increasing pH from 4.0 to 7.0, Zhang et al. [48] also reported that Acidobacteria and Chloroflexi decreased, while Proteobacteria and Gemmatimonadetes increased.
Nitrogen availability was another determining factor regulating the richness, diversity, and composition of the soil bacterial community. Cyanobacteria was one of the phyla most affected by N content, evidenced by the significant negative correlations between the relative abundance of Cyanobacteria and soil total N (r = −0.867, p ≤ 0.01), NH 4 + (r = −0.836, p ≤ 0.01), and NO 3 − (r = −0.911, p ≤ 0.01) contents (Table 5). Ramirez et al. [17] also reported that the abundance of Cyanobacteria in soil was negatively correlated with N levels (r= −0.50 to −0.64, p ≤ 0.05) during 0-800 kg N ha −1 yr −1 .

Conclusions
Century-long organic or inorganic fertilization significantly altered the richness, diversity, and composition of the soil bacterial community compared to the unfertilized control. Bacterial richness and diversity were enhanced by manure addition but reduced by NPK application. Different fertilization regimes did not change the types of dominant phyla, but did affect the relative abundance of bacterial phyla in a community. Of the bacterial phyla detected, Acidobacteria and Proteobacteria were the most abundant phyla in all three communities. Of the three fertilization treatments, manure generally led to enrichment of all phyla except Cyanobacteria, NPK resulted in shrinking of most phyla but thriving of Chloroflexi, while century-long wheat cultivation without fertilization promoted the diazotrophic group Cyanobacteria. Soil pH and NO 3 − availability were the two most dominant parameters governing the richness, diversity, and composition of the soil bacterial community. Soil pH significantly correlated with Acidobacteria, Proteobacteria, Chloroflexi, and Gemmatimonadetes and was major contributors to variations in bacterial richness and diversity among the three treatments. Soil N contents were significantly correlated to Cyanobacteria, whose relative abundance in control was 13.6-51.8 times than those in the fertilized soils. Cyanobacteria, especially Microcoleus paludosus and Leptolyngbya appalachiana, could play an important role in maintaining wheat productivity in the century-long unfertilized control. This study extends our understanding of soil bacterial community structure under long-term organic, inorganic, or no fertilization regimes and provides insight for making improved fertilization strategies.