Cultures of R. faveluksii OR191 were routinely sub-cultured on TYC solid medium (Howieson and Dilworth, 2016) incubated at 28°C for 3–4 days. For long-term maintenance, bacterial strains were grown in TYC broth and preserved in 20% glycerol at –80°C.

R. favelukesii OR191 was streaked onto TYC solid medium and grown at 28°C for 3 days to obtain well grown, well separated colonies. A single colony was selected and used to inoculate 5 ml TYC broth medium. The culture was grown for 48 h on a gyratory shaker (200 rpm) at 28°C. Subsequently 1 ml was used to inoculate 60 ml TYC broth medium and grown on a gyratory shaker (200 rpm) at 28°C until an OD_600nm of 0.6 was reached. DNA was isolated from 60 ml of cells using a CTAB bacterial genomic DNA isolation method (Joint Genome Institute, 2022). Final concentration of the DNA was set to 0.5 mg ml^–1.

The draft genome of Rhizobium favelukesii OR191 was generated at the DOE Joint Genome Institute (JGI) using Illumina data (Bennett, 2004). An Illumina standard shotgun library was constructed and sequenced using the Illumina HiSeq 2000 platform which generated 17,712,488 reads totaling 2,657 Mbp. All raw Illumina sequence data were passed through DUK, a filtering program developed at JGI, which removes known Illumina sequencing and library preparation artifacts (Mingkun, L., Copeland, A. and Han, J., unpublished). The following steps were then performed for assembly: (1) filtered Illumina reads were assembled using Velvet (version 1.1.04) (Zerbino, 2010), (2) 1–3 Kbp simulated paired end reads were created from Velvet contigs using wgsim (Github and GitHub, 2011), (3) Illumina reads were assembled with simulated read pairs using Allpaths–LG (version r39750) (Gnerre et al., 2011). Parameters for the assembly steps were 1) Velvet: –v –s 51 –e 71 –i 2 –t 1 –f “-shortPaired -fastq $FASTQ” –o “-ins_length 250 -min_contig_lgth 500” for Velvet and 2) wgsim: -e 0 -1 76 -2 76 -r 0 -R 0 -X 0. The final draft assembly contained 240 contigs in 240 scaffolds. The total size of the assembly is 7.4 Mbp with an average of 350 × coverage of the genome.

Genes were identified using Prodigal (Hyatt et al., 2010), as part of the DOE-JGI genome annotation pipeline (Mavromatis et al., 2009; Chen et al., 2013). The predicted CDSs were translated and used to search the National Center for Biotechnology Information (NCBI) non-redundant database, UniProt, TIGRFam, Pfam, KEGG, COG, and InterPro databases. The tRNAScanSE tool (Lowe and Eddy, 1997) was used to find tRNA genes, whereas ribosomal RNA genes were found by searches against models of the ribosomal RNA genes built from SILVA (Pruesse et al., 2007). Other non–coding RNAs such as the RNA components of the protein secretion complex and the RNase P were identified by searching the genome for the corresponding Rfam profiles using INFERNAL (Nawrocki and Eddy, 2013). Additional gene prediction analysis and manual functional annotation was performed within the IMG-Expert Review system (Markowitz et al., 2009) developed by the Joint Genome Institute, Berkeley, CA, United States.

We assessed the phylogenetic position of R. favelukesii OR191, relative to other Rhizobium type and non-type strains, using a 1,296 bp internal region of the 16S rRNA gene. The Azorhizobium caulinodans ORS 571^T sequence was used as an outgroup. Phylogenetic analyses were performed using MEGA X (Kumar et al., 2018). The tree was built using the maximum likelihood method with the Tamura 3- parameter model (Tamura, 1992). A discrete Gamma distribution was used to model evolutionary rate differences among sites [5 categories (+ G, parameter = 0.3233)]. The rate variation model allowed for some sites to be evolutionarily invariable [( + I), 42.2407% sites]. Bootstrap analysis (Felsenstein, 1985) with 500 replicates was performed to assess the support of the clusters.

For phylogenetic analysis of the NodA proteins, where more than one nodA allele was present in the genome, we used the translated amino acid sequence of the nodA gene that was immediately upstream of nodBC. The evolutionary history was inferred by using the Maximum Likelihood method based on the Tamura 3-parameter model. The tree with the highest log likelihood (4204.9163) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites [5 categories (+ G, parameter: 3.2381)]. The rate variation model, used for nodA, allowed for some sites to be evolutionarily invariable [( + I), 18.6224% sites]. The trees were drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 41 nucleotide sequences for nodA. All positions containing gaps and missing data were eliminated. In the final dataset there were a total of 588. Evolutionary analyses were conducted in MEGA X (113).

Phylogenetic analysis and evolutionary history of the NifH proteins was inferred by using the Maximum Likelihood method and JTT matrix-based model (Jones et al., 1992). The tree with the highest log likelihood (–8402.67) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the JTT model, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites [5 categories (+ G, parameter = 2.4848)]. The rate variation model allowed for some sites to be evolutionarily invariable [( + I), 14.51% sites]. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. This analysis involved 38 amino acid sequences; Rhizobium mesoamericanum STM3655 and Ensifer meliloti AK53 were excluded from the analysis due to a lack of sequencing information. There were a total of 572 positions in the final dataset. Evolutionary analyses were conducted in MEGA X (Kumar et al., 2018).

gANI comparisons were performed using the IMG gANI analysis tool (Chen et al., 2021).

Cells of OR191 and Ensifer medicae WSM419 were grown to mid-exponential phase at pH 7.0 in JMM minimal salts media (Howieson and Dilworth, 2016), washed in JMM at pH 7.0 and then resuspended in JMM at either pH 7 or 5.8 to an OD_{600 nm} of approximately 0.01 for pH 7.0 media (100 ml) and 0.25 for pH 5.8 media (100 ml). Cells were incubated at 28°C with shaking at 200 rpm for 4 days. Cells were then removed by centrifugation (10 min at 7,000 g) and EPS was precipitated from the supernatant by adding 0.1 vol. of a 2% hexadecyltrimethyl-ammonium bromide solution stored at 28°C. EPS was pelleted (10 min at 7,000 g) and redissolved in 10% (w/v) NaCl. The EPS was precipitated by adding 2 volumes of acetone, pelleted (10 min at 7,000 g), redissolved in sterile water and dialyzed in water. EPS was lyophilized, resuspended in D₂O (99.96%.), lyophilized again, and dissolved in D₂O (99.96%.) to a concentration of 15 mg ml^–1. ¹H-NMR spectra were recorded on a Varian 400-MR spectrometer at 80°C. A 30 μl aliquot of a 2% (w/v) solution of the succinoglycan-binding dye Calcofluor (Fluorescent brightener 28, Sigma) was added to 1 ml samples containing 7.5 mg EPS in saline [0.89% (w/v)] to visualize fluorescence at a wavelength of 365 nm.

OR191 is a rod shaped isolate that is fast growing and forms typically mucoid colonies on conventional media used to isolate rhizobia (Howieson and Dilworth, 2016; Supplementary Figure 1). Flagella are not present in the transmission electron micrograph of OR191 (Supplementary Figure 1B), which is consistent with the description of the R. favelukesii type strain LPU83^T as non-motile (Torres Tejerizo et al., 2016). The classification, general features and genome sequencing project information for OR191 are provided in Supplementary Table 1, in accordance with the minimum information about a genome sequence (MIGS) recommendations (Field et al., 2008) published by the Genomic Standards Consortium (Field et al., 2011).

The draft genome of strain OR191 was generated at the DOE Joint Genome Institute (JGI) using Illumina data (see “Materials and Methods” section). A summary of the genome project features for OR191 is shown in Table 1. The genome is 7,368,160 bp with 59.66% GC content and is comprised of 240 scaffolds. From a total of 7,704 genes, 7,617 were protein encoding and 87 RNA-only encoding genes. The majority of genes (72.12%) were assigned a putative function whilst the remaining genes were annotated as hypothetical. The distribution of genes into COGs functional categories is presented in Supplementary Table 2.

Many of the features of OR191 are indistinguishable from those of R. favelukesii LPU83^T and T1155, including host range, nod and nif genes, insertion-sequence hybridization profiles, PCR-based genomic fingerprints, and plasmid profiles (Eardly et al., 1985, 1992; Del Papa et al., 1999; Segundo et al., 1999; Wegener et al., 2001; Hou et al., 2009; Bromfield et al., 2010; Wibberg et al., 2014; Table 2). Consistent with the plasmid profiles, we identified genes encoding homologs of the Rep proteins (required for the replication and stable maintenance of plasmid replicons) in the OR191 genome that shared 100% identity with those present in the R. favelukesii pLPU83a, pLPU83b, pLPU83c, and pLPU83d plasmids, suggesting that both R. favelukesii strains possess the same number of replicons. Furthermore, the OR191 genome contains two separate regions that encode proteins required for DNA transfer replication and mating pair formation, which share very high identity with homologs in the characterized plasmids pLPU83a (which is conjugative) and pLPU83b (which is mobilizable) (Torres Tejerizo et al., 2010, 2014). We also found genes encoding a separate conjugative system in OR191 that shares very high identity with a conjugative system that is present on the pLPU83d chromid. LPU83^T shares its lipopolysaccharide profile with OR191 (Wegener et al., 2001), while T1155 shares its phage-resistance type with OR191 (Bromfield et al., 2010). Furthermore, the DNA-DNA hybridization values between strains LPU83^T and OR191 were greater than 84%, indicating that these strains are members of the same species (Torres Tejerizo et al., 2016). DNA alignments of OR191 essential “housekeeping” gene sequences (including 16S rRNA, 23S rRNA, recA, and atpD genes) to the corresponding genes of LPU83^T, T1155 and the closely related Rhizobium tibeticum CGMCC 1.7071^T gave 100% sequence identity of OR191 genes to those of LPU83^T and from 98.83 to 100% sequence identity to those of T1155 or CGMCC 1.7071^T (Table 2). In the 16S rRNA phylogenetic tree, OR191 is in a clade with R. grahamii CCGE 502^T, R. favelukesii strains and R. tibeticum CGMCC 1.7071^T (Supplementary Figure 2).

gANI pairwise comparisons were also calculated for Rhizobium genomes deposited in IMG that were most closely related to OR191 in the 16S rRNA phylogenetic tree. The genome of R. favelukesii T1155 could not be used in this analysis as it has not yet been sequenced. gANI values for the most closely related strains are shown in Table 3; scores over the defined species affiliation cut-off value of > 96.5 gANI (Varghese et al., 2015) are shown in bold font. The gANI pairwise comparisons, along with other previously presented phenotypic and genotypic data, confirmed that OR191 is conspecific with R. favelukesii LPU83^T. The gANI results also suggest that R. tibeticum CGMCC 1.7071^T belongs to the same species as LPU83^T and OR191. However, the DND-DNA hybridization values obtained by Torres Tejerizo et al. (2016) indicate that R. favelukesii and R. tibeticum strains are separate species; moreover, R. tibeticum strain CGMCC 1.7071^T shared lower sequence identities of the recA, atpD, and rpoB genes than the R. favelukesii strains and differed in its plasmid profile, fatty acid profile and ability to metabolize several sole carbon substrates.

The genes of the sequenced R. favelukesii strains OR191 and LPU83^T were compared, using the default parameters of 30% protein identity and E-value of 1e-5 in the JGI IMG phylogenetic profiler. Of a total of 7,704 protein coding genes in OR191, 7,449 genes were found to be in common with LPU83^T. The remaining 255 genes were unique to OR191, with nearly half annotated as hypothetical (128 genes representing 50.20%). These hypothetical genes included extensive prophage loci clustered on scaffold 24.25. By using the prophage predicting tool PHASTER (Arndt et al., 2016) these loci were shown to encode a complete resident prophage. Whole genome analysis of OR191 found four additional resident prophages present on scaffolds 4.5, 7.8, 84.85, and 91.92 (Supplementary Figure 3), in addition to the unique prophage on scaffold 24.25.

Strain OR191 was isolated from alfalfa growing in moderately acidic soil (Eardly et al., 1985) and can grow in laboratory culture at pH 5.2 (Eardly et al., 1992). Most P. vulgaris Rhizobium microsymbionts can also grow in laboratory culture at pH 5.2 (Eardly et al., 1992). In contrast, the E. meliloti strains that are the usual microsymbionts of M. sativa are known to be particularly acid-sensitive in comparison to other rhizobial species, and in laboratory culture fail to grow below pH 5.6 (Graham and Parker, 1964; Eardly et al., 1992; Garau et al., 2005). In order to gain insights into the acid tolerance of OR191, we examined the genome for the presence of loci known to be important for the pH adaptation of very acid tolerant strains (such as Rhizobium tropici CIAT899^T) and to less acid tolerant but well characterized strains (such as Ensifer meliloti 1021 and E. medicae WSM419). Although less acid-tolerant than CIAT899^T, WSM419 is the dominant nodule occupant of Medicago species growing in moderately acidic soils in Sardinia and Greece (Garau et al., 2005).

Most of the genes found to be important for acid tolerance in Rhizobium strains had orthologs in WSM419 and Ensifer meliloti 1021 (Supplementary Table 3), including the most acid-activated genes discovered so far, lpiA and acvB (Reeve et al., 2006), encoding lysyl-phosphatidylglycerol synthase and lysyl-phosphatidylglycerol hydrolase, respectively, and required for modulation of lysyl-phosphatidylglycerol homeostasis (Sohlenkamp et al., 2007). However, the fsrR, tcrA, and tcsA regulatory system, which is implicated in the strong acid induction of the lpiA-acvB operon in WSM419, is not present in OR191 or in species other than E. medicae, indicating that an unknown regulatory system governs acid induction of these genes in other rhizobia (Reeve et al., 2006; Tian et al., 2017). In OR191 and several other Rhizobium strains, the lpiA-acvB operon is immediately downstream of a gene encoding a small-conductance mechanosensitive channel, a membrane protein that is important for bacterial survival (reviewed in Naismith and Booth, 2012). The olsC gene, which has been identified as being important for the production of hydroxylated ornithine lipid species that increase stress tolerance in CIAT899^T (Rojas-Jiménez et al., 2005; Vences-Guzmán et al., 2011), was present in both the highly acid tolerant strain R. tropici CIAT899^T and in R. favelukesii strains but was absent from E. medicae WSM419 and E. meliloti 1021. This accords with proteomic and transcriptomic studies of R. favelukesii LPU83^T that revealed changes in proteins and genes associated with lipid metabolism involving the cell envelope in response to acid stress, however, these studies also suggested that some acid-responsive regulatory systems differ in LPU83^T compared with Ensifer strains (Nilsson et al., 2019, 2020).

R. favelukesii OR191 nodulates and fixes nitrogen with both M. sativa and P. vulgaris. To identify features within the OR191 genome that enable it to be a microsymbiont of both these hosts, the IMG online database was used to compare relevant sequenced genomes of microsymbionts of Medicago species, or of common bean, or of both (strains listed in Supplementary Table 4).

We examined these genomes for genes with known roles in rhizobial symbiotic interactions. In addition to meeting or exceeding the cut-off values of > 30% protein identity over > 75% protein coverage, OR191 nod, nif and fix (cco) homologs were validated as orthologous if they were present within the neighborhood of symbiotic gene clusters. A summary of the identified OR191 genes with characterized roles in symbiotic interactions and their phylogenetic profile in microsymbionts of Medicago and P. vulgaris is provided in Supplementary Table 5.

Symbiosis between legumes and rhizobia is initiated by a molecular dialogue, in which the host secretes flavonoid signals from the root cells, triggering expression of the rhizobial nodulation (nod, noe, and nol) genes and production of lipochitoligosaccharide Nod factors, which are required for infection of the host and nodule organogenesis. Expression of the nodulation genes is controlled by the regulatory protein NodD. In addition to the Nod factor core biosynthetic genes (nodABC) and export genes (nodIJ), rhizobial genomes contain accessory nodulation genes, which encode moieties that decorate the basic Nod factor and act as specific host determinants (Perret et al., 2000).

The core and accessory nodulation genes of OR191 are distributed on four scaffolds (Supplementary Figure 4) and share 99–100% sequence identity with those of R. favelukesii LPU83^T. While all rhizobial genomes contained the core nodulation genes nodABCDIJ, there were differences in the accessory nodulation genes present within Medicago and P. vulgaris microsymbionts.

In addition to the core nodABCIJ genes, all Medicago microsymbiont genomes (E. meliloti strains + OR191, LPU83^T, R. tibeticum CGMCC 1.7071^T and R. mongolense USDA 1844^T) contained a specific set of accessory nod genes. All Medicago microsymbiont genomes harbored three nodD alleles, which in E. meliloti appear to be important in optimizing symbiotic interactions with Medicago hosts (Honma et al., 1990; Van Berkum et al., 1998). The OR191 nodD genes are most closely related to orthologs in LPU83^T and the Medicago-nodulating Rhizobium strains R. tibeticum CGMCC 1.7071^T and R. mongolense USDA 1844^T. All Medicago microsymbiont genomes also contained the nodL and nodHPQ genes that code for synthesis of acetylated and sulfated Nod factors, respectively, which in Ensifer strains are required for symbiotic interaction with Medicago hosts (Lerouge et al., 1990). The Nod factors of R. favelukesii LPU83^T include species with trimeric, tetrameric, and pentameric chitin backbones that contain unsaturated acyl chains and may be sulfated or methylated, however, sulfated Nod factors do not appear to be required for this strain to nodulate alfalfa, as mutation of nodH delayed, but did not abolish, nodulation of this host (Del Papa et al., 2007; Torres Tejerizo et al., 2011). Both LPU83^T and OR191 contained the flavonoid-inducible nodL, noeA, and noeB operon, which is required for effective nodulation of Medicago sp. (Ardourel et al., 1995). NoeA and NoeB are hypothesized to modify the Nod factor for accurate recognition by Medicago spp., allowing the microsymbiont to successfully infect the host.

The nodEF genes, which are required for production of Nod factors with α, β-unsaturated acyl chains (Dénarié et al., 1996), were also present in all Medicago microsymbiont genomes. NodE has homology to β-ketoacyl synthases, while NodF is homologous to acyl carrier proteins (Demont et al., 1993). The nodG gene, encoding a 3-oxoacyl-acyl carrier protein reductase, was identified in all Medicago-nodulating E. meliloti genomes. NodG has a postulated role in fatty acid synthesis in the M. sativa microsymbiont strain E. meliloti 1021 (Mao et al., 2016), but is not required for synthesis of the characteristic α, β-unsaturated acyl chains of the Nod factors (Debellé et al., 2001). Although nodG is not essential for nodulation of M. sativa, it does increase nodulation efficiency (Mao et al., 2016). In OR191, the gene immediately upstream of nodP has sequence similarity to the characterized nodG of E. meliloti but encodes a truncated protein of 55 amino acids instead of the 245 amino acids found for E. meliloti NodG. The lack of a functional NodG might partially explain OR191’s less efficient nodulation performance on M. sativa, in comparison to 1021. R. favelukesii LPU83^T also contains a truncated version of nodG (Torres Tejerizo et al., 2011). We identified similarly truncated nodG genes in the same gene neighborhoods within the genomes of R. tibeticum CGMCC 1.7071^T and R. mongolense USDA 1844^T, indicating that nodG pseudogenes may be a feature of these strains. The nolFG genes, encoding an RND-type efflux pump, are also characteristic of Medicago-nodulating Ensifer strains (Baev et al., 1991) and disruption of this region leads to delayed nodulation of alfalfa (Baev et al., 1991).

The OR191 nodABC, nodEF, nodL, nodHPQ, and nolFG genes are most closely related to those of LPU83^T, CGMCC 1.7071^T and USDA 1844^T, and then to orthologs in Ensifer strains that nodulate Medicago. On the basis of its nod genes and symbiotic phenotype, Rogel et al. (2011) have included OR191 in the symbiovar orientale (Silva et al., 2005; Rogel et al., 2011), along with R. favelukesii LPU83^T, R. gallicum CCBAU 01083, R. mongolense USDA 1844^T, and R. tibeticum CGMCC 1.7071^T (Van Berkum et al., 1998; Zhao et al., 2008; Hou et al., 2009; Torres Tejerizo et al., 2016).

A comparison of OR191 to other P. vulgaris microsymbionts revealed that the only nod genes they have in common are the core nodABCDIJ genes. Both the types of nod genes and the gene arrangements varied considerably in the P. vulgaris microsymbionts; for example, from two to six nodD alleles and from one to three nodA alleles were found and nodHPQ and nodEF clusters were present in some genomes but not others. Similarly, nodU, nodZ, noeI, noeT, and nolL genes were variably present or absent. However, all the P. vulgaris microsymbionts except OR191, LPU83^T, R. tibeticum CGMCC 1.7071^T and R. mongolense USDA 1844^T contained nodS, which encodes an S-adenosyl-L-methionine (SAM)-dependent N-methyltransferase. Notably, nodS is absent from the genomes of all Medicago microsymbionts. A nodS-dependent N-methylation on the non-reducing end of the Nod Factor has been shown to be essential for the nodulation of P. vulgaris by the Rhizobium tropici strain CIAT 899^T (Jabbouri et al., 1995). Although the Medicago microsymbionts (including the R. favelukesii, R. tibeticum, and R. mongolense strains) lack nodS, they contain the nodL gene required for production of O-acetylated Nod factors. NodL acetylation of the Nod factor functionally prevents the NodS-dependent transfer of the N-methyl group substituent, however, it has been proposed that NodL-dependent O-acetylation of the Nod factor can compensate for the lack of this N-methyl group (Waelkens et al., 1995; Lopez-Lara et al., 2001). Nodulation of P. vulgaris would appear to require either an acetylated or a methylated Nod factor, which can be accomplished by NodL or NodS, respectively.

The nodEF genes were found in R. favelukesii, R. tibeticum, and R. mongolense, Ensifer fredii GR64, E. meliloti GVPV12, Ensifer sp. strains 4H41 and BR816 and strains closely related to R. tropici CIAT 899^T, but not in the remaining P. vulgaris-nodulating strains. Although the characterized strain R. tropici CIAT 899^T contains nodEF, it does not produce Nod factors with α,β-unsaturated acyl chains (Ormeño-Orrillo et al., 2012). This can be related to the inability of R. tropici NodA to transfer unsaturated acyl chains to the Nod factor, whereas NodA of E. meliloti Medicago microsymbionts specifies the N-acylation of the Nod factor by an unsaturated or hydroxylated fatty acid (Debellé et al., 1996). However, the CIAT 899^T hsnT, nodF, and nodE genes, which are in an operon with nodA2 in this strain, are important for the biosynthesis and specific decoration of the Nod Factors, with resulting impact on host specificity and symbiotic performance in some legume species (Gomes et al., 2019).

Because the transfer of unsaturated fatty acid chains to the Nod factor backbone requires a specific NodA protein (Debellé et al., 1996), we performed a comparative sequence analysis of the NodA proteins of the Medicago and P. vulgaris microsymbionts. The resulting phylogenetic tree grouped NodA into eight clusters (Group A—G, Figure 1) and showed that Cluster B (NodA of Medicago-nodulating R. favelukesii, R tibeticum, and R mongolense strains) formed a sister group to Cluster A of the Medicago-nodulating E. meliloti strains. This illustrates that even though E. meliloti 4H41 and GVPV12, Ensifer sp. BR816 and strains closely related to R. tropici CIAT 899^T contain nodEF and therefore are potentially able to produce α, β-unsaturated acyl chains, their NodA is phylogenetically distinct from the NodA of Medicago-nodulating strains that can transfer unsaturated acyl chains to the Nod factor.

By mapping nod gene neighborhoods to the NodA phylogenetic tree, we identified that particular arrangements of nod genes are associated with each NodA group (Figure 1). The figure illustrates the high degree of synteny between the nod genes of the Medicago-nodulating Ensifer strains (Group A) and the R. favelukesii, R tibeticum and R mongolense strains (Group B). R. tibeticum CGMCC 1.7071^T is the microsymbiont of Medicago ruthenica, while R. mongolense USDA 1844^T is the microsymbiont of Medicago archiducis-nicolai (Van Berkum et al., 1998; Hou et al., 2009). Both these Medicago species belong to the basal section Platycarpae and are distributed in Siberia, Mongolia, Tibet and northern China, which is far to the northeast of alfalfa’s center of origin (Small, 2010; Steele et al., 2010), suggesting that there may have been a change of microsymbiont in Medicago species, due to either biogeographical factors or a change in host preference. As noted previously, all P. vulgaris-nodulating strains (Group B—G) have either nodS or nodL in their genomes. In the strains belonging to Group D to G, nodSU is consistently located between nodC and nodIJ. The truncated nodU seen in the P. vulgaris-nodulating strains of Group C would appear to have arisen from an ancestral nodSU form. In the Medicago-nodulating strains belonging to Group A and B, the deletion of nodSU would account for nodABCIJ being adjacent to each other.

The rhizobial nif, fix and fdx genes are required for the stepwise reduction of dinitrogen gas and incorporation of fixed N (Fischer, 1994). Whereas the characterized rhizobial diazotroph Azorhizobium caulinodans ORS 571^T has 15 nif genes, most rhizobia are unable to fix nitrogen ex planta and their genomes contain fewer than 15 nif genes (Masson-Boivin et al., 2009). We identified 13 nif genes in OR191 (Supplementary Figure 4 and Supplementary Table 5). The genes and the associated processes that they are involved in are: nifA (regulation); nifH (dinitrogenase reductase); nifDK (α and β subunits of dinitrogenase); nifB (synthesis of the iron-sulfur-containing precursor of the FeMo-co); nifEN (assembly of the FeMo-co); nifQSX (FeMo-co biosynthesis), nifZ (nitrogenase maturation) and fixU/nifTnifW (nitrogen fixation accessory and stabilization proteins, respectively). OR191 lacked genes encoding NifU (iron-sulfur cluster scaffolding protein) and NifV (homocitrate synthase). The OR191 genome also contained genes required for production of ferredoxins (fdxBN), putative electron transfer proteins (fixABCX), the specific high-affinity cbb3-type cytochrome c oxidase (fixNOQP), a Cu²⁺-exporting ATPase (fixI), electron transfer components (fixGHS) and a putative FNR/CRP [fumarate and nitrate reductase regulator protein and cyclic AMP (cAMP) receptor protein] transcriptional regulator (fixK) (Supplementary Table 5). The OR191 nif and fix genes are most closely related to those of LPU83^T, CGMCC 1.7071^T and USDA 1844^T and then to those of Rhizobium and Ensifer strains that nodulate P. vulgaris, rather than to the nif and fix genes of Medicago-nodulating Ensifer strain. Phylogenetic analysis of NifH and NifA shows that OR191, LPU83^T, CGMCC 1.7071^T, and USDA 1844^T form a clade in comparison to other microsymbiont strains analyzed in this study (Figure 2A).

Furthermore, a comparison of the nif genes present in the microsymbionts of Medicago spp. and P. vulgaris showed that all genomes contained the nitrogenase genes nifABHDKENX and fixU/nifT or nifT-like genes (encoding the specific NifT pfam06988), fdxBN and fixABCX genes (Figure 2B). All strains contained fixGHIS and fixNOQP, apart from the Paraburkholderia microsymbionts (previously reported in De Meyer et al., 2016). The most striking differences between the suite of nif genes found in the Medicago and P. vulgaris rhizobia was that Ensifer strains that are Medicago microsymbionts all lacked nifQWZSUV, whereas strains that nodulate and fix with P. vulgaris all contained nifQWZS. All the Rhizobium strains contained nifS but lacked nifU and nifV. Ensifer strains that were microsymbionts of P. vulgaris (E. fredii GR64, E meliloti GVPV12 and Ensifer sp. 4H41 and BR816) all contained nifU and nifV, in addition to the 13 nif genes identified in OR191. The nifU gene was absent from the Paraburkholderia strains, but as rhizobial Paraburkholderia strains are able to fix N ex planta (Elliott et al., 2007), nifU may be replaced by an unidentified gene encoding a functionally equivalent protein.