Nitrogen has a greater influence than phosphorus on the diazotrophic community in two successive crop seasons in Northeast China

Fertilizer-induced changes in soil nutrients regulate nitrogen (N) fixation in the terrestrial biosphere, but the influences of N and phosphorus (P) fertilization on the diazotroph communities in successive crop seasons were unclear. In this study, we assessed the effects of N and P (high vs. low doses) on the abundance and structure of N2-fixation communities after wheat and soybean harvest in a long-term (34 and 35 years) fertilization experiment. In both seasons, long-term N addition significantly decreased the abundance of nifH genes and 16S rDNA; in addition, high doses of N and P fertilizer decreased the richness of diazotrophs, whereas low doses did not. The proportion of the dominant genus, Bradyrhizobium, in the soybean season (86.0%) was higher than that in the wheat season (47.9%). Fertilization decreased diazotroph diversity and the relative abundance of Bradyrhizobium in the wheat season, but had insignificant effects in the soybean season. The addition of N, but not P, significantly changed the communities of both diazotrophs (at the genus level) and rhizobia (at the species level) in the two seasons. Soil pH was positively associated with nifH abundance and diazotrophic richness; soil NO3− content was negatively correlated with diazotrophic richness and positively correlated with diversity. Soil pH and NO3− content were the two main drivers shaping the soil diazotrophic community. Overall, long-term inorganic N had a greater influence than P on both diazotrophic abundance and community composition, and diazotrophic diversity was more clearly affected by fertilization in the wheat season than in the soybean season.


Relationship between soil chemical properties and α-diversity and gene copies. A t-test com-
parison across all treatments in both crop seasons revealed a significant decline in soil pH (Table 1). In contrast, the concentrations of NO 3 − , Avail P, Avail K, TN, and OM, which were clearly higher in soil with fertilizer treatments (N 1 , N 1 P 1 , N 2 , and N 2 P 2 ) than in unfertilized soil (Table 1). Ammonium (NH 4 + ) was greater under theN and P fertilizer treatments than CK, whereas no statistically significant differences were observed for plots fertilized during the soybean season, regardless of fertilization (Table 1).
The abundance of the dominant species, Bradyrhizobium diazoefficiens, was negatively correlated with soil NO 3 − (r = − 0.563, P < 0.05) and NH 4 + (r = − 0.764, P < 0.01) contents, and positively associated with soil pH (r = − 0.636, P < 0.05) (Fig. 4A). In addition, the relative abundance of this species was significantly and negatively correlated with the N-level in the wheat season, whereas no significant relationship was observed in the soybean season (Fig. 4D).
The clustering results showed that the nifH communities were significantly separated by N addition treatments in the wheat ( NMDS at the OTU level. The NMDS results revealed that nifH community composition varied significantly (P< 0.05) with respect to N addition but not P addition. Across the two seasons, the phylogenetic structure of the nifH communities shifted in similar ways. The nifH communities in soils in the wheat season did not differ significantly from those in the soybean season ( Figure S2). Three separate groups were clearly observed along NMDS1 (accounting for 67.31% of the variation in nifH community): no fertilized soils (SCK and WCK, orange circles); low N and low N plus P fertilized soils (SN 1 , SN 1 P 1 , WN 1 , and WN 1 P 1 , green circles); and high N and high N plus P fertilized soils (SN 2 , SN 2 P 2 , WN 2 , and WN 2 P 2 , red circles, except for SN 2 and SN 2 P 2 ) (Supplementary Fig. S2). As the amount of added N increased, communities became more different from those in unfertilized soils in both seasons.
Environmental effects on diazotrophs. In the MRT, the dominant lineages were first split by NO 3 − content, which explained 26.0% of the variation in community structure (Fig. 5A). The tree explained 71.1% of the variance in the standardized diversity indices. At the second node, the split was determined by soil pH, which explained 17.8% of the variation. The communities were then split by pH and NO 3 − , accounting for 14.3% and 8.4% of the variation in the data, respectively (Fig. 5A). The results of RDA with Monte Carlo permutation tests showed that NO 3 − , pH, Avail P, and Avail K were significantly (P < 0.05) correlated with the changes in the composition of the N-fixing community, with contributions of 41.4%, 17.1%, 13.2%, and 9.2%, respectively (Supplementary Table S5). All samples were separated into two groups (wheat and soybean season: squares and circles, respectively) ( Fig. 5B) along the RDA2 axis, except for SCK-3. Along the RDA1 axis, nifH communities under low N, and low N plus low P were clearly different from those under CK, and those under high N, and high N plus high P were more different in the two crop seasons. We found that climate factors (temperature and precipitation) were not significantly (P > 0.05 in MRT, P = 0.05 in RDA) correlated with the changes in the composition of the N-fixing community

Discussion
Fertilization decreased bacterial and diazotrophic abundance. The correlation between nifH abundance and fertilizer inputs was controversial in previous studies. In our study, the abundance of nifH in fertilizer treatments decreased in the two seasons, in agreement with conclusions reached by Zhang 24 and Coelho 25 validating Fan's 11 hypothesis that N fixation and fixers would become less abundant over time in fertilized environments. Urea addition in the present study may have been particularly detrimental for obligate N fixers; their ability to downregulate fixation was limited, and thus they exhibited relatively narrow growth tolerance. For the first time, we reported that the diazotroph to bacteria ratio was more sensitive during the growth and development of soybean than that of wheat. The application of N fertilizers is generally expected to decrease the dependence of the ecosystem on free-living N-fixers 26 , and this finding might explain the slight decrease in nifH/16S rDNA in the wheat season. Although symbioses between some N fixing bacteria with soybean can provide N-fixers with an exclusive niche and contribute to their growth 27 , we found a clear decrease in nifH/16S rDNA in the soybean season, thus indicating that the symbiotic ability of N-fixers with the soybean decreased under a 35 year N fertilization regime. We concluded that the clear differences in this ratio in the two seasons could be explained by the interactions among individual plant strategies, in agreement with Sheffer's 28 conclusion. NifH abundance under N 1 and N 2 did not show clear differences from those under N 1 P 1 and N 2 P 2 , respectively, thus indicating that P did not have a significant effect on nifH gene copy number.
Most studies have found that an increase in N decreases soil pH 29,30 , and we found the same result; however, P had no such effect. In this study, we found that with increased acidification of black soil, the growth and reproduction of N-fixers was significantly inhibited. Other studies have found that the copy number of the nifH gene is strongly positively associated with soil NH 4 +31 and available K 32 , but negatively correlated with total N 31 .
Effects of fertilization on diazotrophic diversity. Many studies have found that N addition increases the richness of diazotrophs 33,34 . However, some studies have found that N has no such clear influence 35 . However, we observed significantly lower Chao index values for N fixers under the addition of high N, and high N plus high P in the two seasons ( Fig. 2A). One possible reason for this finding is that N fixers have a strong advantage in N deficient conditions, but some species have difficulty surviving and may even die under intense increases in the concentration of available N in the soil microenvironment 35 .
Coelho's 36 study on another nonlegume crop, sorghum, showed that the diazotrophic Shannon diversity under high levels of N was lower than that under low levels of N, in agreement with the results in the wheat season in our study. However, the higher proportion of Bradyrhizobiaceae in the soybean season than in the wheat season might explain the lower diversity of nifH sequences in the soybean season. Furthermore, nodules formed by Bradyrhizobium and soybean roots increase the tolerance to various stresses, such as salt 37 , acidity 38 , drought 39 , insecticide 40 , and high aluminium 38 , thus potentially also explaining why nifH diversity was not significantly affected by N in the soybean season but was significant in the wheat season. We additionally found that P had clear inhibitory effects on the richness and diversity of diazotrophs.
Shannon indices of diazotrophic and nitrite-dependent anaerobic methane oxidation bacteria 41 were positively correlated with NO 3 − content, thus indicating that the increase in NO 3 − was beneficial to the diversity of N cycling microorganisms. Santoscaton 42 found that NO 3 − loads are associated with bacterial 16S rDNA abundance but not nifH gene abundance, similar to our results. N fertilizer affects the structural composition of N-fixing bacteria. N but not P fertilizer had significant effects on diazotrophic community composition, which indicated that the level of N fertilizer was the most important factor affecting the structural composition of N-fixing bacteria in the black soil of Northeast China. This result was highly consistent with the response of diazotrophic bacterial 25 , ammonia-oxidizing archaeal 43 , bacterial 44 and fungal 5 communities to N fertilization regimes.
The process of nitrification in soil is performed partly by gram-negative bacteria in the family Bradyrhizobiaceae, in a process involving the conversion of NH 4 + into NO 2 and subsequently NO 3 -45 . Therefore, the higher average concentration of NH 4 + in soybean (42.08 mg kg −1 ) than in wheat (37.06 mg kg −1 ) soils, may lead to an increase in Bradyrhizobiaceae in the soybean season. The high abundance of Bradyrhizobiaceae in the soybean season could also be explained by stable symbiosis between leguminous plants (soybean) and rhizobia, although the roots of nonleguminous plants (wheat) can be colonized by rhizobia 25 . Furthermore, linear relationships between the cultivar and the bacterial community have been reported, such as genotype associations of maize with Azospirillum 46 , alfalfa cultivars with Sinorhizobium 47 , and sorghum cultivars with Paenibacillus 48 . Therefore, the results presented here emphasize the importance of cultivar type in selecting N-fixing strains for use as wheat and soybean inoculants. N fertilizer affects the structural composition of rhizobia. Rhizobia, a collective name for the symbiotic N-fixing bacteria associated with legumes, comprise 14 genera 1 , six of which were found in the current study. The community structure of rhizobia was distinguished by N levels. Bradyrhizobium was reported to be more adapted to acidic soils 49 , while we found a lower abundance under N 2 P 2 with a lower pH (4.79). Thus, we propose that Bradyrhizobium may use suitable amounts of available N to support their growth, whereas N fixation and N fixers will become increasingly less important when NO 3 is excessive. This hypothesis is based on the negative correlation between Bradyrhizobium and the soil NO 3 − content (Fig. 4A). Ahmed 50 concluded that soil NO 3 − has a negative effect on the activity of N-fixing rhizobia by inhibiting the function of the enzymes nitrogenase and leghaemoglobin.
The dominance of B. diazoefficiens over Bradyrhizobium sp., B. japonicum, and B. elkanii revealed a unique community structure of soybean rhizobia in the black soil, a finding not consistent with those of Yan 51  www.nature.com/scientificreports/ negative correlation between Bradyrhizobium diazoefficiens and N fertilizer in the wheat season rather than in the soybean season may be explained by the sensitivity of certain bacterial species present in plant types to N fertilizer. In addition, the relatively higher content of NO 3 − under N 2 and N 2 P 2 in the wheat season may have caused Bradyrhizobium diazoefficiens to become increasingly less important. The abundance of Bradyrhizobium diazoefficiens in the wheat season indicated that it is a genospecies whose growth is clearly inhibited by N fertilizer.

Effects of soil properties on the diazotrophic community. Researchers have confirmed that soil
physicochemical characteristics affect the activity of N-fixers 52 . MRT and RDA results confirmed that the soil NO 3 − content was the most important contributor to the soil diazotroph community, a finding consistent with reports by Yang 53 and Zou 54 . Moreover, soil NO 3 − content was identified as an important predictor of 16S rDNA gene abundance and the α-diversity of the diazotroph community ( Supplementary Fig. S1). Neutral or slightly acidic soil conditions are conducive to biological N fixation 18 , while the strong acidity in high N and P may be a severe problem for N fixation, because in such environments legume nodules fail to form, and some rhizobia become inactive 18 . Seminal work by Wang 16 highlighted the importance of soil pH as a fundamental driver of the distribution of the diazotrophic community, and we reached the same conclusion. In this study, soil acidification in black soil in northeast China, caused by high levels of N fertilization, usually leads to problematic nutrient deficiency or mineral toxicity during N fixation 55 . These findings may aid in predicting the response and feedback of the diazotroph community in farmland ecosystems to high levels of N fertilization.
More recently, researchers have shown that diazotroph diversity and richness are mainly influenced by soil available P 53 and available K 16 ; our results again validated these conclusions. Our findings indicated that soil nutrient availability, which was highly responsive to fertilizer input, was crucial for the establishment of the soil diazotrophic community structure 56 in black soil in Northeast China.

Conclusion
Our work provided solid evidence, after 34 and 35 years of experiments, that N fertilization largely influenced diazotroph communities in the soil in two successive crop seasons in northeast China. N is likely to have a greater influence than P on diazotrophic bacteria. The community structure of N-fixing bacteria and rhizobia was clearly associated with the level of N fertilizer. The lower diazotrophic abundance under N fertilizer treatments may have diminished the capacity for biological N fixation in the two seasons. N had greater effects on diazotrophic diversity and the relative abundance of the dominant genus Bradyrhizobium in the wheat season than in the soybean season. The different response patterns of diazotrophic abundance, community composition, and diversity to the soil properties revealed a complicated mechanism underlying the diazotrophic population's adaptation to long-term N and P fertilization in two crop seasons. However, we only conducted research at the DNA level, and future research will determine the impact of N fertilization on the functional diversity of diazotrophs in two seasons based on mRNA profiling of nifH genes.

Materials and methods
Experimental design and sample collection. The experimental site was set up in 1979 in Harbin city, Heilongjiang Province, China (45° 40ʹ N, 126° 35ʹ E and altitude 151 m), which is in the secondary terrace of the Songhua River. The soil type is black soil, and the parent material is flooded loess-like clay. The tillage method is a combination of shallow tillage and deep rotation. The annual crop rotation of wheat, soybean and maize was repeated every 3 years in the field with five fertilization treatments in a completely randomized block design with three replicates: CK (without fertilizer), N 1 (low N), N 2 (high N), N 1 P 1 (low N plus low P) and N 2 P 2 (high N plus high P). Taking into account the difference in N needs for wheat and soybean, we chose different amounts of low N (75 and 150 kg urea ha −1 year −1 for wheat and soybean, respectively) and high N (150 and 300 kg urea ha −1 year −1 for wheat and soybean, respectively) treatments. The detailed types and amounts of fertilizer were shown in Supplementary Table S1. N, P, and K fertilizers are all applied after harvest of the previous crop in autumn (September). We collected soil samples in September 2013 and 2014 after wheat and soybean harvests, respectively before the next fertilization. The annual average soil temperature, 10 cm below the surface of the soil was 7.4 and 5.9 °C, and the annual precipitation was 2262 and 296 mm in 2013 and 2014, respectively 44 . The experimental design and sample collection are shown in Fig. 6 [The map was produced using 'R (i386 3.1.2, https ://www.R-proje ct.org, The R Core Team, 2019)' , with the open source packages 'maps' , 'mapdata' , and 'maptools ' 57 ]. Soil samples were randomly collected from the plow layer of soil (5-20 cm) and stored as described in our previous research 44 . Soil chemical properties. Soil pH was measured at a 1:5 ratio of soil to distilled water (weight/volume). Nitrate (NO 3 − -N) and ammonium (NH 4 + -N) were extracted from 5 g of air-dried soil with 2 M KCl, steam distillation and titration. The available P (Avail P) was extracted with 0.5 M NaHCO 3 and determined with the molybdenum blue method. Available potassium (Avail K) was extracted with 1 M ammonium acetate and determined by flame photometry. The organic matter (OM) and total N (TN) were determined according to Strickland and Sollins (1987).
High-throughput sequencing and bioinformatics analysis. Fast DNA SPIN Kit (MP Biomedicals, Santa Ana, CA, USA) was used to extract DNA from 0.5 g of fresh soil. The nifH gene was amplified using the primer pairs nifH f and nifH r 58 . Primer sets and PCR reactions were as detailed in Supplementary Table S2, and amplification reactions were sequenced on the Illumina MiSeq PE300 platform. The raw reads have been deposited in the National Center for Biotechnology Information Database (SRX 1034826). The nifH nucleotide sequences were analysed with the QIIME-1.9.1 pipeline. Briefly, the low quality sequences were discarded, and Quantitative PCR analysis. The abundance of bacterial 16S rDNA and the nifH gene were analysed with an ABI 7500 Real-Time PCR detection system with primers 515F-806R 3 and nifH f-nifH r. Primer sets, the qPCR amplification system and the thermal programme are detailed in Supplementary Table S2. Plasmid DNA containing 16S rDNA and nifH fragments were used for quantitative PCR standards. The specificity was determined by melting curve analysis and agarose gel electrophoresis 60 . The ratio of N-fixing microorganisms to bacteria was calculated according to the nifH gene and 16S rDNA copy numbers.
Statistical analysis. Analysis of variance was performed with a randomized complete block design in IBM SPSS Statistics 21. Linear regression analysis was performed to test for statistical significance and the strength of associations between soil chemical properties and α-diversity and gene copy numbers (16S rDNA and nifH) in Origin 2020. Phylogenetic tree was visualized using GraPhlAn2 61 with the data of OTU representive sequences and OTU abundant table. On the basis of Bray-Curtis similarity distance, nonmetric multidimensional scaling (NMDS) was used to analyse the nifH community structure at the OTU level. A multivariate regression tree (MRT) analysis was performed with the package "mvpart" in the "R" statistical program to interpret the main relationships between the biological data (at the class level) and environmental factors [soil chemical properties and climate status (temperature and precipitation)] 62 . The correlations between the N-fixing communities (at the genus level) and environmental factors were determined with redundancy analysis (RDA), by using CANOCO 5.0. A logarithmic transformation was performed to normalize the data and the significance (P-value) for the first two canonical axes was evaluated by means of Monte Carlo tests based on 999 permutations.  www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.