A long-term field experiment was established in 2010 with a 4 years crop rotation design located at the Swift Current Research and Development Centre was utilized for this study. The experiment included low (one-pulse phase), medium (two-pulse phases) and high (three-pulse phases) pulse frequency rotations in the long-term rotation design. For this study we selected the following treatments: (i) low pulse frequency rotations, Pea-Wheat-Wheat-Wheat (PWWW) and Chickpea-Wheat-Wheat-Wheat (CWWW); (ii) medium pulse frequency rotations, Pea-Wheat-Pea-Wheat (PWPW) and Chickpea-Wheat-Chickpea-Wheat (CWCW); and (iii) high pulse frequency rotations, Pea-Pea-Pea-Wheat (PPPW), Chickpea-Chickpea-Chickpea-Wheat (CCCW) and Pea-Lentil-Chickpea-Wheat (PLCW) rotations. The first cycle of the rotations was completed in 2013 and repeated again starting in 2014. Due to the development of serious root diseases issues associated with growing chickpea in the same plots for three continuous years during the first cycle, mustard was grown during the second year of the CCCW rotation during the second cycle (i.e., Chickpea-Mustard-Chickpea-Wheat) to limit disease issues. Each treatment was planted into a 4 × 12 m plot and replicated four times in a randomized complete block design. This study was conducted during the 2017 growing season, which was the fourth year of the rotations during the second cycle and all rotations were seeded to wheat (Triticum aestivum L.).

All plots received glyphosate spray at 0.45 kg/ha as a pre-seeding burn-off treatment and Trifluralin 10G at 8.5 kg/ha prior to seeding (except for wheat plots) to eliminate weed effects in each year. At the beginning of each rotation cycle (2010 and 2014) all wheat plots received a total of 80 Kg/ha actual N using 46-0-0 mineral fertilizer, and pulse crop plots received 22 kg/ha of actual P using 11-51-0 mineral fertilizer. All pulse plots also received inoculation prior to seeding with a granular rhizobium inoculant (TeamTag^®, 4kg ha^-1) annually. All plots were under no-till management practices and were planted with a no-till plot seeder.

Soil sampling was conducted at the soft dough stage to assess the soil microbial community and soil chemical properties. Four soil cores (2.5 cm diameter) were randomly collected in each plot to a depth of 15 cm, pooled together and homogenized in a large Ziploc bag to form one composite soil sample per plot. At the same time, four wheat plants were carefully excavated from the same sampling locations as the bulk soil cores in each plot. The loose soil was shaken from the roots and the remaining soil adhering to the roots was carefully brushed off, pooled together and homogenized as one rhizosphere soil sample per plot. There were a total of 28 bulk soil and 28 rhizosphere soil samples collected for this study. Sub-samples were immediately taken from the composite bulk and rhizosphere soil samples and flash frozen in a liquid nitrogen cryo-shipper for molecular analyses. The remaining soil was stored in a cooler with ice for transportation and then stored at 4°C for up to 48 h until further processing. Soil samples were sieved through a 2 mm sieve, and a subsample was used to determine soil moisture content (gravimetric). The remaining soil was air-dried and ground for chemical analyses.

For all soil samples, soil organic carbon and total carbon were determined using the dry combustion method (after acidification with HCl) using an Elementar vario MICRO cube elemental analyzer (Schumacher, 2002). Soil nitrate nitrogen (N), phosphate phosphorus (P) and potassium (K) were determined using sodium bicarbonate extractions followed by colorimetric analysis using a Technicon Autoanalyzer (Harm et al., 1973; Gentry and Willis, 1988). Soil sulfate sulfur (S) was determined using calcium chloride extractions followed by colorimetric analysis using a Technicon Autoanalyzer (Harm et al., 1973). Soil pH and EC were measured in water saturation paste (Habekost et al., 2008) and paste extracts (Miller and Curtin, 2006). Soil moisture was measured using the gravimetric method.

Genomic DNA was extracted from 0.25 g of soil in duplicate for each bulk and rhizosphere soil sample using the DNeasy PowerSoil kit (Qiagen) using the QIAcube automated nucleic acid extraction system (Qiagen). DNA was quantified using a Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific, Waltham, MA) and shipped on dry ice to the Genome Quebec Innovation Center (Montreal, Canada) for amplicon library preparation and Illumina MiSeq sequencing. The bacterial 16S rRNA genes were sequenced using primers 515-F and 806-R (Caporaso et al., 2012) and fungal ITS1 region was sequenced using primers ITS1F and 58A2R (Martin and Rygiewicz, 2005). For a full description of the amplicon library preparation and Illumina Miseq sequencing see Delavaux et al. (2020). The raw amplicon sequencing datasets are available in the NCBI Sequence Read Archive under BioProject ID PRJNA694704, PRJNA694839, PRJNA694691, PRJNA694475, PRJNA717779, and PRJNA694448. Raw reads were processed (assembled, trimmed, quality filtered, dereplicated, and clustered) using the UPARSE pipeline of USEARCH v.9 (Edgar, 2013). Sequences were clustered into operational taxonomic units (OTU) based on 97% similarity. Taxonomic identity was assigned with a confidence threshold of 80% using the RDP classifier (Wang et al., 2007) and 16S rRNA training set (version 16) for bacteria and ITS UNITE database for fungi (Kõljalg et al., 2013). FUNGuild was used to parse the fungal community datasets by functional ecological guild (Nguyen et al., 2016). See Bainard et al. (2020) for a full description of the bioinformatics.

All statistical analyses were performed using R packages (version 3.5.1) vegan, breakaway, ecodist and ggplot2 (R Foundation for Statistical Computing, Vienna, Austria). To determine the effect of pulse frequency on α-diversity we estimated the total OTU richness using the package breakaway (Willis and Bunge, 2015). One-way analysis of variance (ANOVA) was used to test the effect of pulse frequency on the soil properties and OTU richness. Post hoc comparisons of means were completed using Tukey’s HSD procedure at a 5% level of probability. Permutational multivariate analysis of variance (permanova) was used to test the effect of pulse frequency on the wheat rhizosphere and bulk soil microbial communities. Sequence read variability was normalized among samples by randomly resampling to the lowest number of reads for each dataset (i.e., bacteria = 19,682 reads/sample; fungi = 21,468/sample). To eliminate influences of low abundance OTUs which are in many cases artifacts generated through the sequencing process and better compare treatment effects, unclassified OTUs and rare OTUs (OTUs shown in less than 3 soil samples and OTUs with less than 10 reads in total) were removed from downstream analysis. Network analysis was constructed for soil microbial communities related to different pulse frequencies based on their OTU relative abundances, and was performed using the Molecular Ecological Network Analyses (MENA) pipeline: http://ieg4.rccc.ou.edu/mena (Deng et al., 2012) and visualized by Cytoscape 3.6.1. Module-EigenGene network analysis was applied to examine interactions of the soil microbial community with environmental factors. The connectivity of each node was determined based on its within-module connectivity (Zi) and among-module connectivity (Pi), and the node topologies were separated by Zi and Pi with threshold of 2.5 and 0.62 (Guimerà and Nunes Amaral, 2005; Olesen et al., 2007). OTUs with low Zi and low Pi values were considered as specialists (peripherals) which typically only had links within their own modules. OTUs with low Zi and high Pi values (connectors), or high Zi and low Pi values (module hubs) were considered as generalists. This included either module hubs (actively linked with other OTUs within their own module) or connectors (linking with OTUs from other modules). OTUs with high Zi and Pi values were considered as super generalists (network hubs) which act as both connectors among several modules and module hubs within their own modules.

Pulse frequency in the rotations had a minimal impact on the bulk and rhizosphere soil properties. Total C and organic C in rhizosphere soil was significantly higher in rotations with medium pulse frequency compared to rotations with low and high pulse frequency, and soil moisture was highest in rotations with low pulse frequency (Table 1). In the bulk soil, rotations with a medium pulse frequency had a significantly higher soil S and P content compared to the high pulse frequency rotations. All other soil properties in the bulk and rhizosphere soil were unaffected by pulse frequency. Pulse frequency also impacted wheat growth. Tiller number and dry biomass of wheat were significantly affected by pulse frequency, with rotations with medium pulse frequency exhibiting the highest biomass production. However, wheat yield, height or seeds per plant did not significantly differ among the rotations with different pulse frequencies (Supplementary Table 1).

Pulse frequency did not significantly affect the α–diversity (estimated total richness) of the bacterial community in the rhizosphere or bulk soil (Figures 1A,B). However, pulse frequency did significantly affect the estimated total fungal richness in the rhizosphere (highest under low pulse frequency), but not in the bulk soil (Figures 1C,D). Among the fungal guilds, pathotrophs were the only group that were significantly affected by pulse frequency (Table 2). In the rhizosphere, rotations with medium pulse frequency had the highest proportion of fungal pathotrophs. In the bulk soil, rotations with medium and high frequency had the highest proportion of fungal pathotrophs. Of the fungal pathotrophs putatively assigned as plant pathogens in the bulk soil, the rotations with medium and high pulse frequency had nearly twice the proportion of reads compared to the rotations with low pulse frequency.

The permanova indicated that pulse frequency had a significant effect on the fungal community in both the rhizosphere and bulk soil (Table 3). This appears to be primarily linked to shifts in the taxa belonging to the Ascomycota in bulk soil, and taxa belonging to the Ascomycota and Mortierellomycota in rhizosphere soil. The principle coordinate analyses (PCoA) revealed that the fungal community from rotations with low pulse frequency were distinct and showed some separation compared to those from medium and high pulse frequency in both the rhizosphere and bulk soil (Supplementary Figure 1). This trend is particularly evident for the Ascomycete community in both the rhizosphere and bulk soil. Pulse frequency did not have a significant effect on the bacterial community in either the bulk or rhizosphere soil. However, pulse frequency did significantly influence the rhizosphere Firmicutes community, and bulk soil bacterial communities belonging to the Actinobacteria, Armatimonadetes, Bacteroidetes, candidate division WPS-1, Gemmatimonadetes, and Proteobacteria (Table 3 and Supplementary Figure 2). The strongest effect was observed on the Bacteroidetes community in the bulk soil, with the PCoA showing a distinct separation between samples taken from rotations with low pulse frequency and those from medium and high pulse frequency (Supplementary Figure 2).

The network analysis revealed that increasing the pulse frequency in crop rotations modified the interactions of the soil microbial community (Table 4). The top three dominant bacterial phyla (Acidobacteria, Actinobacteria, and Proteobacteria) influenced the structure and interconnectedness of the soil microbial community. The majority of the nodes with high degree in the network analysis in both rhizosphere and bulk soil were linked to these three dominant phyla (Acidobacteria > Actinobacteria > Proteobacteria), and accounted for approximately 50% of the identified dominant nodes (Figures 2, 3). In both rhizosphere and bulk soil, rotations with high pulse frequency had microbial communities with stronger interactions that included a smaller number of nodes, higher number of edges, shorter average path distance, more connectedness, higher betweenness centrality, less modules and more dominant nodes (Table 4). These results suggest that high levels of pulse frequency in crop rotations can enhance the soil microbial community interactions.

Increasing pulse frequency also increased phylogenetic diversity of dominant nodes in the rhizosphere (Figure 2 and Supplementary Figure 3), especially taxa that belonged to the Basidiomycota, Chloroflexi, Thaumarchaeota and candidate division WPS-1. The number of nodes from Ascomycota, Gemmatimonadetes, and Proteobacteria were increased while the number of nodes from Acidobacteria and Actinobacteria were decreased with increased pulse frequency. We also observed a higher number of negative interactions between nodes under high pulse frequency than under low pulse frequency, which may be linked to higher competition between soil microorganisms for limited nutrients and environmental resources, or accumulation of pulse specific pathogens.

In the bulk soil, high pulse frequency led to an increased total phylogenetic diversity of dominant nodes along with a shift in composition compared with low and medium pulse frequency rotations (Figure 3 and Supplementary Figure 4). For example, dominant nodes from Acidobacteria, Actinobacteria, and Gemmatimonadetes were reduced and dominant nodes from Ascomycota, Basidiomycota, Mortierellomycota, Planctomycetes and candidate division WPS-1 were increased in rotations with high pulse frequency. The latter group accounted for approximately half of the total dominant nodes and the majority of the negative interactions with other nodes in the bulk soil microbial community.

Module-EigenGene network analysis grouped soil microorganisms into different modules based on their interactions with environmental factors (Supplementary Figure 5). In both the rhizosphere and bulk soil, module numbers decreased with an increase in pulse frequency, indicating that rotations with higher frequency of pulses can strength interactions within the microbial community. However, the microbial community in bulk soil and rhizosphere soil reacted to the pulse frequency differently.

In the rhizosphere, S, pH and EC had the strongest influence on the soil microbial community (Supplementary Figures 5A–C and Supplementary Table 1). This was particularly evident for the bacterial community under low to medium pulse frequency rotations. The modules that contained the keystone OTUs (black framed modules, Supplementary Figures 5A–C) were more sensitive to EC and pH under low to medium pulse frequency rotations, but more sensitive to TC and TN under high pulse frequency. This indicated that the influence of pulse crops on microbial keystone taxa in the rhizosphere is co-dependent on the pulse frequency and environmental factors.

In the bulk soil, S, N, and K had the strongest influence on the soil microbial community in low pulse frequency rotations (Supplementary Figure 5D), while N had the strongest influence in medium pulse frequency rotations (Supplementary Figure 5E). Under high pulse frequency rotations, EC had the strongest effect (Supplementary Figure 5F), particularly for the bacterial community (Supplementary Table 1). For the modules that include most of the keystone taxa, they had a negative response to S under low pulse frequency, and strong response to N (either negative or positive, depending on the module) under medium pulse frequency, and none of these environmental factors significantly affected these modules under high pulse frequency. These results suggest that in bulk soil, rotations with a high frequency of pulse crops could develop a relatively stable soil microbial community network that is not as responsive to changing environmental conditions.

To evaluate the topological role of identified nodes in the network, we used within-module connectivity (Zi) and among-module connectivity (Pi) values to separate each network into four ecological modules (Figure 4). From the rhizosphere soil (Figure 4A), we identified 22 module hubs and 4 connectors in low pulse frequency, 20 module hubs and 17 connectors in medium pulse frequency, 25 module hubs and 11 connectors in high pulse frequency. From bulk soil (Figure 4B), we identified 23 module hubs and 2 connectors in low pulse frequency, 18 module hubs and 8 connectors in medium pulse frequency, 14 module hubs, 74 connectors, and 9 network hubs in high pulse frequency. The higher number of connectors under high pulse frequency indicated that soil microbial taxa were more active within their own modules under low pulse frequency, but tended to build up interactions with OTUs located at different modules under high pulse frequency. These results suggest that higher levels of pulse frequency in crop rotations can enhance the interactions of soil microbial taxa in different subgroups which may carryout different biological and ecological functions.