The field experiment was established in 2013 at a long-term experimental station of the Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, in Changchun, Jilin Province, China (44°59′ N, 125°23′ E). The mean annual precipitation and air temperature were 614 mm and 6.4°C, respectively. Three tillage practices, NT, RT, and MP tillage, were conducted in the same location, which were arranged in a completely randomized block design with four replicates, and the area for each replicate plot was 25 × 7.8 m. Soybean is cultivated in rotation with corn in alternate years, and this cropping rotation was applied under the three tillage practices. After harvest, the soybean residues were fully and directly returned to the soil surface of the three tillage practices. The maize stalks were manually cut into approximately 30 cm in length and laid on the soil surface of each tillage practice for their corresponding crop residues. NT experienced no soil disturbance except for planting with a two-row John Deere 7200 NT planter (Moline, IL, United States). RT included ridging in June, and no other soil disturbance was conducted until the next year. In contrast to NT and RT, MP treatment involved two tillage applications after planting, with one fall moldboard plowing (approximately 20 cm in depth) after harvest and one spring disking (approximately 7.5–10 cm in depth) with ridge building before planting. For the cultivation of soybean under NT, RT, and MP, 40 kg N ha^–1, 60 kg P ha^–1, and 80 kg K ha^–1 were applied as base fertilizers. For the cultivation of corn under the three tillage practices, 100 kg ha^–1 nitrogen (N), 45.5 kg P ha^–1, and 78 kg K ha^–1 were applied as base fertilizers, and another 50 kg N ha^–1 was used as top dressing at the jointing stage of corn vegetative growth. The soybean yield was determined by hand in 3-m lengths of six interior rows from each plot.

Each tillage treatment was arranged in a randomized block design with four replicate plots. Two repeated soil samples were collected from each plot, and a total of eight repeated soil samples were obtained from each tillage treatment. Specifically, nine individual soil cores (0–20 cm soil depth) were randomly collected and composited together as a repeated soil sample to minimize within-plot variation. Overall, 24 bulk and 24 rhizosphere soil samples (three tillage practices × eight replicates) were collected under soybean cultivation at the beginning pod stage of soybean reproductive growth (on 25 July 2017). Each soil sample was placed in an individual sterile plastic bag and immediately transported back to the laboratory. After passing through a 2-mm sieve, the soil sample was divided into two parts: one was stored at −80°C for DNA extraction and the other was stored at 4°C for determination of soil properties. Soil physical and chemical properties, including soil pH, total carbon (TC), total nitrogen (TN), total phosphorus (TP), total potassium (TK), available phosphorus (AP), available potassium (AK), NH₄⁺-N, and NO₃^–-N, were measured as previously described in Hu et al. (2017).

The total DNA was extracted from 0.5 g of each soil sample (stored at -80°C) using a FastDNA^® SPIN Kit for Soil (MP Biomedicals, United States) according to the manufacturer′s procedures, and the DNA concentrations were measured with a NanoDrop 2000 spectrophotometer (Thermo Scientific, United States). The extracted DNA was stored at −20°C until downstream analysis. The primers PolF (5′-TGC GAY CCS AAR GCB GAC TC-3′) and PolR (5′-ATS GCC ATC ATY TCR CCG GA-3′) were used to amplify the nifH gene (Poly et al., 2001). A sample-specific barcode (6 bp) was added to the forward primer to distinguish the amplified products, which were sequenced on an Illumina MiSeq PE 300 platform at Majorbio BioPharm Technology Co., Ltd., Shanghai, China. The sequences obtained from this research were deposited at NCBI with accession number SRP 226893.

The raw nifH nucleotide sequences were analyzed through Quantitative Insight into Microbial Ecology (QIIME) pipeline¹ (Caporaso et al., 2010). A total of 964,583 high-quality sequences of nifH gene were acquired after removing the barcodes, primers, and low-quality sequences with an average quality score < 20. The obtained sequences were further translated to amino acid sequences on FunGene pipeline² (Fish et al., 2013), and the amino acid sequences that did not match the nifH protein were discarded. The remaining sequences were clustered into operational taxonomic units (OTUs) at 95% similarity level with UPARSE (Edgar, 2013), and the taxonomy of each OTU was aligned against the nifH gene database³. A total of 10,288 OTUs were clustered from 964,583 high-quality nifH gene sequences (from 9,587 to 23,604) across 47 soil samples. One replicate of RT in bulk soil (RTB1) failed in the PCR production procedure, and thus it was discarded in the downstream analysis. To minimize the impact of sequence count variation among samples, a subset of 9,587 sequences per sample, based on minimum sequences, was randomly extracted prior to the downstream analysis. The relative abundance of each diazotrophic group at the phylum, genus, and species taxonomic levels was used to compare the difference and/or similarity of diazotrophic communities in response to tillage practices for the subsequent analysis.

Diazotrophic ecological networks were constructed with sequencing data through Molecular Ecological Network Analysis Pipeline⁴ based on random matrix theory methods, which were described previously (Zhou et al., 2010, 2011; Deng et al., 2012). In brief, eight replicates for each tillage practice for bulk or rhizosphere soils were uploaded into the pipeline, and only OTUs detected in more than half of the replicates were considered. In total, six diazotrophs’ networks associated with three tillage practices for bulk and rhizosphere soils were obtained with automatically generated cutoff values (similarity threshold, S_t). Random networks were generated by keeping the same number of nodes and links as the above empirical networks and were used to examine the statistical significance of the network properties using the Maslov–Sneppen method (Maslov and Sneppen, 2002). The modularity of all the obtained networks was over 0.4, which suggested that they had modular structures (Newman, 2006). A node in a network represents the OTU, and an edge or link between the pair of OTUs was assigned if the correlation between their abundance exceeds the S_t. The “global network properties,” “individual node centrality,” and “module separation and modularity” were calculated on the pipeline. Among these network properties, node degree or connectivity is the number of neighbors of a node, clustering coefficient represents how well a node is connected with its neighbors, density represents how densely the network is populated with edges, and a module is a group of nodes that are highly linked within the group but with few links outside the group (Zhou et al., 2011; Deng et al., 2012). In addition, the percentage of positive or negative edges in each network was calculated to examine the cooperative or the competitive interactions between diazotrophic species. Finally, based on this pipeline, we incorporated soil properties into microbial networks to measure the effects of soil properties under the three tillage practices on diazotrophs’ network structures.

The α-diversity of a diazotrophic community was calculated using the alpha_diversity.py function in QIIME. Significant differences in soil properties, diazotrophic taxa, and α-diversity under different tillage practices were tested by analysis of variance (ANOVA). Pearson’s correlations between the α-diversity or the abundance of diazotrophic taxa and soil properties or network complexity were determined using SPSS ver. 22.0 software. Principal coordinate analysis (PCoA) was conducted based on weighted UniFrac distances to depict the patterns of the β-diversity for diazotrophic communities in R software (ver. 3.5.1), using the “ape” package. The statistical significance of diazotrophic community structures was tested with Adonis analysis in R using “vegan” package. Mantel test was used to link the structure of a diazotrophic community with the soil properties and was conducted in R using “vegan” package. Cytoscape 3.7.2 software was used to visualize the network graphs (Shannon et al., 2003).

Except for NH₄⁺-N, the soil properties were significantly different among NT, RT, and MP in bulk and rhizosphere soils (Table 1). Compared with NT and RT, MP significantly increased soil pH in both bulk and rhizosphere soils. Soil TC, TN, and AK were significantly higher in RT and MP than in NT of both bulk and rhizosphere soils. TP and AP were significantly higher in RT than in MP and NT of bulk soils, while they were remarkably higher in MP and RT than in NT of rhizosphere soils. In addition, MP and RT resulted in higher NO₃^–-N contents in bulk soils compared with that in NT, but there was no significant difference between RT and NT in rhizosphere soils.

The soybean yields were 2,406, 2,826, and 2,592 kg ha^–1 in the NT, RT, and MP, respectively (Supplementary Figure S1). Although the average yield of RT was 17.5 and 9.03% higher than those of NT and MP, respectively, the statistical analysis showed no significant difference in soybean yield among the three tillage practices (Supplementary Figure S1).

Across all samples, Proteobacteria was the dominant phylum, with the relative abundances ranging from 64.6 to 91.9%, followed by high proportions of unclassified diazotrophic sequences ranging from 8.03 to 35.0%. The phyla Actinobacteria, Cyanobacteria, Firmicutes, and Verrucomicrobia occasionally occurred at low frequencies (Figure 1A). Compared with NT, MP significantly increased the relative abundance of Alphaproteobacteria but decreased the relative abundance of Betaproteobacteria, Deltaproteobacteria, Actinobacteria, and Cyanobacteria in bulk soil (Figure 1A). However, RT had little influence on these diazotrophic phyla compared to NT. Except for Verrucomicrobia, the relative abundances of these diazotrophic phyla were not significantly different among the three tillage practices in rhizosphere soil (Figure 1A). At the genus level, Bradyrhizobium dominated the diazotrophic community, and the changes of this genus in response to different tillage practices were consistent with its affiliated phylum (Alphaproteobacteria) in bulk and rhizosphere soils (Figure 1B). In contrast, the three tillage practices had little difference on the relative abundance of Skermanella and Azospirillum in bulk soil, which are both affiliated with Alphaproteobacteria (Figure 1B). Additionally, the correlations between soil properties and diazotrophic taxa showed that soil TC, TN, TP, AP, AK, NH₄⁺-N, and NO₃^–-N significantly influenced the relative abundances of main diazotrophic phyla and genera, expect for Betaproteobacteria. Alphaproteobacteria and Bradyrhizobium were positively correlated with these soil properties, while other diazotrophic phyla and genera had negative correlations with these soil properties (Supplementary Figure S2).

Compared with NT, MP significantly decreased the diazotrophic α-diversity index in bulk soil, while RT showed little difference from NT (Figure 2A). In contrast, the three tillage practices showed no significant difference on α-diversity in rhizosphere soil (Figure 2A). Diazotrophic diversity had significantly negative correlations with the contents of soil TC, TN, C/N, TP, AP, AK, NH₄⁺-N, and NO₃^–-N (Supplementary Figure S3). Additionally, the PCoA plot clearly showed that NT and RT formed a separate group from MP in bulk soil, while no obvious difference in diazotrophic community structures was observed among the tillage practices in rhizosphere soil (Figure 2B). The Adonis analysis revealed that larger differences were detected in NT vs MP (F = 6.728, p = 0.001) and RT vs MP (F = 6.643, p = 0.001) than in NT vs RT (F = 3.698, p = 0.001) as to bulk soil. The Mantel test revealed that all the soil properties examined significantly influenced the variations in diazotrophic community structure, except for TK (Supplementary Table S1).

The diazotrophic empirical networks revealed that the properties, in terms of average clustering coefficient, density, and modularity, were significantly different among the three tillage practices (p < 0.001) (Table 2). Specifically, MP network had fewer nodes and edges than those of NT and RT networks in bulk soils, while little difference in network complexity was detected in rhizosphere soils (Figure 3). The ratio of edge and node, average connectivity, clustering coefficient, and density of RT network were considerably higher than those of NT and MP for both bulk and rhizosphere soils (Table 2). The RT network also had more positive correlations among diazotrophic species and higher connectivity in terms of the distribution of degrees when compared with networks of NT and MP in both bulk and rhizosphere soils (Table 2 and Supplementary Figure S4). The node degree of the diazotrophs’ networks did not increase with the increase of the relative abundance of diazotrophic species, and the nodes with high degrees often had a low relative abundance (Supplementary Figure S5). The variations in soil properties induced by the three tillage practices had more influences on subnetwork in rhizosphere soil than on that in bulk soil (Supplementary Figure S6). Soil pH and AP were the most influential factors in the subnetwork of bulk soil, especially highly correlated with unclassified Proteobacteria and unclassified diazotrophs (Supplementary Table S2). In contrast, soil AK significantly affected the subnetwork in rhizosphere soil and had more correlations with species of Bradyrhizobium (Supplementary Figure S6). Furthermore, the nodes with the highest degree (top five) are displayed in Supplementary Table S3. Unlike networks of NT and MP, the nodes with top degree in the RT network were almost all positively correlated with other nodes in both bulk and rhizosphere soils (Supplementary Table S3). Specifically, OTU10838 (closely related to Bradyrhizobium japonicum), with the highest degree, was simultaneously found in the RT network of both bulk and rhizosphere soils and was positively correlated with other nodes.

No tillage has been commonly reported to increase soil nutrient storage with minimum mechanical soil disturbance compared with conventional tillage (Busari et al., 2015; Dang et al., 2015). However, in most studies, the collection of soil samples from near the surface of soils with crop residue coverage for NT or conventional tillage was conducted with no straw supply (Angers et al., 1997; Melero et al., 2011; Sun B. et al., 2016). Thus, the increased soil nutrients in the NT treatment appeared to be the result of the residue return effect rather than the reduced tillage (Mitchell et al., 2015). In this study, the crop residues were equally returned to the soil for all the three tillage practices, and tillage distributed the soil nutrients uniformly throughout the plow layer. The lower fertility at greater depths of NT resulted in lower contents of detected soil nutrients than in RT and MP (Table 1) because the soil samples of this study were collected from 0 to 20-cm depth rather than near the surface (0–5 cm) and were mixed thoroughly. Despite the reduction in tillage, consistent with the results of Mitchell et al. (2015), RT had similar soil nutrient contents to MP in both bulk and rhizosphere soils (Table 1). Aside from soil nutrients, MP resulted in higher soil pH than NT and RT in both bulk and rhizosphere soils in this research (Table 1). No consistent results have been found in previous studies on changes in soil pH in response to tillage practices (Cookson et al., 2008; Rahman et al., 2008), which might be the reason that tillage had no direct effect on soil pH but depended on soil type, climate conditions, and management factors. Additionally, Martínez et al. (2008) reported that conventional tillage promoted the ability of crop roots to acquire nutrients in deeper soil layers, therefore increasing root penetration and generating higher root mass. Although no significant change in soybean yield was observed among the three tillage practices (Supplementary Figure S1), insufficient access to nutrients might result in lower yields in NT than those in RT and MP, and the reduction in yield with NT remains a major concern (Singh et al., 2011; Pittelkow et al., 2015). In this context, RT is encouraged in terms of its benefits for increasing soil nutrients and improving crop yield compared to NT, and the minor soil disturbance in RT reduced the soil structure deterioration relative to MP.

Biological dinitrogen fixation performed by dinitrogen-fixing microorganisms is an initial process of the nitrogen cycle and the second largest nitrogen source after mineral fertilizers, which contributes up to 16% of the total global N supply (Ollivier et al., 2011). These microbes possess nitrogenase to exercise dinitrogen-fixing functions, in which the nifH gene encodes the subunit of nitrogenase and is highly conserved in prokaryotic microbes. Therefore, it becomes the marker gene of choice for diazotrophic communities in different ecosystems (Zehr et al., 2003; Collavino et al., 2014). In this study, Alphaproteobacteria and the genus Bradyrhizobium within this class dominated the diazotrophic communities (Figure 1), which had been reported to be the major diazotrophs in symbiosis with leguminous plants due to their strong persistence across diverse soil conditions (Pereira e Silva et al., 2013). Although higher soil nutrients were detected in both MP and RT than in NT, MP significantly changed the diazotrophic community composition, while RT had relatively similar diazotrophic communities compared to NT (Figures 1, 2). Compared to NT and RT, MP reduced the abundances of free-living diazotrophs, such as Betaproteobacteria and Cyanobacteria, while it promoted the growth of copiotrophic Alphaproteobacteria due to greater access to soil nutrients under MP (Figure 1) (Zuber and Villamil, 2016; Hutt et al., 2017; Che et al., 2018; Phoma et al., 2018). The differences in nutrient requirements of these diazotrophs were also confirmed by their correlations with soil properties (Supplementary Figure S2). Conversely, NT produced distinct soil stratification structures (Dang et al., 2015), which preserved a higher abundance of Actinobacteria (Figure 1A) that prefer a more stratified environment (Aislabie et al., 2006; Mathew et al., 2012).

Conventional tillage, with periodical soil inversion, leads to the homogenization of nutrient availability and also homogenizes the available microhabitats that have similar and specific microbes across the tillage fraction (Schmidt et al., 2019). Conversely, conservational tillage with no or little soil disturbance increases soil compartmentalization and generates a combination of surface enrichment (due to crop residue coverage) and subsoil nutrient depletion (Dang et al., 2015). These differences in soil stratification resulted in a lower diazotrophic diversity observed in MP than in NT and RT (Figure 2A). The negative influences of conventional tillage on microbial diversity had been reported in previous studies (Diosma et al., 2006; Wortman et al., 2013). The reduction of diazotrophic diversity in MP might be less able to cope with environmental pressures (Wang Y. et al., 2017), which was reflected in negative correlations between diazotrophic diversity and soil properties (Supplementary Figure S3). The changes in microbial community structure had been found to be determined by variations of soil properties in response to tillage practices (Lienhard et al., 2013), which was also detected in the results of the Mantel test of this research (Supplementary Table S1). In contrast to MP, RT did not significantly change the diazotrophic diversity and the community structures relative to NT despite the significant differences in soil properties between NT and RT in this study (Figure 2 and Table 1) or the microbial changes might recover to the status of NT after little disturbance of soils under RT (Wortmann et al., 2008). These results suggested that RT could preserve diazotrophic diversity as high as that of NT and might enhance the potential for microbial nutrient availability in the deeper subsoils, thereby promoting biological dinitrogen fixation (Dang et al., 2015). Noticeably, the changes in diazotrophic communities in response to tillage were mainly observed in bulk soils of this research, and little difference was detected in rhizosphere soils, which might be the stronger determinant of plant roots on diazotrophic communities by the nutrients, exudates, and mucilage released from roots than tillage practices (Philippot et al., 2013).

Microorganisms that coexist in complex interspecies interactions, either negative (e.g., competition) or positive (e.g., cooperation) (Faust et al., 2012), regulate or reflect the structure and the function of ecosystems to some great extent (Fuhrman, 2009). Our study showed that MP simplified the diazotrophs’ network structure relative to NT and RT (Figure 3). The complexity of networks often indicates a coordinated variability of microbial abundance which covary in response to interactions among species or environmental factors (Shi et al., 2016). Although the variations in soil properties showed similar trends in RT and MP, which were different from NT (Table 1), RT resulted in a more complex diazotrophs’ network than MP in bulk soil. This was likely attributed to the diazotrophic species interactions rather than the soil properties controlling the variation of diazotrophs’ networks. Nevertheless, we found that the variations in soil properties induced by the three tillage practices had more effects on diazotrophic interactions in rhizosphere soil than in bulk soil (Supplementary Figure S6). This phenomenon might be expected to spread disturbances quickly through the rhizosphere network once the environment is disrupted and has difficulty in recovering (Fan et al., 2018). In addition, the higher connectivity, clustering coefficient, and density in RT network, than the networks of MP and NT, indicated the much tighter interactions, such that diazotrophic species intensely and efficiently affected each other under RT (He et al., 2017). These network properties could reflect extensively cooperative relationships among members (Banerjee et al., 2018), which was also observed in RT network with more positive correlations than the networks of NT and MP (Table 2). The cooperative characteristic of the diazotrophs’ network in RT improved the fluxes of energy, material, and information with cooperative interactions among diazotrophic assemblages (Rasmann and Turlings, 2016; Zheng et al., 2018). In addition, the extensive cooperation in RT was also reflected in individual nodes, especially in the node with the highest degree (Supplementary Table S3). For example, B. japonicum had more interactions detected in both bulk and rhizosphere soils than any other diazotrophic species (Supplementary Table S3). This diazotrophic species has the ability to adapt to various environments with multiple survival strategies and can improve soybean nodulation and enhance soybean plant growth by co-inoculation of Bacillus strains (Bai et al., 2003). Nevertheless, due to the limited database, the unclassified diazotrophic species should be further examined to better understand the influences of different tillage practices on diazotrophic communities and networks.