Genetic diversity of rhizobial strains isolated from M. polymorpha and M. lupulina nodules
Totally 91 rhizobial isolates were obtained from 39 sampling sites in 7 locations (Fig. S1), in which 70 were from M. polymorpha and 21 were from M. lupulina (Table 1). All of these strains were divided into five rDNA types by IGS–RFLP analysis and displayed 39 BOX-AIR fingerprinting patterns (Table 1). The IGS type A containing 19 BOX-AIR patterns was mainly composed of 63 strains from M. polymorpha (Table 1) collected in five locations, as well as 4 strains from M. lupulina sampled in Yiliang, Kunming. The IGS type C represented 24 strains from both M. polymorpha and M. lupulina with 17 BOX-AIR patterns. In contrast, the IGS types B, D and E each contained only one or two strains, and each strain has its own specific BOX-AIR pattern.
Phylogenies based on housekeeping and symbiotic genes
Based on the BOX-PCR patterns as well as the IGS-RFLP profiles of all the bacterial isolates, 25 representatives were selected for sequencing the 16S rDNA and five housekeeping genes (recA, glnII, atpD, dnaK and glnA) to determine their species affiliation (Fig. S2-S7).
Based on the 16S rDNA phylogeny, the representatives were divided into two clades, with 10 isolates closely related to E. meliloti LMG6133T (>99% similarity) and the remaining 15 clustered with S. medicae WSM A321T (>99% similarity) (Fig. S2). The phylogenetic tree based on MLSA of the 5 concatenated housekeeping genes (Fig. 1) showed a similar topology to that of the 16S rDNA tree, and similarities greater than 97% were observed among the strains within each of the E. meliloti and S. medicae clades (including the type strain in each). Thus, these two clades were identified as E. meliloti and S. medicae, respectively. E. medicae covered 63 test strains belonging to the IGS types A and B, while E. meliloti clade included 28 strains in IGS types C, D and E. E. medicae was predominated in M. polymorpha isolates (occupied 85.7%, 60/70), while E. meliloti was dominant in M. lupulina (occupied 80.9%, 17/21).
Intriguingly, some E. meliloti strains exhibited traces of genetic incongruence, which were evidenced by comparison of the housekeeping gene phylogenies of MLSA. In the phylogenic tree of aptD (Fig. S3), the 15 E. medicae representative strains together with the type strain clustered in one clade that was consistent to the result of MLSA. Nevertheless, the other 10 isolates formed 3 lineages within the clade intermingling with E. meliloti and E. kummerowiae type strains. Similar relationships were also observed in the recA phylogeny (Fig. S5), and the genetic incongruence between aptD/recA phylogeny and MLSA phylogeny suggested that the representative strains SWF67487 and SWF67523 might acquire genetic materials via horizontal gene transfer from S. kummerowiae strains during the natural evolution process. The gene transfer event of E. meliloti isolates from E. kummerowiae can be further supported by the phylogeny of symbiotic genes, since the 25 representatives were divided into 3 groups in the nodC phylogeny (Fig. 2), which were different from the results of MLSA. Particularly, 8 E. meliloti strains formed a monoclade together with S. kummerowiae CCBAU 71714T at 100% identity. Likewise, these 8 strains were also clustered together with S. kummerowiae CCBAU 71714T in the nifH phylogeny (Fig. S8).
Taken together, the phylogenetic analyses based on multiple concatenated housekeeping gene sequences (MLSA) divided the 25 representatives into two groups, corresponding to E. meliloti and E. medicae. Traces of HGT events in some E. meliloti strains were evidenced by the differences between phylogenies of single housekeeping genes, and between the phylogenies of MSLA and symbiotic genes.
Symbiotic performance of representative strains on Medicago polymorpha
Given the colonial morphology, hosts, and the genomic/phylogenetic analyses, 32 representative strains belonged to E. meliloti and E. medicae from both Medicago species were chosen for nodulation assays (Table 1). All the selected strains formed nodules on roots of M. polymorpha. Based on the patterns of BOX-PCR of nodule crushes, most the recovered strains (31/32) were identical to the corresponding inoculants, despite their host of origin. Thus, they were undoubtedly defined as the microsymbionts of M. polymorpha. Only strain SWF67450 was not recovered from nodules on plants inoculated with this strain, which might be explained that the plants in this test were partially or totally contaminated by other strains. Notably, two E. meliloti strains, SWF65114 and SWF65115, originally isolated from M. lupulina plants were capable of inducing efficient nodules on M. polymorpha as well.
Biogeographical patterns and genetic differentiation of Medicago rhizobia
The aforementioned phylogenic analysis indicated genetic incongruence in E. meliloti strains but not in E. medicae representatives. Then, we assessed the gene diversity of the rhizobial isolates by assessing the BOX AIR pattern and Shannon index. E. meliloti populations showed significantly higher gene diversity index than the E. medicae counterparts. Moreover, a high level of genetic differentiation (Fst=0.92054) was observed between E. meliloti and E. medicae, based on the sequences of 5 housekeeping genes (Table 2). To know whether host specificity would affect genetic diversity of rhizobial isolates, we firstly compared the genetic diversity of symbiotic strains within the Ensifer species. The total E. medicae strains from M. polymorpha in both farmland and nature ecosystem showed higher genetic diversity to those from M. lupulina nodules according to the BOX-PCR pattern and Shannon index. Similarly, E. meliloti isolates from M. lupulina showed significantly higher genetic diversity than those from M. polymorpha nodules (Fig. 3). Thus, these data suggested that rhizobial strains from their native host tend to be more genetically diverse.
Then, we evaluated whether soil type would affect rhizobial diversity by performing the PCoA based on Bray-Curtis distance (Table S2, Fig. 4). The amount of multiple nutrient factors was extensively various in natural ecosystem while it tended to be more similar in the farmland sites, as the observation that the component 1 represented 53.4% of relative eigenvalues, and component 2 represented 20.04% of relative eigenvalues, after Cailliez correction. Furthermore, multivariate analysis of variance (MANOVA) revealed that the E. medicae distribution significantly differed between the two habitats (P=0.047), while no overt difference was observed for E. meliloti species between farmland and nature ecosystem. In addition, permutational multivariate analysis of variance (PERMANOVA) corroborated that soil type (farmland/natural ecosystem) accounted for 1.61% (P = 0.15) of variation in the observed beta-diversity of rhizobia (Bray-Curtis distance metric).
On the one hand, all the Ensifer strains originated from natural ecosystem had a higher level of nucleotide diversity than those from farmland (Table 2), suggesting that the soil type might contribute to rhizobial genetic diversity. The greater variation in the content of several nutrients in natural ecosystem might explain the higher genetic diversity of rhizobia. Consistently, moderate level of genetic differentiation was observed within E. medicae (Fst=0.05251) or E. meliloti (Fst=0.0501) between natural ecosystem and farmland. Furthermore, gene flow was detected for both E. meliloti and E. medicae strains between farmland and natural ecosystem. Altogether, these data suggested that Ensifer strains from natural ecosystem were more genetically diverse compared to those from farmland.
Deterministic factors for diversity of rhizobia from two Medicago plants
To identify which specific soil factor could determine the genetic diversity of rhizobial populations, analysis of variance (ANOVA) based on a set of environmental factors was performed. Compared to natural ecosystem, the alfalfa farmland had a significantly higher content of multiple nutrient factors, including organic matter (OM), sodium (Na+), magnesium (Mg2+) and bicarbonate (HCO3-) (Table S3, S4). As shown in the CAP analysis based on Bray-Curtis distance of PA (presence/absence) transformed species data (Fig. 5A), the soil physiochemical data could explain 53.4% of variation in the diversity of five IGS type species. Soil Ca2+, Na+, HCO3-, Cl-, SO42- and pH contributed to a significant partition in the first component (63.51% of total variance). Meanwhile, nitrogen (N), OM, potassium (K), and phosphorus (P) played a minor role in the second component (14.31 % of total variance), showing negative effects on rhizobial distribution. Variation partition analysis revealed that 11.7%, 5.8%, 2.2% and 0.5% could be significantly explained by 6 ions, NPK, pH and OM variables respectively. These results indicated that Ca2+ and Na+ were major factors in shaping the genetic diversity of five IGS types (Fig. 5B).
Rhizobial species correlated with different sets of explanatory variables were further identified. IGS type A was more likely found in the farmland sites with lower HCO3-. IGS type C was more likely found in the nature ecosystem with high soil Ca2+, Na+, HCO3-, Cl-, lower nitrogen, K and OM. The IGS types B, D and E more likely appeared in the nature ecosystem sites with high P and low HCO3-. In summary, IGS types A and C have similar nutrition utilization, and both preferred soils with low potassium, OM, phosphorate and nitrogen contents, while IGS type C preferred soil with high contents of Ca2+ and Na+ (Fig. 5A).
To evaluate the correlation of rhizobial abundance with environment factors, CAP analysis based on Hellinger transformed species data was conducted (Fig. 5C). IGS type A had a significantly higher abundance in farmlands that were positively associated with pH. While more stringent than type A, type C was also negatively correlated with soil N, k, P and OM. IGS types B, D and E were not considered due to their extremely low abundance (Fig. 5C). Variation partition analysis of Hellinger transformed species data (Fig. 5D) exhibited almost the same pattern as the PA transformed species data (Fig. 5B). The range of soil physiochemical content for two rhizobial species identified in this study was given in Table S4.
Taken together, distribution of E. medicae strains was positively correlated with Na+ in soil, and their abundance was also associated with low level of Nit and OM. E. meliloti isolates appeared to prefer good-quality soil conditions with high Ca2+ and Na+ contents, and negatively correlated with soil N and OM. To conclude, the soil contents of N, OM as well as Ca2+ and Na+ were the major soil factors to shape the distribution of rhizobial strains belonging to the two Ensifer species detected in the present study.