A Genomotaxonomy View of the Bradyrhizobium Genus

1Bradyrhizobium is one of the several genera of nitrogen fixing bacteria capable of forming symbiotic nodules in legumes. Bradyrhizobium strains were previously designated the slow growing Rhizobium and recognized as an independent genus in 1982. The number of species in Bradyrhizobium has increased largely in recent years1. Notably, there are photosynthetic bacteria among bradyrhizobia which need no Nod factors to induce nodules and this opened a new research area on plant nodulation. Non-symbiotic bradyrhizobia are dominant in forest soils.

2Different Bradyrhizobium species are the main nodule symbionts of important crop legumes such as soybean, Lima bean or peanuts and have been isolated from nodules of many tropical and temperate legumes. The ancestral symbionts in Phaseolus could have been bradyrhizobia with a later symbiont shift to Rhizobium in nodules of some temperate Phaseolus species.

3A large diversity of bradyrhizobia has been revealed with the sequence of few genes used as molecular markers. Within bradyrhizobia, 16S rRNA genes are known not to provide adequate sequence diversity to clearly recognize distinct species. For this, markers such as recA, glnII, atpD, dnaK, gyrB, and rpoB, have been frequently used to characterize bradyrhizobial species and strains. Novel markers such as ftsA gene provide congruent phylogenies to those derived from recA and glnII gene sequences. Bradyrhizobium diversity knowledge has expanded recently with studies from native legumes in Africa, from Brazilian and Indian Chamaecrista, from threatened native species in Brazil, from Genisteae plants in Poland, indigenous trees in China as examples. In addition, new bradyrhizobial species have been reported as well.

4Additionally, some bradyrhizobial strains are capable of fixing nitrogen as endophytes of some plants. Growth-independent approaches have shown a broad distribution of bradyrhizobia associated with roots of many non-legume plants such as rice, maize and pines and bradyrhizobia have been found in plant tumors and also in earthworm and insect guts.

5Nowadays there are publicly available genomes of 187 bradyrhizobial strains. Bradyrhizobia have characteristic large genomes with few plasmids and no symbiosis plasmids with one exception. In Bradyrhizobium chromosomes there are symbiosis islands which carry nod and nif genes that are responsible for nodulation and nitrogen fixation, respectively. Symbiosis markers that are commonly used toward symbiosis phylogenetic reconstructions are nifH and nod genes. In regard to nodulation genes, there is a large diversity of nodA or nodC genes and the phylogenies from these genes support a vertical and also a horizontal transfer of these genes among bradyrhizobia. Symbiovars in relation to host specificity have been identified for some bradyrhizobial species and novel symbiovars have been described.

6Novel metrics to recognize species are based on genome analysis and average nucleotide identity, ANI and phylogenomics are proving to be very useful toward this goal. A phylogenomic study of Bradyrhizobium strains showed that the presence and type of flagellum are phylogenetically determined. Other phylogenomic-based studies have been reported for the beta-rhizobia highlighting the existence of novel genera. It is the aim of this work to use a similar genomic-based approach and present an up-dated global genomic-based analysis of the Bradyrhizobium genus to further support its taxonomic status. Previously, such a study was performed for Rhizobium, Sinorhizobium and Agrobacterium, shedding light on these genera taxonomical designations.

7All Bradyrhizobium genome sequences available in GenBank as of December 2018 were retrieved. An all-versus-all average nucleotide identity (ANI) matrix was constructed using OrthoANI. Non-metric multidimensional analysis on the ANI matrix was performed with PAST to identify outlier genome sequences. Species-level clusters were defined at a 95% ANI cutoff value and representative genomes from each cluster were selected for further analysis. For the phylogenomic reconstruction, Prodigal was used for de novo prediction of protein-coding genes for each selected genome. Amino acid sequences of 31 conserved phylogenetic markers were retrieved from each genome using AMPHORA2 and aligned with muscle. Each alignment was processed with trimAl to identify and remove poorly aligned regions. All alignments were concatenated and a maximum likelihood phylogeny was constructed with PhyML. Tree node support was evaluated with bootstrap analysis of 1000 replicates.

8The presence of plasmid replication systems in the sequenced genomes was evaluated by looking for homologs to the replication protein B (RepB) using BLASTP searches against a database of all the Bradyrhizobium proteins obtained in the previous section. The number and size of genomic islands in complete and almost complete genomes were determined with the IslandViewer 4 server.

9Bradyrhizobium nucleotide sequences from the recombinase A protein (recA) and glutamine synthetase type II (glnII) genes were retrieved from the GenBank database and aligned using MAFFT. The multiple sequence alignments were trimmed with Bioedit recovering 336 or 454 bp fragments common to most recA or glnII sequences, respectively, shorter sequences were discarded. An all-versus-all identity distance matrix was constructed from the alignments with the dist.seqs command of mothur. Sequences were classified into species-level OTUs with the cluster command of mothur. OTU abundances were used to calculate species richness estimators with SpadesR2 and EstimateS3, and to construct rarefaction curves with PAST.

10One hundred and eighty-seven genomes of strains named as Bradyrhizobium were available in GenBank at the time of writing this manuscript in December 2018 (Supplementary Table 1). A non-metric multidimensional scaling graph of ANI values shared by these genomes allowed the identification of a core set of 180 related genomes plus 7 outliers (Figure 1). Upon examination of selected phylogenetic markers (rrs, recA, and/or glnII), only the 180 related genomes corresponded to bona fide bradyrhizobia. The outliers were misnamed strains belonging to other genera (see Supplementary Table 1 for details).

11Based on a 95% cutoff ANI value, the 180 bona fide bradyrhizobial genomes in GenBank were grouped into 93 species-level clusters (referred here as OTUs), the majority of which (72%) were represented by only one sequenced strain (Table 1). Twenty nine type strains were present among the 180 sequenced bradyrhizobia. The most represented species were Bradyrhizobium japonicum and Bradyrhizobium canariense, each with 15 strains; followed by B. diazoefficiens and B elkanii with 9 and 8 strains, respectively. The abundance of sequences from species able to nodulate soybean indicates the bias toward genomic studies directed at symbionts from this agronomical important legume.

12Besides Bradyrhizobium, other genera of the Bradyrhizobiaceae family with sequenced genomes are Bosea, Afipia, Rhodopseudomonas, Nitrobacter, Tardiphaga, Oligotropha, and Variibacter with 30, 24, 18, 6, 4, 3, and 1 sequences, respectively. A phylogenomic analysis of the family showed that most genera segregated as expected with the sole exception of Oligotropha and Afipia which intermingle (see Figure 2 for a condensed tree and Supplementary Figure 1 for a full tree). Bosea was the most distantly related genus in the family and its clustering with the other genera was not significantly supported by bootstrap analysis.

13Two Bradyrhizobium superclades are recognized based on analysis of the 16S rDNA gene. Superclade I includes B. japonicum and related species while superclade II contains B. elkanii and allied species. The new phylogenomic analysis supports the existence of those clades but revealed a more complex structure within Bradyrhizobium with additional groups (see Figure 3 for a condensed tree and Supplementary Figure 2 for a full tree).

14Superclade III encompassed several photosynthetic strains isolated from Aeschynomene indica, as well as Bradyrhizobium oligotrophicum which also forms nodules in that legume. Bradyrhizobium denitrificans, which comprise photosynthetic bacteria able to nodulate A. indica, also belong to superclade III as shown by MLSA analysis. Members of superclade III were the first to be recognized as having the capacity to induce nodulation in the absence of the canonical nod genes.

15Superclade IV includes strains isolated from several cultivated and wild legumes such as Phaseolus lunatus and Lupinus maria-josephi. Bacteria from this superclade closely group with B. elkanii in 16S rRNA phylogenies but its distinctiveness was first recognized by dnaK sequence analysis and MLSA and is now confirmed by phylogenomics. These bradyrhizobia have the characteristics of extra slow growth forming punctate colonies on YEM medium and usually strong alkali production.

16Superclade V was formed by only two genomes one of them from strain Tv2a.2 from Tachigali versicolor, that was previously shown to occupy instable positions in the Bradyrhizobium phylogeny depending on the marker used. Superclade VI also grouped only two genomes including that of strain Ai1a.2 a representative of neotropical bradyrhizobia which possess a characteristic insertion in their 23S ribosomal gene sequence. Finally, superclade VII included a single strain, GAS401, which was isolated from a forest soil in the United States.

17The wide genotypic and phenotypic diversity within Bradyrhizobium may suggest that it includes several genera. As early as 1990, strains now classified in superclade IV, were proposed to constitute the separate genus Photorhizobium because of their photosynthetic abilities and induction of stem nodulation. The most up to date, although not officially recognized, identity threshold for genus circumscription based on the 16S rRNA gene is 96.4%. When applied to bradyrhizobia, superclade II strains are different enough to be considered a separate genus, however, the same threshold indicates that Bradyrhizobium, Rhodopseudomonas, Nitrobacter, and Afipia may constitute a single genus. Superclade II strains possess recombinant segments in their 16S ribosomal gene that explains their sequence divergence from other bradyrhizobia but this case of localized recombination with other bacteria do not justify their separation from the genus.

18Recently, used the percentage of conserved proteins (POCP) as a genome metric for genus circumscription. According to their proposal, two bacteria may belong to the same genus if they share 50% or more of their proteins. In the case of Bradyrhizobium, all compared strains share >50% of their protein complements supporting a single genus. Different genera of the Bradyrhizobiaceae family had POCPs values between 40 and 50% except for most comparisons between bradyrhizobia and Rhodopseudomonas that can share more than half of their proteins. Thus, in general the POCP metric supports Bradyrhizobium as a single genus but indicate a close relationship with Rhodopseudomonas. Conservation of gene content may indicate not only phylogenetic relationship but also phenotypic similarity. Both Bradyrhizobium and Rhodopseudomonas include nitrogen fixing and photosynthetic strains, and, recently, a study suggest that Rhodopseudomonas can be mutualistic symbionts of some fungi. It will be worth to further explore the common features of both genera based on their genomic sequences.

19Bradyrhizobia are considered bacteria with large genomes, however the range of genome sizes in the genus is still unknown. A graphic showing the size distribution of the 180 genome assemblies of Bradyrhizobium strains is shown in Figure 4. Since only 15 of the genome sequences correspond to strains with closed replicon(s), we determined if the remaining sequences represent partial or whole genomes. A completeness analysis with BUSCO revealed that the four smaller assemblies, all <2 Mbp, represented partially sequenced genomes. Information available from those genomes revealed that they are derived from single cells or metagenomes likely explaining their incompleteness.

20The smallest sequence representing a complete genome, was 6.1 Mbp in size and corresponded to strain GAS165 isolated from a forest soil. This size is similar to that found in Rhizobium and Sinorhizobium genomes. The largest complete assembly (11.7 Mbp) was also from an isolate obtained from soil (strain GAS478). Interestingly, both GAS 165 and GAS478 strains which belong to superclade IV, lacked symbiosis genes which may indicate that genome size in Bradyrhizobium is not related to its ability to engage in symbiosis with legumes. It is worth noting that the largest bradyrhizobial genome is close in size to that of Sorangium cellulosum (13 Mbp) the bacteria with the largest known genome.

21Most Bradyrhizobium strains had genome sizes between 7 and 10 Mbp with a mean size of 8.6 Mbp (Figure 4). When compared to other genera in its family, bradyrhizobia have the largest genomes (Figure 5) followed by Tardiphaga and Bosea, while the smallest genomes were those from Nitrobacter. The three genera with the biggest genomes can interact with plants. On the other hand the smallest genomes are found in metabolically limited bacteria like Nitrobacter, Oligotropha or Variibacter or in intracellular pathogens like Afipia. Thus, genome size in the Bradyrhizobiaceae seems to be related to lifestyle.

22Given the large genome size in bradyrhizobia we wondered which the proportion of plasmids and genomic islands (GIs) was. Bradyrhizobia are typically regarded as unireplicon bacteria, however, strains with plasmids have been reported including a single case of a symbiosis plasmid. The presence of extra-chromosomal replicons was evaluated in the genome-sequenced strains by searching for genes coding for homologs of the plasmid partitioning protein RepB. We found repB genes in 35 genomes (Supplementary Table 2) with almost half of the strains (n = 17) possessing a single homolog which suggested the presence of a single plasmid. The remaining strains may harbor up to 6 plasmids, although it should be noted that a single plasmid could harbor two repABC operons. The few completely-sequenced plasmids ranged in size from 136 to 285 kbp which is a size range similar to that found by pulse-field gel electrophoresis analysis of plasmids in a diverse population of Bradyrhizobium. Plasmids can represent from 2.7 to 6.6% of the genome in a single strain. The sequence of the single Bradyrhizobium symbiosis plasmid known to date was reported as a scaffold that is larger than other plasmids in the genus (736 kbp) and which represents 9.4% of the corresponding genome.

23Mobile elements in unireplicon bacteria are typically present as clusters of genes known as genomic islands (GIs). As mutualistic symbionts of legumes, bradyrhizobia typically possess nodulation and nitrogen fixation genes grouped in a symbiosis island (SI). The percentage of a Bradyrhizobium chromosome devoted to GIs was calculated using the Island Viewer 4 server applied on completely sequenced genomes (Table 2). It was found that from 5.2% and up to 17.8% of the chromosome can be regarded as the GI mobilome of individual strains.

24Another aspect related to the wide diversity found among bradyrhizobia is related to their species richness. Up to December 2018, forty one Bradyrhizobium species have been described. As previously mentioned, there would be 93 species among the genome-sequenced Bradyrhizobium strains, i.e., more than twice the number of presently described species. A quick survey of some of the studies published on Bradyrhizobium diversity in different regions like North America, South America, Africa, Asia, and Australia suggest that the number of bradyrhizobial species must be larger than a hundred but the total number is yet unknown.

25To estimate the number of potential Bradyrhizobium species we decided to treat the 5678 recA and 3575 glnII sequences available in the GenBank database as meta-samples of the worldwide population of bradyrhizobia. Although other phylogenetic markers have also been used to characterize bradyrhizobia, numbers of their available sequences were much lower in comparison to recA and glnII (2411, 1771, 1758, 1465 for dnaK, rpoB, atpD, gyrB, respectively). We refrain to use the 16S rDNA gene despite having a large number of sequences in the databases because it is already known that sequences of this gene are too conserved in bradyrhizobia to discriminate between species and in some cases are even unable to discriminate between Bradyrhizobium and closely related genera.

26Pair wise comparisons between sequences from Bradyrhizobium type strains revealed that nucleotide identities of 98.2% for recA and 98.8% for glnII can be used as cutoff values to discriminate between currently described bradyrhizobial species. Using those cutoff levels, the worldwide sampled bradyrhizobia with sequenced recA genes can be clustered into 648 species-level OTUs while glnII data revealed 431 potential species. Rarefaction curves of both meta-samples did not leveled off (Figure 6) suggesting that there are still species to be discovered. Non-parametric richness estimation applied to both meta-samples indicated that there could be from 750 to 880 species of Bradyrhizobium in nature (Figure 7). It is worth noting that both meta-samples are composed primarily of root nodule isolates. Non-symbiotic bradyrhizobia inhabiting niches such as soil or rhizosphere, or as endophytes of non-legume plants may encompass additional species.

27The large genomic and phenomic diversity, as well as the high number of predicted species may indicate that the Bradyrhizobium genus is old. Interestingly, up till now this is the only nodule bacteria genus that contains photosynthetic bacteria and rhizobia that do not need Nod factors for nodulation. It is notable that it contains some very efficient strains for nitrogen fixation but also many generalists that form nodules in tropical legumes which in general are considered to precede temperate legumes. Several authors have noted that more “primitive” legume plants form nodules predominantly with Bradyrhizobium, thus lending support to the hypothesis that these rhizobia are their ancestral symbionts. Similarly, in Phaseolus the ancestral symbionts were bradyrhizobia with a later symbiont shift to Rhizobium in nodules of some temperate species. The large diversity and number of different nod genes in bradyrhizobia may be suggestive that nodulation arose in bradyrhizobia. This has been a subject of discussion with alternative views placing the origin of nod genes in beta-proteobacteria based on Nod factor exporter gene phylogenies. Since these exporters are not strictly required for nodulation, we can suppose that accessory nod genes such as those for transporting Nod factors may be later additions after the emergence of common nodulation genes in Bradyrhizobium, then we may conciliate other possible origins for accessory nod genes.

28During this article reviewing process, two papers that we must mention were published. The first one by stated that Bradyrhizobiaceae is an illegitimate name that needs to be replaced by Nitrobacteraceae because the latter contains Nitrobacter and was proposed earlier than Bradyrhizobiaceae and therefore takes precedence. Since this nomenclatural change did not affect our conclusions and because Bradyrhizobiaceae is how the family is still known, we choose to retain the name throughout this paper. However, we recognize that future studies will probably have to use Nitrobacteraceae as the proper family name. The second paper reported a phylogenomic analysis of bradyrhizobia and related taxa. Similar to our findings, confirmed the distinctiveness of the japonicum, elkanii and photosynthetic superclades and described additional infrageneric groups, albeit with some differences probably due to the use of different sets of genes for the phylogenomic reconstructions. In their analysis, our superclade IV was scattered in three groups named jicamae, soil 1 and soil 2; and our superclade VI was included in the elkanii group. In the first case we recognized that the larger number of genes used by in comparison to us, may have better resolved strains from our superclade IV, however, in the latter case we consider that superclade VI can be rightly segregated from the elkanii group by the presence of an insertion in their 23S ribosomal genes. Unlike us, tried to relate symbiotic and photosynthetic lifestyles with the infrageneric structure found within Bradyrhizobium but no clear relationship was found probably because those lifestyles are coded by accessory genes that can be gained and loss. On the other hand, our study includes aspects that were not covered by like species richness in Bradyrhizobium and the mobilome and putative plasmid content. We also provide a more thoroughly consideration of genome size in the Bradyrhizobiaceae and whether or not bradyrhizobia may be split into different genera. Thus, our analysis and those of provide complementary views on the highly diverse Bradyrhizobium genus.

29EO-O and EM-R contributed conception and design of the study. EO-O performed the bioinformatics analysis. Both authors analyzed the results and wrote manuscript, and contributed to manuscript revision, read and approved the submitted version.

30The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.