BNF plays an important role in N supply during rice growth
The rice yields of three fertilization treatments (NPK, NPKM, NPKS) were significantly higher than CK, while the yields of CK kept consistent since 2015 (Fig. 1a). Non-symbiotic BNF supplied 37.32% crop N for rice [26]. Thus, we inferred that nitrogen-fixing bacteria of CK played an important role to maintain rice production. The 15N isotope dilution method was used to evaluate the contribution of endophytic diazotrophs in the nitrogen nutrition supply for rice growth (Table 1). The root length, plant height and fresh weight of rice of C20 and T20 had no significant difference but were significantly lower than T10. The N concentration among the three treatments had no significant difference, but δ15N of C20 (33.26%±1.10%) and T20 (33.85%±0.25%) were significantly higher than T10 (29.43%±0.49%). It certified that endophytic nitrogen-fixing bacteria could supply N for rice by BNF under low nitrogen nutrition (%Ndfa=11.51).
Twenty-five isolates isolated from surface-sterilized rice roots including Proteobacteria (80%), Firmicutes (12%) and Actinobacteria (8%) (Additional file 2: Table S2). The nifH gene fragments were amplified from five isolates, then their nitrogenase activities were measured by ARA (Additional file 3: Table S3). The nitrogenase activity of strain BV6 (Burkholderia) was highest (25.71 µmol/h/ml).
The diversity of total rhizosphere and endophytic diazotrophs were shaped by growth period, year and fertilization
Illumina sequencing was performed to investigate the community compositions of rhizosphere and endophytic diazotrophs and their response to different fertilization managements in different growth periods. 318 samples collected from the field (120 rhizosphere soil DNA samples, 119 root DNA samples, 48 rhizosphere soil RNA samples and 31 root RNA samples) were sequenced, and after normalization, each sample had 22552 sequences. These sequences were clustered into 706 operational taxonomic units (OTUs) at 90% identity.
For alpha diversity analyses, richness index (Chao1), evenness index (Shannon) and phylogenetic diversity index (PD) of rhizosphere diazotrophs were significantly higher than endophytic diazotrophs at the DNA level (Additional file 4: Table S4), which was consistent with the results of rarefaction curves (Additional file 5: Figure S1). Interestingly, total rhizosphere diazotrophs had the highest α-diversity at the tillering stage, while total endophytic diazotrophs had the lowest α-diversity at this stage (Additional file 4: Table S4). The effects of growth period on all tested α-diversity indices were statistically significant, but fertilization only had significant effects on the evenness index (Table 2). For soil and root DNA samples, almost all α-diversity indices of diazotrophs among four fertilization treatments (CK, NPK, NPKM and NPKS) had no significant difference (Additional file 6: Table S5, Additional file 7: Table S6).
The phylogenetic trees of 50 dominant OTUs based on the sequences of nifH genes were generated (Fig. 2, Additional file 8: Figure S2). The relative abundances of dominant OTUs of total rhizosphere diazotrophs were notably different from total endophytic diazotrophs. Moreover, the whole community structure of total rhizosphere diazotrophs was notably different from and more stable than total endophytic diazotrophs (Additional file 9: Table S7, Fig. 1b). We visualized and quantified the differences between diazotrophic communities (β-diversity) using non-metric multidimensional scaling (NMDS) and permutational multivariate analysis of variance (PERMANOVA). The growth period was a primary driver of the total rhizosphere (R2=0.22, P<0.001) and total endophytic diazotrophic β-diversity (R2=0.33, P<0.001) (Fig. 1c-d, Additional file 10: Figure S3a-b, Table 2), revealing a stronger temporal effect on total endophytic diazotrophs than total rhizosphere diazotrophs. However, total rhizosphere diazotrophs were more sensitive to the year and fertilization management (R2=0.19 and 0.11, P<0.001) than total endophytic diazotrophs (R2=0.06 and 0.05, P<0.001). Year, growth period and fertilization management showed significant interactions in pairs.
Taxonomic composition and fertilizer-sensitive taxa
Fourteen phyla were determined by pooling the nifH sequences, and Proteobacteria was the predominant phylum, contributed approximately 57.0% and 76.1% of soil and root samples, respectively (Additional file 11: Figure S4). The top three classes of soil samples were Alphaproteobacteria (25.29%), Deltaproteobacteria (20.65%) and Betaproteobacteria (8.83%) (Fig. 1e, Additional file 12: Figure S5a). However, Betaproteobacteria (34.95%), Gammaproteobacteria (18.51%) and Alphaproteobacteria (15.37%) were the three dominant classes of root samples. At the genus level, Bradyrhizobium (24.80%), Geobacter (7.72%) and Burkholderia (5.71%) were the top three genera of soil samples (Additional file 12: Figure S5b-c), but for root samples, they were Burkholderia (26.74%), Thiorhodospira (8.88%) and Bradyrhizobium (8.25%).
In these four cluster trees (Fig. 1e, Additional file 12: Figure S5), soil DNA samples and root DNA samples were divided into two clusters, and samples from the same growth period were frequently clustered together except for soil samples in 2016. However, PERMANOVA and ANOSIM suggested that the growth period significantly affected the total rhizosphere diazotrophic communities in 2016 (data not shown). Generally, the communities of either total rhizosphere or endophytic diazotrophs of NPK and NPKS treatments were clustered together, and they were separated from CK and NPKM treatments (Fig. 1e, Additional file 12: Figure S5a). Furthermore, the soil samples of CK were quite different from the other three fertilization treatments, and differences between CK and NPKM were smaller than CK and NPK or CK and NPKS (Additional file 13: Table S8). We noted that the relative abundances of Pelobacter of NPKM were lower than that of CK and higher than that of NPK, but Azoarcus was the opposite (Additional file 14: Figure S6). Furthermore, for Bradyrhizobium, Geobacter, Heliobacterium, Rubrivivax and Desulfomonile, their abundance of of NPKM were between that of CK and NPK at most sampling points.
The rhizosphere soil physicochemical indices OM, TN, TK, TK, AN, AP, AK, pH, NO3-N and NH4-N were tested (Additional file 15: Table S9). Soil pH value of CK was highest among four fertilization treatments. NPKM increased the content of organic matter (OM), total N (TN), total P (TP), available N (AN) and available P (AP) and NO3-N. Fertilization management, year and period had significant effects on all tested physicochemical parameters (P<0.01, Additional file 16: Table S10), except that the period had no significant effects on TN and NH4-N (P>0.05). Most rhizosphere and endophytic diazotrophic α-diversity indices had negative correlation with TP and AN at six sampling points (Additional file 17: Table S11). The Mantel test showed that AK, TK and NH4-N were significantly negatively correlated with the community structures of total rhizosphere and endophytic diazotrophs (Additional file 18: Table S12). Correlation analysis showed that pH had a significant positive correlation with Pelobacter, but had a significant negative correlation with Thermodesulfovibrio in rhizosphere soil (Additional file 19: Figure S7).
The community structures of total endophytic diazotrophs were more stabilized in the reproductive stage than the vegetative stage
The similarity distances (1-Bray dissimilarity) of diazotrophic community structures among three rice growth periods were calculated at each fertilization treatment (Fig. 3a-b). For root samples, the similarity distances between the tillering and heading stages (red nodes) were higher than the tillering and mature stages (green nodes), but lower than the heading and mature stages (blue nodes). However, this trend was not observed in soil samples.
Tracking diazotrophic community changes at the class level revealed that many classes showed distinct temporal dynamics in root samples (Fig. 3c-d). For example, Alphaproteobacteria and Deltaproteobacteria enriched as time went on, but Betaproteobacteria decreased. The endophytic diazotrophic abundance of “others” (mainly unclassified taxa) remarkably increased, and was close to rhizosphere diazotrophs at the mature stage. The abundance of shared OTUs between soil and root samples in soil samples (Pink flows) apparently increased during rice cultivation, but slightly decreased in root samples.
The linear discriminate analysis (LDA) effect size (LEfSe) was used to identify the distinguishing diazotrophs during different growth periods in paddy soil and rice roots [27]. The numbers of the discriminating taxa of rhizosphere diazotrophs, whose LDA scores between different growth periods greater than 2, were 172 and 96 in 2015 and 2016, and of the endophytic diazotrophs were 149 and 97 (data not shown). The discriminating genus significantly increasing at the same sampling period over two years were filtered and heatmaps were drawn to demonstrate more particular information of the responses of these genera to the growth period (Fig. 3e-f). For the rhizosphere soil samples, the relative abundances of Treponema, Chloroherpeton, Paludibacter, Coraliomargarita and Desulfobulbus were significantly higher at the tillering stage than at the heading and mature stage (Fig. 3e), and they were correlated positively with AK (Additional file 19: Figure S7a-b); the relative abundances of Oscillatoria, Cyanothece and Coleofasciculus were significantly higher at the heading stage, and Oscillatoria was correlated negatively with NO3-N; the relative abundances of Heliobacterium, Sinorhizobium and Frankia were significantly higher at the mature stage, and they were correlated positively with NO3-N. For root samples, the relative abundance of Tolumonas was significantly higher at the tillering stage than at the other stages (Fig. 3f); the relative abundances of Thiorhodospira, Sulfurospirillum, Pelosinus and Selenomonas were significantly higher at the heading stage, and Thiorhodospira was correlated negatively with pH and NO3-N; the relative abundance of Rhizobium, Desulfovibrio, Geobacter, Frankia and Acetobacterium were significantly higher at the mature stage, and they were correlated negatively with AK (Additional file 19: Figure S7c-d).
Diversity and community structures of active nitrogen-fixing bacteria
Consistent with the total root-associated diazotrophs, the community structures of active rhizosphere diazotrophs were notably different from active endophytic diazotrophs (Additional file 9: Table S7); growth period was a primary driver of the active rhizosphere (R2=0.15, P<0.001) and active endophytic diazotrophic β-diversity (R2=0.17, P<0.001); fertilization management had a higher effect on the β-diversity of active rhizosphere diazotrophic (R2=0.14, P<0.001) than endophytic diazotrophs (R2=0.10, P>0.05) (Additional file 20: Table S13, Additional file 10: Figure S3c-d). The rhizosphere diazotrophic community structures of CK and NPKM had no significant difference, and both of them had significant differences with the rhizosphere diazotrophic community structures of NPK and NPKS. (Additional file 13: Table S8). The relative abundances of Geobacter of NPKM were higher than that of CK and lower than that of NPK in soil RNA samples (Additional file 14: Figure S6). Furthermore, TP, and AN, and had negative correlation with α-diversity indices of active rhizosphere and endophytic diazotrophs, and NH4-N and NO3-N had significant negative correlation with active rhizosphere diazotrophs (Additional file 17: Table S11). The Mantel test showed that TK, AN, AP, pH and NO3-N were significantly correlated with the community structures of active rhizosphere diazotrophs, and NO3-N had a significant correlation with the community structures of active endophytic diazotrophs (Additional file 18: Table S12).
Richness, Evenness and PD indices of all DNA samples were higher than corresponding RNA samples (Additional file 4: Table S4). The β-diversity of the total and active diazotrophic communities was investigated as well (Additional file 10: Figure S3e). The soil DNA samples clustered together, whereas the soil RNA samples, root DNA samples and root RNA samples were more dispersed. To compare their aggregation difference, the similarities (Z scores) in between-group diazotrophic community composition were visualized in bar plots. Z scores of the soil RNA samples, root DNA samples and root RNA samples had not significant differences but they were significantly lower than the soil DNA samples (Fig. 4a).
The rhizosphere DNA and RNA samples separated on the cluster tree (Fig. 4c), but root DNA samples and their corresponding RNA samples frequently clustered together (Fig. 4d). The correlation coefficient of the total and active diazotrophic communities in soil (Spearman correlation: r= 0.244, P<0.05) was slightly smaller than that in roots at the OTU level (Spearman correlation: r= 0.271, P<0.05). In Fig. 4b, Z scores of the soil DNA and RNA samples at the heading stage were highest, but there was no significant differences between Z scores of the soil DNA and RNA samples at tillering and mature stages, and Z scores of the root DNA and RNA samples at heading and mature stages.
The abundance changes of diazotrophs at the RNA and DNA levels were calculated in the rhizosphere soil and rice roots (Additional file 21: Figure S8a-b). Burkholderia, Anaeromyxobacter, Bradyrhizobium and Pelobacter predominated in the total and active rhizosphere diazotrophic communities. NifH gene expressions in Burkholderia and Anaeromyxobacter were somewhat higher than that of the total rhizosphere diazotrophic community, whereas nifH gene expressions in Bradyrhizobium and Pelobacter were opposite. Burkholderia, Azoarcus, Thiorhodospira and Bradyrhizobium predominated in the total and active endophytic diazotrophic communities. Growth period and fertilization appreciably affected the nifH gene expressions of endophytic diazotrophs. For instance, the abundance changes of Azoarcus at the RNA and DNA levels at the heading and mature stage were the opposite. The abundance changes of Burkholderia, Thiorhodospira and Bradyrhizobium at the RNA and DNA levels of CK and the other three fertilization managements were the opposite. Due to high between-sample variations, which might be caused by environmental differences, there was a few or no statistically significant difference of diazotrophic abundance at the RNA and DNA level in the soil or roots.
Root-enriched nitrogen-fixing bacteria
To investigate the distribution of nitrogen-fixing bacteria in rice roots and rhizosphere soil, the relative abundances of 39 dominant genera in these two niches were compared by the STAMP software. The size of the bubbles in the chart (Fig. 5) represents the relative abundance of 39 dominant genera in rice root samples, and the color of the bubbles represents the difference between the square root of the relative abundances of these genera in the roots and soil.
At the DNA and RNA levels, the relative abundances of Thiorhodospira, Rhizobium and Burkholderia in all root samples at three growth periods were almost significantly higher than that in corresponding soil samples (P<0.05); the relative abundances of Paludibacter, Magnetospirillum, Clostridium and Rhodospirillum in root samples at the heading and mature stages were significantly higher than that in corresponding soil samples; the relative abundances of Sulfurospirillum at the heading stage, Azoarcus and Frankia at the mature stage in root samples were significantly higher than that in corresponding soil samples. However, Bradyrhizobium, Pelobacter, Heliobacterium, Geobacter and Anaeromyxobacter had significantly higher relative abundance in the rhizosphere soil than that in rice roots.
Co-occurrence patterns of total and active root-associated diazotrophs
Co-occurrence networks were built to explore potential interactions and niche-sharing of diazotrophs in paddy soil and rice roots during three key growth periods (Fig. 6, Additional file 22: Figure S9). In each of these networks, nodes represent OTUs (relative abundance>0.01%) and edges represent significant co-occurrence relationships (Spearman |r|>0.65 and P<0.01). Topological properties of co-occurring network were calculated (Additional file 23: Table. S14). Overall, the ecological networks of rhizosphere and endophytic diazotrophs had different network sizes and degree of connectivity, and showed distinct time related successions at three key rice growth periods. Consistent with the α-diversity analyses (Additional file 4: Table. S4), networks of rhizosphere diazotrophs were more complex in comparison with endophytic diazotrophs. For instance, the numbers of nods (237.5±51.64) and edges (1444±840.69) in soil DNA samples were greater than the numbers of nods (117.67±48.55) and edges (201.83±136.28) in root DNA samples. Moreover, soil RNA samples had greater numbers of nods (224.67±43.04) than root RNA samples (150±36.77), and had greater numbers of edges (704.67±306.95) than root RNA samples (485±417.19).
Compared with networks at the tillering stage, the rhizosphere and endophytic diazotrophic assemblages formed larger and more complex networks at the heading and mature stages. Networks of total rhizosphere and endophytic diazotrophs at the heading and mature periods exhibited similar network topological structures (nodes, links, average path distance, average clustering coefficient, average degree and modularity) in 2016. However, for active diazotrophs, soil RNA network at the mature stage was most complex, and root RNA network at the heading stage was more intricate. The percent of positive correlations in total rhizosphere diazotrophic networks (57.31%±3.63%) was much less than that in total endophytic diazotrophic networks (83.05%±7.52%), active rhizosphere diazotrophic networks (88.81%±9.22%) total endophytic diazotrophic networks (88.43%±15.64%).