The diazotrophic community in oat rhizosphere: effects of legume intercropping and crop growth stage

In this study, the abundance, diversity and structure of the diazotrophic community in oat rhizosphere soil in three cropping systems and at two oat growth stages were investigated using real-time PCR and Illumina MiSeq sequencing. The nifH gene abundance in oat-soybean intercropping (OSO) and oat-mungbean intercropping (OMO) was significantly greater than that in sole oat (O), but the nifH gene abundance significantly decreased at the later stage in all the treatments. Alpha diversity indices in OSO and OMO were higher at the heading stage, but lower at the maturity stage than that in O. Bradyrhizobium and Skermanella were the dominant genera identified in all samples, with an average proportion of 35.8% and 12.4%, respectively. The proportion of dominant genera showed significant differences and varied with cropping system and growth stage. Principal component analysis showed that growth stage had a stronger effect than intercropping on the diazotrophic community structure. However, Mantel test and redundancy analysis showed there was no environmental factor significantly correlated to the diazotrophic community structure. Our results demonstrate that intercropping had a weaker effect than growth stage on the abundance, diversity and structure of the diazotrophic community in oat rhizosphere soil.


Introduction
Intercropping, a long-established farming technique of cultivating two or more crops in the same space at the same time, is an important component of agricultural production systems in most developing countries [1,2] . The advantages of intercropping include higher overall productivity, better utilization of land and resources, and the establishment of soil microbial diversity [3][4][5] . Cereal-legume intercropping is a commonly used cropping system in northern China [6] . It can benefit from the biological N 2 -fixing of legumes, which reduces the use of synthetic N fertilizer and offers an alternative for the development of sustainable agriculture [7] .
Biological nitrogen fixation (BNF), the enzymatic reduction of N 2 to ammonia, is important in the global nitrogen cycle [8] and supplies approximately 110 Tg$yr -1 N to the land ecosystems and 140 Tg$yr -1 N in the oceans [9] . The BNF process is catalyzed by a highly conserved nitrogenase enzyme which is made up of two multisubunit metalloproteins; the heterotetrameric core encoded by the nifD and nifK genes, and the dinitrogenase reductase subunit encoded by the nifH gene [10] . The nifH gene exists in almost all diazotrophs and a good phylogenetic correlation has been found between the nifH and 16S rRNA genes [10] , indicating the nifH gene could be considered as a suitable marker to investigate the diversity and composition of diazotrophic organisms [11,12] .
Soil diazotrophs are sensitive to agricultural practices such as cropping system, plowing and fertilization [13,14] . Intercropping increased the population and diversity of the soil bacterial community in the rhizosphere soil more than sole cropping [15,16] . Xiao et al. [17] reported that the abundance and diversity of nifH genes showed variations in soil with continuous and rotational soybean cropping. Sampling time also affected the diazotrophic structure in agricultural soil [18] . In short, both cropping system and sampling time had great impacts on the abundance and structure of soil diazotrophs. However, the abundance, diversity and structure of the diazotrophic community in the rhizosphere soil of cereal crop as affected by intercropping with legume and cereal crop growth stage in cereal-legume intercropping system has not been sufficiently investigated.
The aim of this study was to evaluate the effects of intercropping with legume and crop growth stage on the abundance, diversity and structure of the diazotrophic community in oat rhizosphere soil using real-time PCR and Illumina MiSeq sequencing.

Field sites and soil sampling
A field experiment was conducted in Baicheng Academy of Agricultural Sciences (45°37′ N, 122°48′ E), Baicheng City, Jilin Province, China. The area has a temperate semihumid climate with an average annual temperature and precipitation of 4.9°C and 407.9 mm, respectively. The experiment had three treatments: sole oat (O), oat-soybean intercropping (OSO) and oat-mungbean intercropping (OMO), with three replicates for each treatment in a random plot design (6 m Â 6 m). O was sown at row distance of 30 cm, OSO consisted of three oat rows at 30 cm with two soybean rows at 60 cm with 45 cm between the oat and soybean rows, and OMO consisted of three oat rows at 30 cm with two mungbean rows at 50 cm, with 40 cm between the oat and mungbean rows. P and K fertilizers were applied at 55 kg$hm -2 P 2 O 5 as (NH 4 ) 2 PO 4 and 45 kg$hm -2 K 2 O as K 2 SO 4 . No other N fertilizer was applied. Oat, soybean and mungbean were sown on May 9, and harvested on August 7, August 13 and September 7, respectively.
Rhizosphere soil samples were collected at the heading (July 6) and maturity (August 7) stages of oat in 2013, with samples coded as O_H, OSO_H, OMO_H, O_M, OSO_M and OMO_M, respectively. For each plot, five groups of oat plants (5 cm length) were removed from the soil by hand and the soil adhering to the root hairs after gentle shaking was sampled as rhizosphere soil. The rhizosphere soil samples from each replicate were mixed together to give a composite sample, which was then divided into two parts. One part was stored at 4°C for soil properties analysis and the other part stored at -80°C for DNA extraction.

Soil chemical properties analysis
Soil pH was determined with a soil to water ratio of 1:2.5, using a pH-meter (FE28, Mettler, Toledo, OH, USA). Total organic carbon and total nitrogen were determined using the K 2 Cr 2 O 7 oxidation-reduction titration method and Kjeldahl digestion method [19] , respectively. Soil nitrate and ammonium were extracted from fresh soil with 2 mol$L -1 KCl at a solution to soil ratio of 2.5:1 and determined by a SAN ++ Continuous Flow Analyzer (Skalar Analytical B.V., Breda, The Netherland).

DNA extraction and real-time PCR
DNA was extracted from 0.3 g of fresh soil using the E.Z. N.A. ® Soil DNA Kit (Omega Bio-tek Inc., Norcross, GA, USA) according to manufacturer's instructions. DNA from each sample was extracted in three replicates and pooled to form one mixed DNA sample. The extracted DNA was checked on a 1% agarose gel and the concentration was determined using a NANO Quant (Tecan, Männedorf, Switzerland).
Real-time PCR was performed on an ABI 7300 thermocycler (Applied Biosystems, Waltham, CA, USA) using the SYBR ® Premix Ex Taq TM (Takara, Dalian, China) according to the manufacturer's instructions, with 1 mL DNA extract (1-10 ng) as the template in 25-mL reaction mixtures using the primers PolF and PolR [20] . Plasmid DNA contained the nifH gene was used as the template for the standard curve construction. Briefly, the nifH gene was cloned into the pMD18-T vector (Takara, Dalian, China) according to the manufacturer's instructions. The positive plasmid containing correct insert was extracted and the concentration of plasmid DNA was determined using a NANO Quant (Tecan, Männedorf, Switzerland) to calculate the nifH gene copy number. At last, 10-fold serial dilutions of a known copy number of the plasmid DNA were subjected to real-time PCR assay in triplicate for the generation of an external standard curve.
Amplicons were extracted from 2% agarose gels and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer's instructions and quantified using a QuantiFluor TM -ST (Promega, WI, USA). Purified amplicons were pooled in equimolar amounts and paired-end sequenced (2 Â 350) on an Illumina MiSeq platform according to the standard protocols.

Processing of sequencing data
Raw fastq files were de-multiplexed, quality-filtered using QIIME (version 1.17) with the following criteria: (1) The 350 bp reads were truncated at any site receiving an average quality score < 20 over a 10 bp sliding window, discarding the truncated reads that were shorter than 50 bp; (2) exact barcode matching, 2 nucleotide mismatch in primer matching, reads containing an ambiguous characters were removed; (3) only sequences that overlapped longer than 10 bp were assembled according to their overlap sequence. Reads that could not be assembled were discarded. Operational taxonomic units (OTUs) were clustered with 97% similarity cutoff using UPARSE (version 7.1) [22] and chimeric sequences were identified and removed using UCHIME [23] .
The phylogenetic affiliation of nifH gene was analyzed by ribosomal database project (RDP) Classifier (RDP website: RDP Release 11) [24] against the FGR functional gene database (FunGene website). The read numbers in all the samples were normalized to the same sequencing depth and ACE, Chao1 and Shannon indexes of community diversity were chosen to evaluate community diversity using mothur (version v.1.30.1, Mothur website: Alpha_ diversity) [25] .

Statistical analysis
The copy numbers of nifH gene were log-transformed prior to statistical analysis and mapped using Sigmaplot software (version 12.5) (Systat Software Inc. San Jose, CA, USA). The analysis of variance (ANOVA) and Pearson correlation analysis of soil properties and nifH gene abundance were performed using SPSS software (version 20.0) for Windows (IBM Corp., Armonk, NY, USA). Mantel test, principal component analysis (PCA) and redundancy analysis (RDA) were performed using the vegan package in the R Language and Environment for Statistical Computing [26] .

Rhizosphere soil chemical properties
Soil pH values ranged from 7.00 to 7.85, and significantly higher pH values were observed in the samples taken at the maturity stage than at the heading stage (P < 0.05, Table 1). Significant variations were observed in total N, NH 4 + -N and NO 3 --N concentration between different rhizosphere soil samples (P < 0.05). The highest total N content was observed in O_H (1.02 g$kg -1 ) and the lowest was detected in O_M (0.95 g$kg -1 ), and no significant difference was observed among the remaining rhizosphere soil samples. Significantly higher NH 4 + -N concentrations were observed in the rhizosphere soil samples at the heading stage than at the maturity stage (P < 0.05). The highest NO 3 --N concentration was detected in OSO_M (10.62 mg$kg -1 ) which is 6.32 times more than that in OMO_M (1.45 mg$kg -1 ). The C/N ratios in O_H and OSO_H were significantly lower than in the rest rhizosphere soil samples. However, no significant difference was observed in organic matter content (17.21-18.50 g$kg -1 ) between all the rhizosphere soil samples.

Abundance of nifH gene
The abundance of nifH gene ranged from 1.99 Â 10 8 to 3.28 Â 10 9 copies per gram d.w.s and significant differences were observed in the nifH gene abundance among different rhizosphere soil samples (P < 0.05, Fig. 1). The highest nifH gene abundance was detected in OMO_H and was 15.5 times more than that in O_M. The nifH gene abundances in O_H, OSO_H and OMO_H were 2.54, 2.70 and 7.96 times more than that in O_M, OSO_M and OMO_M, respectively. Furthermore, the nifH gene abundances in OSO and OMO were significantly higher than that in O at both stages (P < 0.05). Correlation analysis indicated the nifH gene abundance was negatively correlated with soil pH (r 2 = 0.820, P < 0.001). Additionally, cropping system (P < 0.001), crop growth stage (P < 0.001) and their interaction (P < 0.001) had significant impacts on the nifH gene abundance.

Sequencing results and diversity indices
Illumina MiSeq sequencing was used to investigate the diazotrophic diversity and composition. A total of 67959 qualified reads were obtained from all the rhizosphere soil samples and each library had 7168 reads with 99 to 139 OTUs obtained at the 97% similarity after normalization. All rarefaction curves tended to approach the saturation plateau (Fig. S1) and Good's coverage was more than 99.8% for each library, indicating the sequencing data were reliable and the nifH gene sequences identified represented Table 1 Chemical properties of the rhizosphere soil samples  Table 2). The alpha diversity evenness (ACE), species richness (Chao1) and Shannon index confirmed different trends among the three treatments at the two stages of oat ( Table 2). The diversity indices of ACE, Chao1 and Shannon index were higher in OSO_H and OMO_H than in O_H, but lower in OSO_M and OMO_M than in O_M. The ANOVA results showed significant differences with P < 0.0001 for Shannon index, but no significant difference was observed for Chao1 (P = 0.145) and ACE (P = 0.387) indices.

Taxonomy composition of nifH gene in different rhizosphere soil samples
Sequences that could not be classified into any known group were assigned as unclassified. Taxonomy groups that accounted for less than 1% in average proportion of total reads were combined into one part as "Others" in drafting. All the OUTs were assigned into 6 different phyla, 16 families and 19 genera, with 3 phyla, 10 families and 10 genera shared by all samples. Proteobacteria was the most abundant group which made up 92.6% to 98.5% of the total reads at the phylum level in different samples. Proteobacteria_unclassified, Bradyrhizobiaceae and Rhodospirillaceae were the three dominated groups which accounted for 20.3%-53.9%, 26.8%-46.3% and 4.81%-20.6%, respectively, of the total reads at the family level in different samples. At the genus level, each library had 11 to 15 genera, and Proteobacteria_ unclassified,  Bradyrhizobium and Skermanella were the three most abundant groups in all samples, with a relative abundance of 20.3%-53.9%, 26.8%-46.3% and 4.81%-20.4%, respectively (Fig. 2). Moreover, 98.0% and 52.6% of the total reads were assigned to the phylum and genus levels, respectively, while only 15.9% of the total reads could be classified to the species level (data not shown).

Correlation between the diazotrophic groups and rhizosphere soil properties
The correlation analysis results showed that significant correlations were observed in soil pH, organic matter, total nitrogen, NH 4 + -N and C/N ratio with the proportion of some dominant diazotrophic groups at the genus level ( Table 3). The relative abundance of Skermanella was negatively correlated with soil pH (r 2 = 0.825, P < 0.05), organic matter (r 2 = 0.848, P < 0.05) and C/N ratio (r 2 = 0.886, P < 0.05). The Proteobacteria_unclassified abundance was negatively correlated with soil pH (r 2 = 0.910, P < 0.01) and C/N ratio (r 2 = 0.832, P < 0.01). The Bradyrhizobium abundance was positively correlated with soil pH (r 2 = 0.840, P < 0.05), but negatively correlated with NH 4 + -N concentration (r 2 = 0.945, P < 0.01). The Azohydromonas abundance was negatively correlated with soil pH (r 2 = 0.876, P < 0.01), but positively correlated with C/N ratio (r 2 = 0.884, P < 0.05). However, no significant correlation was obtained between NO 3 --N concentration and all diazotrophic genus groups.
Mantel tests were performed to analyze the environmental factors responsible for the structural changes of the diazotrophic community in different rhizosphere soil samples. The results showed that no environmental factor was significantly correlated to the diazotrophic structure  (data not shown). Furthermore, redundancy analysis (RDA) was performed to identify the rhizosphere soil properties affected the diazotrophic structure (Fig. 3). Parameters that could explain the variance of community structure were ranked by RDA results (Table 4), which showed that no environmental factor was significantly correlated to the diazotrophic community structure.
3.6 Comparison of the diazotrophic structure among different rhizosphere soil samples Principal component analysis (PCA) was used to identify the diazotrophic structure differences in different rhizosphere soil samples (Fig. 4). The resolution of PCA at the genus level showed that O_H, OSO_H and OMO_H located to the left, whereas O_M, OSO_M and OMO_M were distributed to the right. Additionally, OSO_H and OMO_H were together and far from O_H, while OSO_M was close to O_M but far from OMO_M, indicating that the diazotrophic structures were differential at different stages and cropping systems. These results demonstrated that both cropping system and crop growth stage influenced the diazotrophic structure in oat rhizosphere soil samples, and cropping system had a weaker impact than crop growth stage on the diazotrophic structure.

Differences between the diazotrophic genera among different rhizosphere soil samples
At the genus level, the differences between the diazotrophic communities from different rhizosphere soil samples were analyzed via Venn diagrams (Fig. 5). A total of 19 genera were generated, 99.8% and 97.6% of which belonged, respectively, to the shared genera at the heading and maturity stages (Fig. S2). The OSO_H alone contained more diazotrophic variations (4 genera) than O_H alone (1 genus) and OMO_H alone (1 genus). Moreover, only one genus (Rhodomicrobium) was shared by only O_H and OSO_H and three genus groups (Trichormus, Nostocaceae_unclassified and Rhodospirilla-ceae_unclassified) were shared by only O_M and OSO_M, respectively.

Discussion
Nitrogen fixation is an important procession for global N cycle and crop production and is performed by diazotrophs possessing the nitrogenase enzyme. The nifH gene is the best-known gene for encoding the dinitrogenase reductase subunit [10] and is wildly used to investigate the diazotrophic community in the rhizosphere of legume [27] , tobacco [28] , rice [29,30] and wheat [31] . In this study, we combined real-time PCR and Illumina MiSeq sequencing to reveal the effects of cropping system and crop growth stage on the abundance, diversity and structure of the diazotrophic community in oat rhizosphere soil. The nifH gene abundance across the three cropping systems and two stages were found to be much higher than that observed in potato fields [18,32] , wheat fields [33] and sorghum rhizosphere soils [34] . The presence of symbiotic N-fixers in oat [31,35] and intercropped with leguminous crops might be the main reasons for the higher abundance of nifH gene in the soil samples assayed. Soil pH is an important factor influencing the nifH gene abundance and high pH value favors the potential for BNF [36] . The pH  value was negatively correlated with nifH gene abundance in this study, which is inconsistent with previous studies [18,32] . The pH values in our study ranged from 7.0 to 7.8, while pH values in the study of Pereira e Silva et al. [18,32] ranged from 4.3 to 7.7. Variations in soil acidity and alkalinity may thus explain the differences in nifH gene abundance. Additionally, organic carbon has been shown to have a positive effect on the nifH gene abundance [32] ; however, no significant correlation was observed between organic carbon and nifH gene abundance in this study. This inconsistency can be explained by the soil type differences in these studies. The soil studied previously was mostly sandy while soil in those studies were high in clay content and it was observed that an increase of clay content from 60% to 80% drastically decreased the abundance of soil diazotrophs [37] . In addition, the stronger effect of soil pH on nifH gene abundance may weaken the effect of organic carbon.
Significantly higher nifH gene abundance in the intercropping treatments was observed than in sole oat (P < 0.001) and significantly lower nifH gene abundance in samples at the maturity stage was observed than at the heading stage (P < 0.001), indicating both cropping system and crop growth stage had significant influences on the diazotrophic abundance. This was in line with previous studies that significantly higher abundance of soil microbes was observed in intercropping rhizosphere soil samples than that in sole soybean or maize [15,38] . Sampling time was shown to influence the soil diazotrophic abundance [18] , and Hai et al. [34] showed that sampling stage had a more significant effect than fertilization on the nifH gene abundance in sorghum rhizosphere soil. Root exudates have also been shown to affect soil diazotrophic populations [35,39] . In addition, cropping system, crop growth stage and their interaction had significant impacts on the nifH gene abundance. All these results indicated there was no single causative factor but soil pH could be a potentially important factor influencing the diazotrophic abundance in the soil samples assayed.
In accordance with previous observations [40,41] , the diazotrophs belonging to the Alphaproteobacteria, Betaproteobacteria and Gammaproteobacteria classes were observed in all samples and accounted for 92.6% to 98.5% in this study. Cyanobacteria, Firmicutes and Verrucomicrobia diazotrophs were also detected in some samples, at less than 3.2% on average. The high abundance of Alphaproteobacteria (50.6%) could be explained by the presence of legumes [40] . Although Rhizobium is well known to be detected in a broad range of non-legumes [42] , no Rhizobium was found even though Rhizobiales is the major order in Alphaproteobacteria. In addition, the high proportion of Bradyrhizobium is noteworthy, as it is known as a symbiotic N-fixer [18] and it may benefit the growth of oat. The differences in relative abundance of Rhizobium and Bradyrhizobium showed different adaptations of the two genera to the studied soil. The abundance of Azohydromonas, Bradyrhizobium and Skermanella showed differences between the three cropping system treatments, suggesting cropping system influenced the diazotrophic structure. The relative abundance of Azohydromonas, Bradyrhizobium and Skermanella increased at the later stage for oats, indicating crop growth stage also greatly influenced the diazotrophic structure. This is in agreement with earlier reports that crop growth stage (or sampling time) affected the abundance and structure of the diazotrophic community [18,32,43] . Moreover, the presence of inter-cropped soybean and/or mungbean and root exudates from soybean and/or mungbean could be another reason for the diazotrophic structure changes in this study.
Of the soil properties, only pH was significantly correlated with the nifH gene abundance. However, significant correlations were observed with soil pH, organic matter, NH 4 + -N and C/N ratio with some dominant diazotrophic genera, indicating these properties affected the diazotrophic structure. However, Mantel tests and RDA results showed no soil property was significantly correlated with the diazotrophic structure. These results suggested that the change in diazotrophic structure was likely to have been caused by multiple factors rather than a single soil property. Furthermore, the PCA profile implied that intercropping had a weaker effect on the diazotrophic structure than crop growth stage, which is in line with previous findings [33,44] .

Conclusions
In conclusion, our results demonstrate that both intercropping and crop growth stage can have significant impacts on the abundance and composition of the diazotrophic community in oat rhizosphere soil. Intercropping with legumes increased the abundance of nifH gene, but this abundance decreased from the heading to maturity stage of oat. Bradyrhizobium and Skermanella were the two dominant diazotrophic genera identified in oat rhizosphere soil samples. Both intercropping and crop growth stage changed the abundance and composition of the diazotrophic community and intercropping had a weaker effect than crop growth stage on the diazotrophic community.