Symbiotic and Nonsymbiotic Members of the Genus Ensifer (syn. Sinorhizobium) Are Separated into Two Clades Based on Comparative Genomics and High-Throughput Phenotyping.

Rhizobium-legume symbioses serve as paradigmatic examples for the study of mutualism evolution. The genus Ensifer (syn. Sinorhizobium) contains diverse plant-associated bacteria, a subset of which can fix nitrogen in symbiosis with legumes. To gain insights into the evolution of symbiotic nitrogen fixation (SNF), and interkingdom mutualisms more generally, we performed extensive phenotypic, genomic, and phylogenetic analyses of the genus Ensifer. The data suggest that SNF likely emerged several times within the genus Ensifer through independent horizontal gene transfer events. Yet, the majority (105 of 106) of the Ensifer strains with the nodABC and nifHDK nodulation and nitrogen fixation genes were found within a single, monophyletic clade. Comparative genomics highlighted several differences between the "symbiotic" and "nonsymbiotic" clades, including divergences in their pangenome content. Additionally, strains of the symbiotic clade carried 325 fewer genes, on average, and appeared to have fewer rRNA operons than strains of the nonsymbiotic clade. Initial characterization of a subset of ten Ensifer strains identified several putative phenotypic differences between the clades. Tested strains of the nonsymbiotic clade could catabolize 25% more carbon sources, on average, than strains of the symbiotic clade, and they were better able to grow in LB medium and tolerate alkaline conditions. On the other hand, the tested strains of the symbiotic clade were better able to tolerate heat stress and acidic conditions. We suggest that these data support the division of the genus Ensifer into two main subgroups, as well as the hypothesis that pre-existing genetic features are required to facilitate the evolution of SNF in bacteria.

Significance

The bacterial genus Ensifer contains ecologically important N2-fixing symbionts of leguminous plants, as well as nonsymbiotic species. However, the evolutionary dynamics of symbiotic nitrogen fixation within this genus are unclear, and it remains an open question of whether the gain of classical symbiotic N2-fixation genes is sufficient to allow a bacterium to fix nitrogen. Our results suggest that the symbiotic species of the genus Ensifer predominately group separately from the nonsymbiotic species, but that symbiotic abilities were likely acquired multiple times within this group. This study provides new insight into the evolution of symbiotic N2-fixation in a bacterial genus, while supporting the hypothesis that genetic features aside from the classical symbiotic N2-fixation genes contribute to the evolution of symbiotic potential.

Introduction

Symbioses are pervasive phenomena present in all Eukaryotic forms of life (López-García et al. 2017). These includes facultative symbiotic interactions, obligate symbioses, and the evolution of organelles (Douglas 2014), with symbiotic nitrogen fixation (SNF) being a paradigmatic example of the latter (Masson-Boivin and Sachs 2018). SNF (the conversion of N2 to NH3) is performed by a polyphyletic group of bacteria from the Alphaproteobacteria and Betaproteobacteria (whose nitrogen-fixing members are collectively called rhizobia) and members of the genus Frankia (Masson-Boivin et al. 2009; Wang and Young 2019) while intracellularly housed within specialized organs (nodules) of specific plants in the family Fabaceae and the genus Parasponia, as well as the actinorhizal plants (Werner et al. 2014; Griesmann et al. 2018; van Velzen et al. 2018). The advantages and evolutionary constraints to SNF have long been investigated in the conceptual framework of mutualistic interactions and the exchange of goods (see, for instance, Heath and Tiffin 2007; Werner et al. 2015; Sørensen et al. 2019), and quantitative estimations with metabolic reconstructions have also been performed (Pfau et al. 2018; diCenzo et al. 2020).

The establishment of a symbiotic nitrogen-fixing interaction requires that the bacterium encode several diverse molecular functions, including those related to signaling and metabolic exchange with the host plant, nitrogenase and nitrogenase-related functions, and escaping or resisting the plant immune system (Oldroyd et al. 2011; Haag et al. 2013; Poole et al. 2018). In general, the primary genes required for SNF (i.e., the nod, nif, and fix genes) are located within mobile genetic elements that include symbiotic islands and symbiotic (mega)-plasmids (Checcucci et al. 2019; Tian and Young 2019; Geddes et al. 2020), facilitating their spread through horizontal gene transfer (HGT) (Sullivan et al. 1995; Barcellos et al. 2007; Pérez Carrascal et al. 2016). Emphasizing the role of HGT in the evolution of rhizobia, rhizobia are dispersed across seven families of the Alphaproteobacteria and one family of the Betaproteobacteria, and most genera with rhizobia also contain nonrhizobia (Garrido-Oter et al. 2018; Wang 2019).

An interesting area of investigation is whether the evolution of mutualistic symbioses, such as SNF, depends on metabolic/genetic requirements (“facilitators,” as in Gerhart and Kirschner [2007]) aside from the strict symbiotic genes (Long 2001; Sanjuán 2016; Zhao et al. 2018). In other words, 1) is the acquisition of symbiotic genes present in genomic islands or plasmids sufficient to become a symbiont? or, 2) are metabolic prerequirements or adaptation successive to HGT required? A comparative genomics study of 1,314 Rhizobiales genomes identified no functional difference between rhizobia and nonrhizobia based on Kyoto Encyclopedia of Genes and Genomes annotations (Garrido-Oter et al. 2018), suggestive of an absence of obvious facilitators. In contrast, experimental studies are generally consistent with an important role of nonsymbiotic genes in establishing or optimizing rhizobium–legume symbioses. Several studies have shown that effective symbionts are not produced following the transfer of symbiotic plasmids from rhizobia of the genera Rhizobium or Ensifer (syn. Sinorhizobium) to closely related nonrhizobia from the genera Agrobacterium or Ensifer (see, for instance, Hooykaas et al. 1982; Finan et al. 1986; Rogel et al. 2001; reviewed in diCenzo et al. [2019]). Similarly, the same symbiotic island is associated with vastly different symbiotic phenotypes depending on the Mesorhizobium genotype (Nandasena et al. 2007; Haskett et al. 2016). Further supporting the need for additional adaptations to support SNF, symbiosis plasmid transfer coupled to experimental evolution can lead to the gain of more advanced symbiotic phenotypes (Doin de Moura et al. 2020).

The genus Ensifer provides an ideal model to further explore the differentiation, or lack thereof, of symbiotic bacteria from nonsymbionts. This genus comprises rhizobia such as Ensifer meliloti and Ensifer fredii, as well as nonrhizobia like Ensifer morelense and Ensifer adhaerens, and many members have been extensively studied producing an abundant set of experimental and genomic data (for a recent review, see diCenzo et al. [2019]). The genus Ensifer, as currently defined, resulted from the combination of the genera Sinorhizobium and Ensifer based on similarities in the 16S rRNA and recA sequences of the type strains and the priority of the name Ensifer (Willems et al. 2003; Young 2003). Multilocus sequence analysis supported the amalgamation of these genera (Martens et al. 2007), although it was subsequently noted that E. adhaerens (the type strain) is an outgroup of this taxon based on whole genome phylogenomics (Ormeño-Orrillo et al. 2015). A more recent taxonomy approach based on genome phylogeny suggests that the genus Ensifer should again be split, with the initial type strains of Ensifer and Sinorhizobium belonging to separate genera (Parks et al. 2018).

In this article, we report an extensive comparative genomic analysis, and initial phenotypic characterization, of legume symbionts and nonsymbionts of the genus Ensifer. We identified that SNF likely evolved multiple times through independent HGT events; even so, most symbionts were found in a single clade, consistent with a requirement for pre-existing genetic features to facilitate the evolution of SNF. Moreover, the symbiotic and nonsymbiotic clades differed in their pangenome composition, and tests with a subset of strains suggested they also differed in their substrate utilization and resistance phenotypes as measured by the Phenotype MicroArray platform. We suggest that the data support the division of the genus Ensifer into two subgroups, corresponding to the genera Ensifer and Sinorhizobium of the Genome Taxonomy Database (Parks et al. 2018).

Materials and Methods

Genome Sequencing, Assembly, and Annotation

Prior to short-read sequencing, all strains were grown to stationary phase at 30 °C in TY medium (5 g l−1 tryptone, 3 g l−1 yeast extract, and 0.4 g l−1 CaCl2). Total genomic DNA was isolated using a standard cetyltrimethylammonium bromide method (Perrin et al. 2015). Short-read sequencing was performed at IGATech (Udine, Italy) using an Illumina HiSeq2500 instrument with 125-bp paired-end reads. Two independent sequencing runs were performed for E. morelense Lc04 and Ensifer psoraleae CCBAU 65732, whereas E. morelense Lc18 and Ensifer sesbaniae CCBAU 65729 were sequenced once. For long-read sequencing, E. sesbaniae was grown to midexponential phase at 30 °C in MM9 minimal medium (MOPS buffer [40 mM MOPS, 20 mM KOH], 19.2 mM NH4Cl, 8.76 mM NaCl, 2 mM KH2PO4, 1 mM MgSO4, 0.25 mM CaCl2, 1 µg ml−1 biotin, 42 nM CoCl2, 38 µM FeCl3, 10 µM thiamine-HCl, and 10 mM sucrose). Total genomic DNA was isolated as described elsewhere (Cowie et al. 2006). Long-read sequencing was performed in-house with a Pacific Biosciences Sequel instrument.

Reads were assembled into scaffolds using SPAdes 3.9.0 (Bankevich et al. 2012; Vasilinetc et al. 2015); in the case of E. sesbaniae, long reads were corrected and trimmed using Canu 1.7.1 (Koren et al. 2017) prior to assembly. Scaffolds returned by SPAdes were parsed to remove those with <20× coverage or with a length <200 nucleotides. Using FastANI (Jain et al. 2018), one-way average nucleotide identity (ANI) values of each assembly were calculated against 887 alpha-proteobacterial genomes available through the National Center for Biotechnological Information (NCBI) with an assembly level of complete or chromosome. Based on the FastANI output, each draft genome assembly was further scaffolded using MeDuSa (Bosi et al. 2015) and the reference genomes listed in supplementary table S1, Supplementary Material online. For most assemblies, scaffolds under 1 kb in length were discarded. The exception was for S. sesbaniae, for which case scaffolds <10 kb were discarded. Genome assemblies were annotated using Prokka 1.12-beta (Seemann 2014), annotating coding regions with Prodigal (Hyatt et al. 2010), tRNA with Aragon (Laslett and Canbäck 2004), rRNA with Barrnap (https://github.com/tseemann/barrnap; last accessed October 18, 2020), and ncRNA with Infernal (Kolbe and Eddy 2011) and Rfam (Kalvari et al. 2018).

Species Phylogenetic Analyses

All Ensifer (and Sinorhizobium) genomes were downloaded from the NCBI Genome Database regardless of assembly level. Strains that either 1) lacked a RefSeq assembly, 2) had genome sizes <1 Gb, or 3) appeared to not belong to the Ensifer clade based on preliminary phylogenetic analyses were discarded, leaving a final set of 157 strains (supplementary data set S1, Supplementary Material online). Eight complete Rhizobium genomes (supplementary data set S2, Supplementary Material online) were downloaded to serve as an outgroup. Genomes were reannotated with prokka to ensure consistent annotation. All genomes were downloaded on November 12, 2018, and associated metadata are available as supplementary data sets S1 and S2, Supplementary Material online.

To construct an unrooted, core gene phylogeny, the pangenome of the 157 Ensifer strains was calculated using Roary 3.11.3 (Page et al. 2015) with a percent identify threshold of 70%. As part of the running of Roary, the nucleotide sequences of the 1,049 core genes (identified as those found in at least 99% of the genomes; supplementary data set S3, Supplementary Material online) were individually aligned with PRANK (Löytynoja 2014) and the alignments concatenated. The concatenated alignment was trimmed using TRIMAL 1.2rev59 (Capella-Gutiérrez et al. 2009) with the automated1 option, and used to construct a maximum likelihood phylogeny (the bootstrap best tree following 100 bootstrap replicates, as determined by the extended majority-rule consensus tree criterion) using RAxML 8.2.9 (Stamatakis 2014) with the GTRCAT model as recommended (https://cme.h-its.org/exelixis/resource/download/NewManual.pdf; last accessed October 18, 2020). All phylogenies prepared in this study were visualized with the online iTOL webserver (Letunic and Bork 2016).

To construct a rooted phylogeny, the AMPHORA2 pipeline (Wu and Scott 2012) was used to identify 31 highly conserved bacterial proteins in each Ensifer and Rhizobium proteome, based on HMMER 3.1b2 (Eddy 2009) and the 31 hidden Markov models (HMMs) that come with AMPHORA2. Customs Perl scripts were then used to remove proteins that were either found in <95% of genomes or were found in multicopy in at least one genome, leaving a set of 30 proteins (Frr, InfC, NusA, Pgk, PyrG, RplA, RplB, RplC, RplD, RplE, RplF, RplK, RplL, RplM, RplN, RplP, RplS, RplT, RpmA, RpoB, RpsB, RpsC, RpsE, RpsI, RpsJ, RpsK, RpsM, RpsS, SmpB, Tsf). Orthologous groups were aligned using MAFFT 7.310 (Katoh and Standley 2013) with the localpair option, following which the alignments were trimmed using TRIMAL 1.2rev59 with the automated1 option. Alignments were concatenated and used to construct a maximum likelihood phylogeny (the bootstrap best tree following 304 bootstrap replicates, as determined by the extended majority-rule consensus tree criterion) using RAxML with the PROTGAMMAJTTDCMUT model. This model was chosen as a preliminary run using RAxML with the automatic model selection indicated that the best scoring tree was obtained with the selected model.

ANI and AAI Calculations

Pairwise ANI values were calculated for all Ensifer strains using FastANI (Jain et al. 2018) with default parameters; a value of 78% was used in cases where no value was returned by FastANI. Pairwise average amino acid identity (AAI) values were calculated with the compareM workflow (https://github.com/dparks1134/CompareM; last accessed October 18, 2020). Results were visualized and clustered using the heatmap.2 function of the gplots package in R (Warnes et al. 2016), with average linkage and Pearson correlation distances.

Pangenome Calculation

All proteins of the reannotated Ensifer strains were clustered into orthologous groups using CD-HIT 4.6 (Li and Godzik 2006) with a percent identity threshold of 70% and an alignment length of 80% of the longer protein. The output was used to determine core and accessory genomes using a prevalence threshold of 90% as many of the genomes were draft genomes. Gene accumulation curves were produced using the specaccum function of the vegan package of R (Oksanen et al. 2018), with the random method and 500 permutations. Principal component analysis (PCA) was performed with the prcomp function of R, and was visualized with the autoplot function the ggplot2 package (Wickham 2016).

Identification and Phylogenetic Analysis of Common Nod, Nif, and Rep Proteins

The proteomes were collected for the 157 Ensifer strains, as well as all strains from the genera Rhizobium, Neorhizobium, Agrobacterium, Mesorhizobium, and Ochrobactrum with an assembly status of Complete or Chromosome (supplementary data set S4, Supplementary Material online). Additionally, the seed alignments for the HMMs of the nodulation proteins NodA (TIGR04245), NodB (TIGR04243), and NodC (TIGR04242), the nitrogenase proteins NifH (TIGR01287), NifD (TIGR01282), NifK (TIGR01286), and the replicon partitioning proteins RepA (TIGR03453), and RepB (TIGR03454) were downloaded from TIGRFAM (Haft et al. 2012). Seed alignments were converted into HMMs with the HMMBUILD function of HMMER 3.1b2 (Eddy 2009). Each HMM was searched against the complete set of proteins from all 157 reannotated Sinorhizobium and Ensifer strains using the HMMSEARCH function of HMMER. The amino acid sequences for each hit (regardless of e-value) were collected. Each set of sequences was searched against a HMM database containing all 21,200 HMMs from the Pfam (Finn et al. 2016) and TIGRFAM databases using the HMMSCAN function of HMMER, and the top scoring HMM hit for each query protein was identified. Proteins were annotated as NodA, NodB, NodC, NifH, NifD, NifK, RepA, or RepB according to supplementary table S2, Supplementary Material online.

The nodA, nodB, and nodC genes are generally found as an operon. Thus, the NodA, NodB, and NodC proteins were putatively associated to operons based on identifying proteins that are encoded by adjacent genes in their respective genomes; orphan proteins not encoded by adjacent genes were discarded as the subsequent phylogenetic analysis was based on concatenated NodA, NodB, and NodC alignments. Each set of orthologs were aligned using MAFFT with the localpair option, and alignments trimmed using TRIMAL and the automated1 algorithm. Alignments were concatenated so as to combine alignments for proteins encoded by adjacent genes, producing a NodABC alignment. The same procedure was followed to produce NifHDK and RepAB alignments. Maximum likelihood phylogenies were built on the basis of each combined alignment using RAxML with the PROTGAMMAJTT (NodABC, NifHDK) or the PROTGAMMALG (RepAB) models. These models were chosen as preliminary runs using RAxML with the automatic model selection indicated that the best scoring trees were obtained with the selected models. The final phylogenies are the bootstrap best trees following 352 bootstrap replicates, as determined by the extended majority-rule consensus tree criterion.

Plant Assays

Phaseolus vulgaris (var. TopCrop, Mangani Sementi, Italy) seeds were surface sterilized in 2.5% HgCl2 solution for 2 min and washed five times with sterile water. Seeds were germinated in the dark at 23 °C, following which seedlings were placed in sterile polypropylene jars containing vermiculite:perlite (1:1) and nitrogen-free Fåhraeus medium, and grown at 23 °C with a 12-h photoperiod (100 µE m−2 s−1). One-week-old plantlets were inoculated with 100 µl of the appropriate rhizobium strain (suspended in 0.9% NaCl at an OD600 of 1); five plants were inoculated per strain and then grown for 4 weeks at 23 °C with a 12-h photoperiod (100 µE m−2 s−1). Plant growth assays were repeated three independent times. Nodules were collected and surface sterilized as described elsewhere (Checcucci et al. 2016), crushed in sterile 0.9% NaCl solution, and serial dilutions were plated on TY agar plates and incubated at 30 °C for 2 days. PCR amplification of the 16S rRNA gene was performed using crude lysates from single colonies recovered from root nodules, as in Barzanti et al. (2007). Sequencing of the PCR amplified 16S rRNA gene was performed from both the 27f and the 1495r primers using BrilliantDye Terminator Cycle Sequencing chemistry (Nimagen, Nijmegen, The Netherlands) on a 3730xl DNA Analyzer (ThermoFisher Scientific, Waltham, MA).

Phenotype MicroArray

Phenotype MicroArray experiments using Biolog plates PM1 and PM2A (carbon sources), PM9 (osmolytes), and PM10 (pH) were performed as described previously (Biondi et al. 2009). Data were collected over 96 h with an OmniLog instrument. Data analysis was performed with DuctApe (Galardini et al. 2014). Activity index (AV) values were calculated following subtraction of the blank well from the experimental wells. Growth with each compound was evaluated with AV values from 0 (no growth) to 9 (maximal growth), following an elbow test calculation. Phenotype MicroArray experiments were performed once as results for these experiments are highly repeatable (Johnson et al. 2008; Bochner et al. 2010; Dunkley et al. 2019).

Biofilm Assays

Overnight cultures of strains grown in TY and LB (10 g l−1 tryptone, 5 g l−1 yeast extract, 5 g l−1 NaCl) media were diluted to an OD600 of 0.02 in fresh media, and six replicates of 100 µl aliquots were transferred to a 96-well microplate. Plates were incubated at 30 °C for 24 h, after which the OD600 was measured with a Tecan Infinite 200 PRO (Switzerland). Each well was then stained with 30 µl of a filtered 0.1% (w/v) crystal violet solution for 10 min, and then the medium containing the planktonic cells was gently removed from the wells. Next, the wells were rinsed three times with 200 μl of phosphate-buffered saline (PBS; 0.1 M, pH 7.4) and allowed to dry for 15 min. About 100 µl of 95% (v/v) ethanol was added to each well and then incubated for 15 min at room temperature. The OD540 of each well was measured (Rinaudi and González 2009), and biofilm production reported as the ratio of the OD540/OD600 ratio. Biofilm assays consisted of six replicates, and were performed two independent times.

Growth Curves

Overnight cultures of each strain were grown in the same medium to be used for the growth curve. For minimal media, either 0.2% (w/v) of glucose or succinate was added as the carbon source. Cultures were diluted to an OD600 of 0.05 in the same media, and triplicate 150-µl aliquots were added to a 96-well microplate. Microplates were incubated without shaking at 30 °C or 37 °C in a Tecan Infinite 200 PRO, with OD600 readings taken every hour for 48 h. Growth rates were evaluated over 2-h windows during the exponential growth phase. All growth curves were performed in triplicate and repeated two independent times.

To evaluate bacterial growth when provided root exudates as a nitrogen source, root exudates were produced from Medicago sativa cv. Maraviglia as described elsewhere (Checcucci et al. 2017). Single bacterial colonies from TY plates were resuspended in a 0.9% NaCl solution to an OD600 of 0.5. Then, each well of a 96-well microplate was inoculated with 5 µl of culture, 75 µl of nitrogen-free M9 with 0.2% (w/v) succinate as a carbon source, and 20 µl of root exudate as a nitrogen source as done previously (Checcucci et al. 2017). Triplicates were performed for each strain. Microplates were incubated without shaking at 30 °C in a Tecan Infinite 200 PRO, with OD600 readings taken every hour for 48 h. Growth rates were determined as described above.

Results

Genome Sequencing of Four Rhizobiaceae Strains

Draft genomes of E. morelense Lc04, E. morelense Lc18, E. sesbaniae CCBAU 65729, and E. psoraleae CCBAU 65732 (Wang et al. 1999, 2013) were generated to increase the species diversity available for our analyses. Summary statistics of the assemblies are provided in supplementary table S3, Supplementary Material online. The genome sequences confirmed the presence of nodulation and nitrogen-fixing genes in E. morelense Lc18, E. sesbaniae CCBAU 65729, and E. psoraleae CCBAU 65732, whereas these genes appeared absent in the E. morelense Lc04 assembly. Strains Lc04, CCBAU 65729, and CCBAU 65732 were confirmed to belong to the genus Ensifer, as one-way ANI comparisons revealed that the most similar alpha-proteobacterial genomes were from the genus Ensifer. However, the genome of strain Lc18 was most similar to genomes from the genera Rhizobium and Agrobacterium, consistent with an earlier 16S rRNA gene restriction fragment length polymorphism analysis (34). Thus, we propose renaming E. morelense Lc18 to Rhizobium sp. Lc18. As this strain does not belong to the genus Ensifer, it was excluded from further analyses.

Symbiotic and Nonsymbiotic Ensifer Strains Segregate Phylogenetically

An unrooted, core gene phylogeny of 157 Ensifer strains was prepared to evaluate the phylogenetic relationships between the symbiotic and nonsymbiotic strains (fig. 1). A rooted phylogeny based on a multilocus sequence analysis was also prepared (supplementary fig. S1, Supplementary Material online). Each of the 157 strains was annotated as symbiotic or nonsymbiotic based on the presence of the common nodABC nodulation genes and the nifHDK nitrogenase genes. Consistent with previous work (Garrido-Oter et al. 2018), both phylogenies revealed a clear division of the symbiotic and nonsymbiotic strains into two well-defined clades. Nevertheless, a few exceptions were noted. Ensifer sesbaniae was found within the nonsymbiotic clade; however, E. sesbaniae was reported to be a symbiont of legumes such as P. vulgaris (Wang et al. 2013), and the ability of E. sesbaniae to nodulate P. vulgaris was confirmed in this study (supplementary fig. S2, Supplementary Material online). Similarly, at least one of the six symbiotic proteins were not detected in five strains of the symbiotic group, although we cannot rule out that these are false negatives due to incomplete genome assemblies or genome assembly errors. ANI (genospecies threshold: 95%) and AAI (genospecies threshold: 96%) calculations suggested the presence of 12 and 20 genospecies within the nonsymbiotic and symbiotic groups, respectively (fig. 1 and supplementary figs. S3 and S4, Supplementary Material online), confirming that the nonsymbiotic clade is not an artifact of low species diversity. Thus, we conclude that the genus Ensifer consists of two well-defined clusters, each consisting predominately of either symbiotic or nonsymbiotic strains.

. 1.

Unrooted phylogeny of the genus Ensifer. A maximum likelihood phylogeny of 157 strains was prepared from a concatenated alignment of 1,049 core genes. Nodes with a bootstrap value of 100 are indicated with the gray dots. The scale represents the mean number of nucleotide substitutions per site. The white and gray shading is used to group strains into genospecies on the basis of ANI and AAI results (supplementary figs. S3 and S4, Supplementary Material online), using ANI and AAI genospecies threshold of 95% and 96%, respectively. From outside to inside, rings represent the genome assembly level (black, finished; white, draft), and the presence (black) or absence (white) of NodA, NodB, NodC, NifH, NifD, and NifK. Gray boxes indicate the presence of a truncated version of the corresponding gene (as a result of incomplete genome assembly) detected through inspection of the RefSeq annotations. Strains are named as recorded in NCBI at the time of collection.

SNF Likely Arose Multiple Times within the Genus Ensifer

A possible explanation for the phylogenetic segregation of SNF within the genus Ensifer was that the symbiotic genes were gained once through a single HGT event. To test this hypothesis, the phylogenetic relationships of the NodABC and NifHDK proteins of the order Rhizobiales were examined (fig. 2). SNF genes are situated on megaplasmids in the genus Ensifer; thus, a phylogeny of RepAB partitioning proteins of the order Rhizobiales was prepared as a proxy of the evolutionary relationships among the symbiotic megaplasmids (fig. 2). We predicted that the NodABC, NifHDK, and RepAB proteins of the genus Ensifer would form a single, monophyletic clade in each phylogeny, if the above hypothesis were true. This was not observed. Instead, all three phylogenies were inconsistent with a single origin of SNF within the genus Ensifer as the Ensifer strains were predominantly split into three clades: 1) E. meliloti and Ensifer medicae, 2) E. fredii and related strains, and 3) Ensifer americanum and related strains. As the same clades are observed in the species tree (fig. 2), this observation suggests SNF was independently acquired through HGT in each clade. The relationships between the SNF genes of the remaining Ensifer species (e.g., Ensifer aridi and E. psoraleae) was not clear; however, the most parsimonious solution was that there were three additional acquisitions of SNF (fig. 2). Overall, the phylogenetic analyses of the NodABC, NifHDK, and RepAB proteins support the hypothesis that there were multiple, independent acquisition of symbiosis genes (hence SNF) by lineages within the genus Ensifer; however, the gain (and/or maintenance) of symbiosis genes preferentially occurred within one monophyletic group of species.

. 2.

Evolution of SNF within the genus Ensifer. Maximum likelihood phylogenies of concatenated alignments of (A) NodABC nodulation proteins, (B) NifHDK nitrogenase proteins, and (C) RepAB replicon partitioning proteins of the order Rhizobiales. Branches corresponding to proteins from the genus Ensifer are indicated with color. (D) A subtree of the core gene species phylogeny of figure 1. Colors denote taxa whose symbiotic proteins are predicted to have been vertically acquired from a common ancestor. The scale bars represent the mean number of amino acid (A–C) or nucleotide (D) substitutions per site.

The Genomic Features of the Symbiotic and Nonsymbiotic Clades Differ

The pangenome of the 157 Ensifer strains was calculated to evaluate if there were global genomic differences between the symbiotic and nonsymbiotic clades. Both clades had open pangenomes (supplementary fig. S5, Supplementary Material online). A PCA based on gene presence/absence revealed a clear separation of the two clades (fig. 3), suggesting a divergence of the pangenomes of these clades. The symbiotic clade was subdivided into two groups along the second component of the PCA (fig. 3), which may suggest further levels of genomic separation. About 2,130 genes were found in the core genomes of both clades, whereas 20% (542 genes) and 40% (1,377 genes) of the core genomes of the symbiotic and nonsymbiotic clades, respectively, were absent from the core genome of the other clade; of these, about a third were completely absent from the other clade’s pangenome (fig. 3). Of the 14,514 accessory genes (defined as genes found in at least 10% of at least one clade, excluding the 2,130 Ensifer core genes), only 2,352 (16%) were found in the pangenomes of both the symbiotic and the nonsymbiotic clades. Moreover, a statistically significant difference (Wilcoxon rank-sum test, P < 0.0001) in the genome sizes of the two clades was observed (fig. 3); strains of the nonsymbiotic clade carried 325 more genes, on average, than strains of the symbiotic clade (median difference of 470). Finally, based on the limited number of strains with finished genomes, strains of the symbiotic clade appear to generally have three copies of the rRNA operon whereas strains of the nonsymbiotic clade appear to have a norm of five copies of their rRNA operon. Together, these multiple lines of data are consistent with there being a broad genomic divergence of the symbiotic and nonsymbiotic clades of the genus Ensifer.

. 3.

Global genome properties of the genus Ensifer. (A) A PCA plot based on the presence and absence of all orthologous protein groups in each of the 157 Ensifer strains. (B) Box-and-whisker plots displaying the number of genes per genome in the symbiotic and nonsymbiotic Ensifer clades. (C) A Venn diagram displaying the overlap in the core genomes of the symbiotic and nonsymbiotic Ensifer clades. (D) A Venn diagram displaying the overlap in the accessory genomes of the symbiotic and nonsymbiotic Ensifer clades.

Phenotypic Features of the Symbiotic and Nonsymbiotic Clades Differ

A subset of ten strains (table 1), five each from the symbiotic and nonsymbiotic clades, were subjected to a panel of assays to investigate how phenotypes vary across the genus Ensifer. These ten strains were chosen so as to provide broad phylogenetic coverage of the genus, while excluding strains for which extensive phenotypic characterizations have been previously published. No statistically significant differences were observed in the ability of members of the two clades to form biofilms (supplementary fig. S6 and table S4, Supplementary Material online). However, the tested strains clearly differed in their ability to grow in LB media; whereas the tested strains of the nonsymbiotic clade displayed robust growth in LB, the tested strains of the symbiotic clade largely failed to grow (fig. 4). Tested strains of the nonsymbiotic clade also displayed a slightly faster specific growth rate, on average, than the tested strains of the symbiotic clade in TY media (nonsymbiotic clade: 0.54 ± 0.03 h−1; symbiotic clade: 0.44 ± 0.08 h−1; P = 0.03 from an ANOVA followed by Tukey’s test; supplementary fig. S7A and table S5, Supplementary Material online). On the other hand, tested strains of the symbiotic clade were, on average, better able to withstand heat stress (37 °C) in TY media (fig. 4). No statistically significant difference in the average ability of the tested strains of the two clades to grow in minimal media with succinate or glucose as a carbon source, or with M. sativa (a symbiotic partner of E. meliloti and E. medicae) root exudates as a nitrogen source, was detected (supplementary fig. S7 and table S5, Supplementary Material online).

Table 1

Ensifer Strains Phenotypically Characterized in This Study

Strain	Original Source	SNFa	Ensifer Cladeb	Reference
Ensifer adhaerens Casida A	Isolated from a Pennsylvania soil sample	No	Nonsymbiotic	Casida (1982)
Ensifer adhaerens OV14	Isolated from the rhizosphere of Brassica napus	No	Nonsymbiotic	Wendt et al. (2012)
Ensifer sp. M14	Isolated from arsenic-rich sediments of a gold mine	No	Nonsymbiotic	Drewniak et al. (2008)
Ensifer morelense Lc04	Isolated from root nodules of Leucaena leucocephala	No	Nonsymbiotic	Wang et al. (1999)
Ensifer sesbaniae CCBAU 65729	Isolated from root nodules of Sesbania cannabina	Yes	Nonsymbiotic	Wang et al. (2013)
Ensifer fredii NGR234	Isolated from root nodules of Lablab purpureus	Yes	Symbiotic	Trinick (1980)
Ensifer sojae CCBAU 05684	Isolated from root nodules of Glycine max grown in saline-alkaline soils	Yes	Symbiotic	Li et al. (2011)
Ensifer americanum CFNEI 156	Isolated from root nodules of Acacia acatlensis	Yes	Symbiotic	Toledo et al. (2003)
Ensifer psoraleae CCBAU 65732	Isolated from root nodules of Psoralea corylifolia	Yes	Symbiotic	Wang et al. (2013)
Ensifer medicae WSM419	Isolated from root nodules of Medicago murex	Yes	Symbiotic	Howieson and Ewing (1986)

This column indicates if the strain can (Yes) or cannot (No) form nitrogen-fixing nodules on legumes.

This column indicates if the strain belongs to the symbiotic or nonsymbiotic clade of the genus Ensifer as defined in this study.

. 4.

Growth properties of phylogenetically diverse Ensifer strains. Ensifer strains were grown in microplates without shaking. Data points represent the average of triplicate samples, whereas the error bars indicate the SD. Shades of pink are used for strains of the symbiotic clade, whereas shades of blue are used for strains of the nonsymbiotic clade. (A) Growth in LB medium at 30 °C. (B) Growth in TY medium during heat stress (37 °C).

The phenotypic properties of the genus Ensifer were further examined through evaluating the ability of the same ten strains to catabolize 190 carbon sources, and to grow in 96 osmolyte and 96 pH conditions, through the use of Biolog Phenotype MicroArrays. Clustering the strains based on growth properties largely separated the tested strains of the symbiotic clade and nonsymbiotic clade into distinct groups (fig. 5 and supplementary fig. S8, Supplementary Material online). The exception was Ensifer sojae, which formed its own intermediate group in the phenotype data. To aid in identifying which conditions best separate the tested strains of the symbiotic clade (including E. sojae) from the tested strains of the nonsymbiotic clade, a linear discriminant analysis (LDA) was run over the AV values summarizing growth in each condition (supplementary data set S5, Supplementary Material online). In general, tested strains from the nonsymbiotic clade better tolerated high pH (pH 9.0–9.5) than did the tested strains from the symbiotic cluster. In contrast, tested strains of the symbiotic clade had better tolerance to low pH conditions (pH 3.5–4.5). In addition, the tested strains from the two clades clearly differed in their overall metabolic abilities with tested strains of the nonsymbiotic clade generally having a broader metabolic capacity than those of the symbiotic clade (supplementary table S6, Supplementary Material online). Whereas tested strains of the symbiotic clade displayed robust growth on 65 carbon sources on average, the tested strains of the nonsymbiotic clade grew on an average of 81 carbon sources (P < 0.05, Student’s t-test). Overall, these initial experiments provide support for the hypothesis that a variety of phenotypes, not just the ability to nodulate legumes, differ between the symbiotic and nonsymbiotic clades of the genus Ensifer.

. 5.

Phenotypic properties of phylogenetically diverse Ensifer strains. Ten Ensifer strains were screened for their ability to catabolize 190 carbon sources, and to grow in 96 osmolyte and 96 pH conditions using Biolog Phenotype MicroArray plates PM1, PM2, PM9, and PM10. Growth in each well was summarized on a scale of 0 (dark blue) through 9 (dark red), with higher numbers representing more robust growth. A larger version of this figure, in which each condition is labeled along the Y axis, is provided as supplementary figure S8, Supplementary Material online.

Discussion

We investigated the evolution of SNF within the genus Ensifer, which includes a mix of nitrogen-fixing bacteria and bacteria unable to fix nitrogen, as a model for the evolution of interkingdom mutualisms. Our results indicate that, despite SNF having likely evolved multiple times within the genus Ensifer, the symbiotic and nonsymbiotic strains are largely separated into two phylogenetic clades reminiscent of the general division of pathogenic and environmental strains between the genera Burkholderia and Paraburkholderia (Sawana et al. 2014). Although it is possible that this result will fail to remain true as more genome sequences are published, we believe the result to be robust as the current set of genome sequences are of strains isolated from diverse legumes and other diverse environments including the rhizosphere, pristine caves, and an abandoned mine. In addition to the prevalence of SNF, the two clades differ with respect to their genomic composition (pangenome content and genome size) and phenotypic properties (metabolic capacity and stress tolerance) based on initial tests of a subset of strains. There have been several revisions to the taxonomy of the genus Ensifer. Recently, a genome-based approach was proposed to standardize bacterial taxonomy (Parks et al. 2018) that splits the genus Ensifer into two genera: Sinorhizobium and Ensifer. The symbiotic and nonsymbiotic clades identified here correspond with the genera Sinorhizobium and Ensifer, respectively, supporting the proposal to divide the genus Ensifer into two genera.

Our analyses revealed a complex evolutionary history of SNF within the genus Ensifer. In addition to SNF emerging a predicted six times or more, we detected possible losses of SNF and allele switches. Between one and six of the NodABC and NifHDK proteins were not detected in five of the strains in the symbiotic clade (fig. 1). Although this may be indicative of multiple losses of symbiosis, we cannot rule out that these are false negatives due to incomplete genome assemblies or genome assembly errors; five of the six genomes were draft genomes, and the one strain with a complete genome (E. meliloti M162) can nodulate ten of 27 tested Medicago truncatula genotypes suggesting it does contain nod and nif genes (Sugawara et al. 2013). Based on the RepAB phylogeny (fig. 2), the symbiotic megaplasmid of Ensifer arboris likely shares common ancestry with the symbiotic megaplasmids of the sister species E. meliloti and E. medicae. Yet, the nodulation and nitrogen fixation genes appeared distinct. Thus, we hypothesize that there was a recent replacement of the symbiotic genes in E. arboris, or alternatively, in E. meliloti and E. medicae. This hypothesis is supported by the observation that the host ranges of these rhizobia differ; unlike E. meliloti and E. medicae, E. arboris does not nodulate plants of the genus Medicago (Zhang et al. 1991).

The reason for the phylogenetic bias in the evolution of SNF within the genus Ensifer remains unclear, especially considering that strains from both clades are plant-associated (Bai et al. 2015). One hypothesis is that each clade occupies distinct niches within the soil and plant-associated environments. Indeed, analysis of a subset of strains suggested that species of the nonsymbiotic clade have a broader metabolic capacity (supplementary table S6, Supplementary Material online), as well as a larger average genome size (fig. 3), which is consistent with these species being more capable of adapting to fluctuating nutritional environments. This is further supported by the apparently higher number of rRNA operons in strains of the nonsymbiotic clade, which is generally thought to allow bacteria to more quickly respond to changing nutrient conditions (Stevenson and Schmidt 2004; Roller et al. 2016). Moreover, the nonsymbiotic clade could be differentiated from the symbiotic clade based on its pangenome content (fig. 3), which leads us to hypothesize that strains of these clades acquire genes from distinct gene pools, further supporting the hypothesis that they belong to distinct gene-cohesive groups and ecological niches. This hypothesis may then be interpreted in the framework of the stable ecotype model (sensu Cohan 2006), where the symbiotic and nonsymbitic clades represent two, ecologically distinct and monophyletic groups and where periodic selection events (e.g., fitness for SNF) are recurrent.

An alternate, but not mutually exclusive, hypothesis is that the symbiotic Ensifer clade contains “facilitator” genes required to support SNF, similar to the theory that ancestral legumes contained a genetic “predisposition” necessary for the eventual evolution of rhizobium symbioses (Soltis et al. 1995; Doyle 2011; Werner et al. 2014). Conversely, evolution of SNF may also require the absence of “inhibitor” genes, such as the absence of virulence factors (Marchetti et al. 2010). As we did not evaluate cause-and-effect relationships, our data set does not definitely address these hypotheses. However, we observed numerous genotypic and likely phenotypic differences between the symbiotic and nonsymbiotic clade, providing some support for these hypotheses. For example, the tested strains of the symbiotic clade appeared to have higher tolerance to low pH (supplementary fig. S8 and data set S5, Supplementary Material online), which is notable as the curled root hair is an acidic environment (Hawkins et al. 2017). At the genomic level, 231 of the core genes of the symbiotic clade were absent from the pangenome of the nonsymbiotic clade and thus are good candidates as possible facilitators and follow-up studies. However, facilitators and inhibitors could also take the form of polymorphisms within highly conserved genes, as shown for bacA and the Sinorhizobium–Medicago symbiosis (diCenzo et al. 2017).

In summary, we show that the legume symbionts and nonsymbionts of the genus Ensifer are largely segregated into two phylogenetically distinct clades that differ in their genomic and phenotypic properties. We suggest that these observations, which follow the guidelines recently reported for rhizobia and agrobacteria (de Lajudie et al. 2019), support the division of the genus into two genera: Ensifer for the nonsymbiotic clade and Sinorhizobium for the symbiotic clade. However, formal descriptions and publication of the genera in the International Journal of Systematic and Evolutionary Microbiology (IJSEM) are still required. We also provide evidence that SNF genes were likely acquired several independent times within this genus, but predominately within one monophyletic clade. These observations suggest that the presence or absence of other genomic features (“facilitators” or “inhibitors”) aside from the core symbiotic genes could be required for the establishment of an effective symbiosis. This suggestion is supported by the ability to differentiate the strains of the two clades based on their pangenome content and, at least for the tested subset, their phenotypic properties. However, as cause-and-effect relationships were not examined, follow-up study is required to more directly test this facilitators hypothesis.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.