Soybean Root Nodule and Rhizosphere Microbiome: Distribution of Rhizobial and Nonrhizobial Endophytes

As anticipated, the most dominant sequences (93% to 99%) retrieved from each of the 193 root nodules were related to genus Bradyrhizobium (Fig. 1 and 2). Regardless of the location (tap or lateral roots) or size of the root nodules (Fig. S4A and B), the relative abundance of Bradyrhizobium sequences was consistently high (>90%) in all nodules (Fig. 1 and 2; Fig. S2; Table 1). Previous studies (4, 12, 43, 53–56) reported that Bradyrhizobium spp. are major microsymbionts within soybean root nodules. However, members of other rhizobial genera such as Mesorhizobium (6), Rhizobium (7, 8), and Ensifer (9–11) were also reported to nodulate soybean. Interestingly, we observed the presence of DNA sequences of other soybean-compatible rhizobia both within root nodules and in the rhizosphere. However, the relative abundance of other rhizobial endophytes compared to Bradyrhizobium was very low (Fig. 3). This observation suggests a preferred selection of Bradyrhizobium strains by the host plant, which was not influenced by location of the root nodule on the host. Likely, soil surrounding each of 193 root nodules would not be identical in terms of factors such as temperature, O₂ concentration, pH, macro- and micronutrients, etc. It is important to mention that we could not measure the geochemical characteristics of soil surrounding each root nodule due to the insufficient amount of soil associated with each root nodule (approximately <0.5 g). For the geochemical analysis, we would need approximately 20 to 30 g of soil. In addition, the collection of soil directly associated with each root nodule was technically not feasible.

Previous studies (10, 53, 57) have reported a high dominance of Ensifer spp. in soybean root nodules in alkaline soils, whereas Bradyrhizobium spp. were detected as major rhizobial endophytes under acidic soil conditions (5). In general, the soils at the study site were highly acidic (pH 5.0 to 5.5), and the pH at this field location was improved through agricultural lime (calcium carbonate) amendment before the sowing of soybean seeds. The final pH of the rhizosphere at the time of plant harvesting (R6 plant growth stage) was slightly acidic (pH 6.5 to 6.8) (Table S4). In the present study, we did not detect any Ensifer-related sequence in 193 nodules, and only 62 Ensifer-related sequences were retrieved from the rhizosphere (Fig. 3).

In the rhizosphere, we observed a relatively high abundance of Bradyrhizobium-related sequences (21,746; 85% of total soil rhizobial sequences) followed by Rhizobium-related sequences (2,130; 8%). Sequences of Mesorhizobium and Neorhizobium were detected in very low abundance (<5%) in the rhizosphere (Fig. 3). Our results suggest that Bradyrhizobium japonicum not only was preferred as the endosymbiont in root nodules but also may have been selectively favored in the rhizosphere. This observation was consistent with the work of Minami et al. (58) and Moawad et al. (59), who suggested that soybean root exudates selectively favored the growth of Bradyrhizobium spp. within the rhizosphere, compared to other rhizobial genera. However, it is important to note that we did not collect bulk soil samples separate from the rhizosphere; therefore, it was not possible to determine whether Bradyrhizobium spp. were selectively enriched in the rhizosphere or were naturally abundant in our soils compared to other rhizobial genera.

We also observed variations in the geochemical characteristics of rhizosphere soils at the three sampling locations. The soil at location 3 had significantly higher concentrations of NO₃-N (26.7 ± 1.8 ppm) and K (134 ± 7 mg/liter) than soil the other two locations (Table S4). Location 1 had a higher concentration of NH₄-N compared to the other two locations. Location 2 had significantly lower concentrations of phosphorus (P) than the other two locations. Despite the significant variations in soil characteristics (Table S4), we did not observe any specific clustering of endophyte community structure based on the sequences originating from the same plant or same sampling location (Fig. 4; Fig. S4). We observed a selection of B. japonicum strains in root nodules at all three sampling locations. The consistent high abundance of Bradyrhizobium in all 193 root nodules suggests that the influence of specific microenvironments such as different nodules on a plant or nodules across three different locations did not affect the selection of rhizobia. Despite differences in the geochemical characteristics, all three locations were in the same field, and therefore, the climatic conditions were the same. In contrast, previous studies suggested a very strong role of soil characteristics in the selection of rhizobial endophytes within soybean root nodules (4, 6, 9, 21, 22). These previous studies were conducted across diverse geographic regions marked by significant differences in soil pH, moisture content, temperature, precipitation, etc. Our finding on the lack of influence of geochemical characteristics on the selection of bacterial endophytes is based on a limited number of samples that were collected from a single field. Further studies that include amending the soil with various concentrations of macronutrients and transplanting different soybean genotypes could improve our understanding of the potential role of variable nutrient concentrations on the selection of bacterial endophytes.

Among the NREs, Nitrobacter and Tardiphaga were detected in almost all individual root nodules (Fig. 1 and 2). Overall, we retrieved 103,707 sequences related to Nitrobacter, and their relative abundance ranged from 1% to 5% within each individual root nodule (Fig. 1). The average abundance of Nitrobacter in the root nodules was <1% of total bacterial sequence; however, a subset of root nodules contained slightly greater levels of Nitrobacter sequences. For example, nine nodules of plant 1 had 3 to 4% Nitrobacter sequences (Fig. 1). This variation did not appear to be related to nodule size, position on root system (tap versus lateral roots), or distance from the soil surface (Fig. 1; Fig. S1 and S4). Tardiphaga-related sequences were detected in 97% of total root nodules with relative densities of 0 to 1%. Both Nitrobacter and Tardiphaga belong to the family Bradyrhizobiaceae and have shown close genetic similarity with the genus Bradyrhizobium based on 16S rRNA (60, 61) as well as five other housekeeping genes, including atpD, dnaK, gyrB, recA, and ropB (61).

We postulate that the consistent presence of Nitrobacter within soybean root nodules may be due to the following reasons. First, Nitrobacter strains are selected within root nodules due to their beneficial role in plant growth. Previously, Ibiene et al. (62) reported that Nitrobacter can solubilize phosphate and have a positive effect on tomato plant growth through inoculation-based studies. Second, the presence of Nitrobacter within soybean root nodules is a parasitic interaction rather than mutualistic. Nitrobacter strains may overcome the host plant defense because of their high genetic similarity to Bradyrhizobium and potentially survive within the root nodules as chemoorganotrophs. Genome sequencing by Starkenburg (60) reported that 41% of genes (1,300 genes) of Nitrobacter strains are identical to those in Bradyrhizobium.

To assess the variation within the Nitrobacter sequences, a maximum-likelihood phylogenetic analysis was conducted. The phylogenetic analysis of all Nitrobacter-related sequences revealed the presence of six major phylogenetic clusters (Fig. S5). Most of the sequences from nodules and rhizosphere samples were grouped into cluster 6 along with other Nitrobacter species sequences from GenBank (Fig. S5).

As with Nitrobacter, we also observed a consistent presence of Tardiphaga-related sequences in almost all root nodules (190 of 193 nodules) (Fig. 2; Table 1). Overall, Tardiphaga-related sequences (20,949) were grouped into three phylogenetic clusters (Fig. S6). Most of the nodule and rhizosphere sequences were related to three reference sequences of Tardiphaga that were isolated by De Meyer et al. (61) from the root nodules of Robinia pseudoacacia (Fig. S6). Previously, Tardiphaga-related strains were isolated from the root nodules of Glycine max (34), R. pseudoacacia (61), and Vavilovia formosa (63) plants. Tardiphaga spp. are genetically similar to the Bradyrhizobium spp. (61) and reported to have nodM and nodT genes (63), which are involved in signaling to the host plant and root nodule formation (64). Moreover, Tardiphaga isolates have been reported to have a role in the N cycle by carrying out dissimilatory nitrate reduction under anaerobic conditions (65).

In the present study, Bradyrhizobium, Nitrobacter, and Tardiphaga were preferred within almost all root nodules, suggesting selection by the host plant or their enhanced ability to penetrate the root nodules. How Nitrobacter and Tardiphaga penetrate and survive within soybean root nodules is an interesting question. We speculate that these bacteria positively affect plant growth, though evidence for this remains to be provided. In all 193 root nodules, Bradyrhizobium-related sequences were the most abundant. Furthermore, none of the 193 root nodules showed high abundance or presence of Nitrobacter and Tardiphaga over Bradyrhizobium. This suggests that Nitrobacter and Tardiphaga enter the root nodules along with Bradyrhizobium and/or that they are not able to overcome the host plant defense in the absence of Bradyrhizobium. The selection, the penetration into host cells or presence within intracellular spaces, and the potential role of these two genera will be the focus of future studies using fluorescence in situ hybridization microscopy and coinoculation-based approaches.

In general, both Nitrobacter- and Tardiphaga-related strains are rarely isolated as plant endophytes in different crops (34, 61–63) (Table S1). In previous culture-based studies on soybean, it was shown that the lack of Nitrobacter may be due to either their slow growth as chemoorganotrophs or lack of NO₂ in the growth medium required for growth as chemolithotrophs. Similarly, Tardiphaga species have a growth rate of up to a 10-day doubling time in culture (63). Our study, along with other culture-independent studies, reports the presence of Nitrobacter and Tardiphaga within root nodules of soybean (34) and alfalfa plants (14, 66).

One of the major goals of this research was to assess the influence of size and location of root nodules on the preferential selection, diversity, and relative abundance of frequently reported NREs (Table S1) within different root nodules of a single plant. In the last few years, several studies have identified many NREs, such as Pseudomonas (14, 16, 31, 32, 39, 40, 51, 67), Variovorax (33, 40), Novosphingobium (40), Flavobacterium (5, 16), Achromobacter (36, 38, 40), Bacillus (5, 16, 27, 32, 35, 37, 46, 67), etc., in soybean root nodules. These studies focused on the isolation and characterization of these NREs for several plant growth-promoting attributes, such as N fixation, phytohormone production, or potential as biocontrol agents (29–32, 40–43, 45, 46, 68). Most of these studies also demonstrated the beneficial role of the isolated NREs in plant growth through subsequent coinoculation studies (27, 30, 44, 67). In the present study, we detected a small number of the sequences (1,635) related to these frequently reported NREs in soybean root nodules (<1% of total sequences). Their presence was inconsistent among different root nodules of a single plant (Table S2). For example, Novosphingobium-related sequences were retrieved from approximately 25% of the 193 nodule samples, and Pseudomonas-related sequences were found in only 16% of root nodules. This was true for almost all other frequently reported (Table S1) NRE genera (Table 1). The overall low abundance and inconsistent presence of these genera among different root nodules of a single plant suggest that these frequently reported NREs (Table S1) are not preferentially selected as endophytes by host plants and most likely have a very limited role in plant growth as endophytes.

The variation in the NRE sequences detected across different root nodules seems to be the result of a random process rather than due to the influence of the location of the root nodule, as root nodules in close proximity showed the presence of very different NREs and vice versa. For example, Pseudomonas-related sequences were detected within three root nodules of plant 1 (root nodule 2 [N2], N6, and N16), and these nodules were apart from each other (Fig. S1; Table S2). N16 was on a lateral root, and the other two nodules (N2 and N6) were on taproot systems. N3, N4, N5, and N7 were in close proximity to N6, but none of them showed the presence of Pseudomonas-related sequences. Overall, this was true for almost all other NRE-related sequences, implying that their presence was not influenced by the size or location of the root nodules (Fig. S1 and S4A and B; Table S2).

It is important to mention that detection of rare NREs present in different root nodules most likely reflects true sequences rather than potential contamination from DNA extraction or PCR reagents. The following three lines of evidence suggest that the observed NRE sequences are true sequences. (i) None of our negative-control samples showed any amplification for first or second PCRs. (ii) These rare NREs were detected in only a few nodule samples (Table 1); if these sequences had originated from contaminating reagents, they would be consistently present across all root nodule samples. (iii) We did detect a high abundance of the identical sequences from rhizosphere soil samples as well. Specifically, more than 25,000 Pseudomonas sequences retrieved from rhizosphere soil were identical to 40 Pseudomonas sequences detected in various root nodules. This was also true for nearly all other rare NREs. Furthermore, we only used the DNA sequences that can be classified at a >80% confidence cutoff level, indicating that they are most likely true sequences rather than sequencing artifacts.

Almost all of the prokaryotic sequences that we detected in root nodules were also found in the 276,338 sequences from nine rhizosphere samples. In contrast to root nodules, there was a higher abundance (2% to 25% of total rhizosphere sequences) of the commonly reported NREs, such as Pseudomonas (11.8% of total rhizosphere sequences), Achromobacter (5.6%), Arthrobacter (3.7%), Planococcaceae (3.5%), Flavobacterium (3.5%), Sphingobacterium (2.9%), Paenibacillus (2.4%), and Stenotrophomonas (2%), in the rhizosphere (Table 1). Based on reports of the positive influence of these NREs on the growth of the host plant, these NREs likely have a significant role in plant growth as free-living plant growth-promoting rhizobacteria (27, 28, 30–32, 43).

In the present study, we have carefully surface-sterilized each root nodule and removed the epidermis. The rare NRE-related sequences (excluding Nitrobacter and Tardiphaga sequences) could originate from surface contamination of root nodules. Even though bacterial cells on the root nodule surface may not be culturable on nutrient agar, their DNA can be amplified through PCR. The surface contamination effect could be more pronounced when multiple root nodules from different plants are pooled for the DNA extraction (14, 16). Likewise, frequent isolation of NREs by culture-based studies could be due to their copious growth on a nutrient-rich growth medium, even if they are naturally present in low abundance in root nodules. Therefore, detecting these commonly reported NREs does not necessarily indicate that they play a significant role in plant growth as endophytes.

In summary, we observed a high abundance of Bradyrhizobium species within 193 root nodules of nine plants (>90% sequences within each root nodule). The two nonrhizobial endophytes (Nitrobacter and Tardiphaga) were also consistently detected within almost all root nodules, though they were in low abundance (1 to 3%) compared to Bradyrhizobium. The consistent presence of these three bacterial genera (one rhizobial genus and two NREs) within almost all root nodules suggests their preferential selection by the host plant that was not influenced by the size, age, or location of the root nodules on a root system. In contrast, we detected low frequency and inconsistent presence of other frequently reported soybean NREs (Table S1) such as Bacillus, Pseudomonas, Flavobacterium, and Variovorax species within very few nodules. In other words, no single NRE that was frequently reported as a soybean endophyte in previous studies (Table S1) was detected consistently or in high abundance in most of the root nodules The overall low abundance and inconsistent presence of frequently reported soybean NREs (Table S1) among different root nodules of a single plant suggest that these frequently reported NREs are not preferentially selected as endophytes by host plants and most likely have a very limited role in plant growth as endophytes. Simple isolation and characterization of frequently reported NREs from root nodules alone may not be sufficient to establish their roles as beneficial endophytes. Their relative abundance within root nodules and preferential selection by the host plant are equally important.

Soybean seeds were obtained from Hancock Seed Co., Dade, FL, USA. The seeds were grown at the Kindrick Farm, Green County, Missouri (latitude, 37°21′44.28″N; longitude, −93°40′91.96″W). Seeds were planted in May 2018, and plants were harvested in September 2018 for the microbiome analysis. We harvested plants at the R6 plant growth stage. Soybean seeds were planted with no tillage at row widths of 15 in. Seeds were planted with planter’s disks, which made small slits for seed placement. The no-tillage approach was adopted to reduce disturbance and soil erosion and to increase available moisture compared to conventionally tilled systems. Seeds were not inoculated with any bacterial culture. Seeds were treated with Acceleron standard seed applied solutions (69). No irrigation was applied, and rainwater was the sole source of water for the crop. The total rainfall from May to September 2018 was 19 in. Rainfall averaged 3.8 in. per month and was lowest during July (1.04 in.) and highest during August (6.71 in.). The average monthly temperature for the 5 months was 25°C ± 2°C. Kindrick Farm has a history of corn cropping for the past 7 years. Three locations within the farm were sampled: location 1 (toe slope), location 2 (sinkhole 1), and location 3 (sinkhole 2). The three soil locations and their geochemical characteristics are described in detail in the supplemental material (Fig. S7; Table S4).

Three replicate plants from three separate rows at each location were excavated, and intact root systems containing root nodules and rhizosphere were brought to the laboratory. We selected all plants of the approximately same height and the same age. Rhizosphere soil was separated from the root systems by gently shaking the plant roots on a sterilized bag. Soil detached in this way was used for geochemical analysis, and soil that remained attached to roots after gentle shaking was used for DNA extraction. For the microbial analysis, root-adhering soil of each plant was washed into separate beakers containing water. The resulting soil suspension was centrifuged at 14,000 rpm for 5 min, and 0.5 g of soil pellet was used for DNA extraction. Roots were then thoroughly washed with distilled water and photographed. Approximately 20 to 24 individual root nodules per plant were selected across the entire root system (Fig. S1). The samples were labeled L1 to L3 for locations, P1 to P9 for plants, and N1 to N24 for nodules per plant (Table S3). We collected the root nodules from both taproots and lateral roots. A taproot, or primary root, is a thick large central root from which lateral roots sprout. The root nodules on the taproot could easily be distinguished from the ones on lateral roots (Fig. S1). A total of 193 individual root nodules from these nine plants were collected for microbiome analysis. Each individual nodule was thoroughly washed in distilled water and surface sterilized in 75% ethanol. The root nodule epidermis was aseptically removed using sterilized forceps while keeping roots submerged in 75% ethanol (Fig. S8). Surface-cleaned root nodules were crushed individually in 1 ml of sterile water using a sterilized mortar and pestle. Crushed root nodule material was used for DNA extraction. Complete surface sterilization of the nodule was confirmed by rolling a few randomly selected nodules on tryptic soy agar plates.

DNA was extracted from each root nodule as well as from the rhizosphere samples using the DNeasy PowerLyzer PowerSoil DNA extraction kit (Qiagen, MoBio, USA). Briefly, the crushed root nodule material was centrifuged at 10,000 rpm for 10 min, and the supernatant was discarded. The pellet containing bacterial cells and plant material was suspended in the Power bead solution (provided with the DNA extraction kit), and DNA was extracted according to the manufacturer’s instructions. For soil samples, 0.4 g of rhizosphere soil was used for DNA extraction. Extracted DNA was eluted in 25 μl of sterilized water and stored at −20°C until used in PCR.

To assess the bacterial communities in the root nodules and rhizosphere, the Illumina MiSeq paired-end DNA sequencing approach was used. Briefly, we used a two-PCR approach. During the first PCR, the partial bacterial 16S rRNA gene fragment of the V3–V5 region was amplified using primers F515 (5′-GTGCCAGCMGCCGCGG-3′) and R907 (5′-CCGTCAATTCMTTTRAGTTT-3′). These universal bacterial primers also contained the sequences for Illumina DNA sequencing primers (Fig. S9). PCR was performed in a 25-μl PCR mixture containing 1× buffer, a 0.2 μM concentration of each primer, 2.0 mM MgSO₄, a 0.2 μM concentration of each deoxynucleoside triphosphate (dNTP), 1.0 μl of template, and 1.0 unit of high-fidelity Platinum Taq polymerase (Invitrogen, USA). The PCR conditions were as follows: 3 min at 96°C, followed by 30 cycles of denaturation at 94°C for 45 s, primer annealing at 56°C for 45 s, extension at 72°C for 45 s, and a final extension at 72°C for 7 min. As a standard PCR procedure, we ran positive and negative controls along with each set of PCRs. Amplified PCR products along with the negative controls were cleaned using the ExoSAP-IT PCR cleanup system (Invitrogen, USA), according to the manufacturer’s instructions.

Cleaned PCR products from the first PCR along with the negative controls were diluted (10 times) in PCR-grade water and used as the template in the second PCR. All the reagents used in the second PCR were as described above with the exception of PCR primers. The primers used in the second PCR contained the Illumina sequencing adapters A and B (A, 5′-CAAGCAGAAGACGGCATACGAGAT-3′; B, 5′-AATGATACGGCGACCACCGAGATCTACA-3′) and standard unique multiplex identifier sequences (listed in Table S5) for each sample (Fig. S9; Table S5). The PCR conditions for this second PCR were as follows: initial denaturation for 3 min at 90°C, followed by 10 cycles of denaturation at 94°C for 30 s, annealing at 60°C for 30 s, and extension at 72°C for 30 s, with a final extension step at 72°C for 7 min. PCR-amplified uniquely indexed PCR products were pooled in equimolar concentrations. Pooled, amplified DNA products were purified with an Agencourt AMPure magnetic bead system (Beckman Coulter, Brea, CA). Purified PCR products were sequenced using the Illumina MiSeq paired-end DNA sequencing platform.

Paired-end Illumina reads were assembled to generate consensus sequences, and initial quality filtration was done using mothur software (70). Initial quality filtrations such as min-length = 370, max-ambig = 3, max-homp = 8, etc., were performed as described previously (71). Chimeric DNA sequences were identified and removed using the Ribosomal Database Project (RDP) platform. Plant-associated cyanobacterium/chloroplast sequences were also identified and removed. High-quality DNA sequences (3,657,224) were classified using the RDP naive Bayesian classifier 2.5 (72). These classified DNA sequences were also aligned and clustered at 97% DNA similarity for multivariate analysis.

DNA sequences related to rhizobial genera were retrieved from root nodules (3,530,886 sequences) and rhizosphere (276,338), and were extracted using mothur. The rhizobial sequences along with the reference sequences from GenBank were aligned and clustered into different groups at 97% DNA similarity. The representative sequences from each cluster were selected for the maximum-likelihood (ML) phylogenetic analysis. The phylogenetic analysis was performed to assess the distribution of rhizobial genera within soybean root nodules and rhizosphere. MEGA software version 5.2 (73) was used to construct an ML phylogenetic tree, and branches representing the same genera were collapsed into subgroups. For the phylogenetic tree, the Tamura-Nei model of uniform nucleotide substitution rate was selected, and confidence in tree topology was assessed by 10,000 bootstrap iterations. Bootstrap values greater than 50% were reported in the phylogenetic tree. Similarly, phylogenetic trees for Nitrobacter and Tardiphaga were also constructed.

The variations in the endophytic bacterial community in soybean root nodules of different plants were accessed by nonmetric multidimensional scaling (NMDS) as described earlier (74). Briefly, we used two different data sets for NMDS analysis. The first data set contained the distribution of all classified bacterial genera in root nodules (containing 3,657,224 sequences), whereas the second data set contained only the distribution of 1,635 sequences of NREs, excluding the sequences related to Bradyrhizobium, Nitrobacter, and Tardiphaga. Within almost all root nodules, these three bacterial genera were consistently detected at a relatively high abundance compared to other commonly reported NREs. For both data sets, we have used the abundance-based Bray-Curtis similarity metric, and data were square root transformed prior to calculating the resemblance metric. The Bray-Curtis similarity metric was calculated using the function vegdist in the R package vegan (75), and the significance of differences across different plants as well as at different locations was tested by PERMANOVA (76) using the function adonis in R (75).

Briefly, the rhizosphere soils of each plant (P1 to P9) were collected separately by gently shaking the plant roots on a sterilized bag. Approximately 20 to 30 g of soil per plant was collected and sent to the University of Missouri Soil and Plant Testing Laboratory for geochemical analysis. The soil samples were tested for pH, organic matter, nitrate concentration, ammonium concentration, phosphorus, potassium, calcium, magnesium, and potential cation exchange capacity using the procedures described in reference 77.

The 16S rRNA gene sequences were deposited in the NCBI Sequence Read Archive (SRA) under BioProject ID PRJNA625756.