Comparative genomic and phylogenomic analyses of the Bifidobacteriaceae family

Background Members of the Bifidobacteriaceae family represent both dominant microbial groups that colonize the gut of various animals, especially during the suckling stage of their life, while they also occur as pathogenic bacteria of the urogenital tract. The pan-genome of the genus Bifidobacterium has been explored in detail in recent years, though genomics of the Bifidobacteriaceae family has not yet received much attention. Here, a comparative genomic analyses of 67 Bifidobacteriaceae (sub) species including all currently recognized genera of this family, i.e., Aeriscardovia, Alloscardovia, Bifidobacterium, Bombiscardovia, Gardnerella, Neoscardovia, Parascardovia, Pseudoscardovia and Scardovia, was performed. Furthermore, in order to include a representative of each of the 67 (currently recognized) (sub) species belonging to the Bifidobacteriaceae family, we sequenced the genomes of an additional 11 species from this family, accomplishing the most extensive comparative genomic analysis performed within this family so far. Results Phylogenomics-based analyses revealed the deduced evolutionary pathway followed by each member of the Bifidobacteriaceae family, highlighting Aeriscardovia aeriphila LMG 21773 as the deepest branch in the evolutionary tree of this family. Furthermore, functional analyses based on genome content unveil connections between a given member of the family, its carbohydrate utilization abilities and its corresponding host. In this context, bifidobacterial (sub) species isolated from humans and monkeys possess the highest relative number of acquired glycosyl hydrolase-encoding genes, probably in order to enhance their metabolic ability to utilize different carbon sources consumed by the host. Conclusions Within the Bifidobacteriaceae family, genomics of the genus Bifidobacterium has been extensively investigated. In contrast, very little is known about the genomics of members of the other eight genera of this family. In this study, we decoded the genome sequences of each member of the Bifidobacteriaceae family. Thanks to subsequent comparative genomic and phylogenetic analyses, the deduced pan-genome of this family, as well as the predicted evolutionary development of each taxon belonging to this family was assessed. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3955-4) contains supplementary material, which is available to authorized users.


Background
The Bifidobacteriaceae is the sole family member of the Bifidobacteriales order, and has been shown to represent the deepest branch within the Actinobacteria phylum [1]. Currently, the Bifidobacteriaceae family includes 55 (sub) species of the genus Bifidobacterium [1,2] and members of eight additional genera, i.e., Aeriscardovia, Alloscardovia, Bombiscardovia, Gardnerella, Neoscardovia, Parascardovia, Pseudoscardovia and Scardovia, which together encompass 12 species [1,3]. Furthermore, a novel, yet unculturable species was identified from termites and included in the Bifidobacteriaceae family with the taxonomic denomination of 'Candidatus Ancillula trichonymphae' [4]. The name of this latter organism originates from a variety of flagellates of the genus Trichonympha, of which this strain is symbiont [4].
Bifidobacteriaceae are chemoorganotrophs with a fermentative type of metabolism, Gram-positive, nonspore-forming, non-motile, and anaerobic or facultative anaerobic bacteria [5]. They reside in different ecological niches, such as the human and animal gastrointestinal tract (GIT), oral cavity and the (social) insect gut [6], while they may also be found in blood and sewage, possibly due to environmental contamination. Many bifidobacteria are appreciated for their purported healthpromoting activities as well as their relevance in early life colonization and contributions to the infant gut glycobiome [7]. Conversely, members of the other genera of the Bifidobacteriaceae family are generally associated with human and animal dental caries, and are commonly isolated from human clinical samples of tonsil abscesses and bacterial vaginosis [8,9]. Furthermore, in contrast to bifidobacteria, which mainly include strict anaerobes with some exceptions, such as Bifidobacterium animalis subsp. lactis and Bifidobacterium asteroides [10,11], other members of the family can grow under aerobic conditions and possess DNA with a lower G + C content [12][13][14].
The most controversial species belonging to this family is Gardnerella vaginalis [9], originally described by Leopold in 1953 and named Haemophilus vaginalis [15]. Subsequently, taxonomic studies mixed with data obtained from biochemical analyses and electron microscopic examinations, supported the need for its re-classification as a new genus [16,17]. Currently, G. vaginalis is described as an opportunistic pathogen whose presence is tightly associated with bacterial vaginosis [18,19]. Furthermore, Parascardovia denticolens and Scardovia inopinata, which were classified in 1996 as Bifidobacterium denticolens and Bifidobacterium inopinatum [20], respectively, together with Scardovia wiggsiae and Bifidobacterium dentium are associated with human dental caries [21]. While members of these species are present at high numbers in the saliva of adults, their presence strongly correlates with other caries-associated organisms [22,23].
Notably, comparative genome analyses of the genus Bifidobacterium have been targeting the entire genus [2,24] or one specific bifidobacterial taxa, i.e., B. bifidum [25], B. adolescentis [26], B. breve [27], B. longum [28] or the B. animalis subsp. lactis taxon [11]. In contrast, the genomics of the other eight genera belonging to the Bifidobacteriaceae family have not yet been investigated in detail. Here, we decoded the genomes of 11 species belonging to the Bifidobacteriaceae family for which there was no prior genomic data. Furthermore, we performed an in depth comparative genomic analysis, as well as a phylogenetic reconstruction of the 67 (sub) species currently assigned to the Bifidobacteriaceae family.

Bifidobacteriaceae strains
We retrieved the complete and partial genome sequences of 56 Bifidobacteriaceae strains from the National Center for Biotechnology Information (NCBI) public database (Table 1). Additionally, we sequenced and analyzed the genome sequences of 11 Bifidobacteriaceae strains deposited in the GenBank sequence database (Table 2).

Bacterial strains and growth condition
Bifidobacteriaceae pure cultures were inoculated in de Man-Rogosa-Sharpe (MRS) medium (Scharlau Chemie) supplemented with 0.05% (wt/vol) L-cysteine hydrochloride and were grown in an anaerobic atmosphere (2.99% H 2 , 17.01% CO 2 , and 80% N 2 ) in a chamber (Concept 400, Ruskin) at 37°C for 16 h. DNA was extracted as described previously [29] and subjected to further phenol-chloroform purification using a previously described protocol [30].

Genome sequencing and assemblies
DNA extracted from the various Bifidobacteriaceae strains was subjected to whole genome sequencing using MiSeq (Illumina, UK) at GenProbio srl (Parma, Italy) following the supplier's protocol (Illumina, UK). Fastq files of the paired-end reads obtained from targeted genome sequencing of the isolated strains were used as input for genome assemblies through the MEGAnnotator pipeline [31]. The MIRA program (version 4.0.2) was used for de novo assembly of each Bifidobacteriaceae genome sequence [32].

Sequence annotation
Protein-encoding open reading frames (ORFs) were predicted using Prodigal [33]. Transfer RNA genes were identified using tRNAscan-SE v1.4 [34], while ribosomal RNA genes were detected using RNAmmer v1.2 [35]. Results of the gene-finder program were combined with data from RAPSearch2 analysis (Reduced Alphabet based Protein similarity Search) [36] of a non-redundant protein database     [37]. The combined results were inspected by Artemis [38], which was used for manual editing purposes aimed at verifying and, where necessary, redefining the start of each predicted coding region, and to remove or add coding regions.

Bifidobacteriaceae pan-genome analysis
For the 67 genome sequences of each member of the Bifidobacteriaceae family, a pan-genome calculation was performed using the PGAP pipeline [39]. The ORF contents from all genomes used in this study were organized in functional clusters using the GF (Gene Family) method involving comparison of each protein to all other proteins using BLAST analysis (cutoff E-value of 1 X 10 −5 and 50% identity over at least 50% of both protein sequences), followed by clustering into protein families, named Bifidobacteriaceae-specific clusters of orthologous groups (BaeCOGs), using MCL (graph-theory-based Markov clustering algorithm) [40]. A pan-genome profile was built using an optimized algorithm incorporated in PGAP software, based on a presence/absence matrix that included all identified BaeCOGs in the analyzed genomes. Following this, unique protein families for each of the 67 Bifidobacteriaceae genomes were classified. Protein families shared between all genomes, named core BaeCOGs, were defined by selecting the families that contained at least one protein member for each genome.

Functional analysis
The prediction of genes encoding enzymes that possess structurally-related catalytic and carbohydrate-binding modules catalyzing hydrolysis, modification or synthesis of glycoside bounds was performed by means of the CAZy database [44]. Functional annotation of each gene was performed employing the eggNOG database [45]. A survey of complete pathways involved in both primary and secondary metabolism was performed by means of the MetaCyc metabolic pathways database [46]. Gene function was predicted using a cutoff E-value of 1 × 10 −10 to identify the best hit from each database.

Gene gain/loss through evolution reconstruction
Predicting gene acquisition or gene loss as a result of evolution of the bacterial species with at least four available genomes was performed with Count software [47] using Dollo's parsimony.

Results and discussion
General genome features of Bifidobacteriaceae genomes  (Table 1). Nevertheless, the Bifidobacteriaceae strains that do not belong to the Bifidobacterium genus possess at least one tRNA gene for each of the 20 amino acids (Additional file 1: Table S1). Additionally, a deeper screening of the anticodon sequences for each (sub) species does not display major differences, except for a lower abundance of the anticodon GGG in Proline tRNA of the Bifidobacteriaceae strains that do not belong to the Bifidobacterium genus (Additional file 1: Table S1). Consequently, it may be argued that the lower number of tRNA genes among certain members of the Bifidobacteriaceae family is not associated with a simplification of the codon usage of these strains. Furthermore, while the bifidobacterial genomes contain between one and six rRNA loci, with an average of 3.2 per genome [2], the genomes of the Aeriscardovia, Alloscardovia, Bombiscardovia, Gardnerella, Neoscardovia, Parascardovia, Pseudoscardovia and Scardovia genera exhibit a lower average number of rRNA loci, i.e., 2.6, which is consistent with the less extensive ORFome and tRNA arsenal identified in the corresponding genomes (Table 2). Furthermore, in silico analyses of the 12 genomes of the Bifidobacteriaceae family based on the Virulence Factor Database (VFDB) [53], did not reveal the occurrence of any virulence genetic determinants. Such results confirm previously reported findings for the genome of B. dentium Bd1 [23].
Pan-genome, core genome and unique genes of the Bifidobacteriaceae family Comparative genome analyses involving the 67 (sub) species belonging to the Bifidobacteriaceae family were performed to unveil the corresponding pan-genome, core genome and unique genes of this bacterial family. All genomes were subjected to identical ORFs finding and annotation protocols [31] in order to generate comparable data sets for each Bifidobacteriaceae taxa. A total of 25,744 BaeCOGs (Bifidobacteriaceae-specific clusters of orthologous genes) were identified in the 67 Bifidobacteriaceae (sub) species, of which 8359 had members present in at least two genomes. The pangenome size, when plotted versus the number of included genomes, clearly shows that the power trend line has yet to reach a plateau (Fig. 1). Actually, the number of new genes discovered by sequential addition of genome sequences was reduced from 839 to 636 BaeCOGs in the first three genomes additions to a number that ranged from 274 to 272 BaeCOGs in the final three additions, demonstrating the existence of an open pangenome within Bifidobacteriaceae family. This finding suggests that more members of the Bifidobacteriaceae family have yet to be identified, especially members of the family that do not belong to the Bifidobacterium genus, as these remain poorly characterized in various environments compared to the currently recognized bifidobacterial (sub)species.
Pan-genome analysis of the Bifidobacteriaceae family allowed the identification of 353 COGs shared by all 67 (sub) species, representing the core genome of currently sequenced Bifidobacteriaceae representatives (core Bae-COGs). An examination of the functional annotation of Core BaeCOGs employing the eggNOG database [45] shows that the most conserved core genes specify housekeeping functions such as replication, transcription and translation, or functions related to adaptation such as carbohydrate, nucleotide and amino acid metabolism as well as cell envelope biogenesis (Fig. 1).
The pan-genome analysis also allowed the identification of truly unique genes (TUGs) of the Bifidobacteriaceae family, i.e., those genes that are presented in one particular strain yet absent in any of the other examined representative of the Bifidobacteriaceae family. The number of TUGs range from 42 for B. indicum LMG 11587 to 585 for Bifidobacterium cuniculi LMG 10738. EggNOG analysis showed that the majority of TUGs (59%) have no functional annotation (Additional file 1: Table S2). Nevertheless, taking into account the classified genes through the eggNOG analysis excluding the hypothetical and no-function genes, the highest number of genes fall in carbohydrate metabolism and cell envelope biogenesis together with replication and transcription. As mentioned above, these are the same four categories identified as containing the highest numbers of core BaeCOGs. Interestingly, the functional annotation of TUGs revealed that bifidobacterial genomes exhibit a higher abundance of TUGs involved in carbohydrate metabolism and cell envelope biogenesis compared to those of other members of the Bifidobacteriaceae family, reflecting a 38% and 78% of additional TUGs based on the average numbers between groups, respectively [2,54] ( Fig. 1 and Additional file 1: Table S2). These results are consistent with previous genomic and functional analyses based on the reference strains for each (sub) species of the Bifidobacterium genus indicating that bifidobacteria are under strong selective pressure to acquire and retain accessory genes for carbohydrate utilization in order to be competitive in the specific ecological niches in which they reside [2,54].

Phylogenomic analyses of members of the Bifidobacteriaceae family
The availability of genome sequences for each member of the Bifidobacteriaceae family allowed an in-depth analysis of the evolutionary development followed by each member of this extensive family. A phylogenetic supertree was constructed based on the concatenation of 314 protein sequences that represent the core BaeCOGs with the exclusion of paralogs identified in each genome (Fig. 2).  The generated phylogenetic supertree showed that ten strains of the Bifidobacteriaceae species that do not belong to the Bifidobacterium genus represent the deepest branches of this supertree, being separated from the 55 bifidobacterial (sub) species (Fig. 2). Consequently, the bifidobacterial species were positioned at the top of the supertree, reflecting their highest average gene ratio (1.24) as compared to the non-bifidobacterial members of the Bifidobacteriaceae family. Thus, these data clearly indicate that evolution of currently known bifidobacterial species involved a relatively limited number of ancestral gene loss incidences, yet an extensive number of gene acquisition events, corroborating previously published data [2]. Interestingly, two members of the Bifidobacteriaceae family that do not belong to the Bifidobacterium genus, i.e., G. vaginalis ATCC 14018 and Bo. coagulans DSM 22924, appear to possess a higher level of phylogenetic relatedness with bifidobacterial strains as compared to other, non-bifidobacterial members of the Bifidobacteriaceae family (Fig. 2). While G. vaginalis ATCC 14018 shares the same phylogenetic branch as that of Bifidobacterium subtile LMG 11597, Bo. coagulans DSM 22924 is positioned within the deepest branch of the Bifidobacterium genus, i.e., within the B. asteroides group [2,10], exhibiting higher genome relatedness with Bifidobacterium actinocoloniiforme DSM 22766, which was also isolated from Bombus. In order to validate the branching of these two Bifidobacteriaceae strains, a phylogenetic tree based on the 16S rRNA genes was constructed, as well as a tree based on five housekeeping genes including hsp60, rpoB, dnaJ, dnaG and clpC (Additional file 2: Figure S1). While Bo. coagulans DSM 22924 shares the same phylogenetic branch with B. actinocoloniiforme DSM 22766 in both trees, G. vaginalis ATCC 14018 occupies different positions within these phylogenetic trees. Nonetheless, the housekeeping-based tree confirmed the position of G. vaginalis ATCC 14018 within the Bifidobacterium genus, while in the 16S rRNA-based tree it is placed between bifidobacterial species and the other non-bifidobacterial species. These findings cast doubts on the correct taxonomical classification of Bo. coagulans DSM 22924, and reinforce the importance of a phylogenomic approach as a tool for taxonomic validation [55].
Furthermore, we investigated the occurrence of genes predicted to encode enzymes for anaerobic respiration in the pangenome of the members of the Bifidobacteriaceae as previously described for the genome of B. asteroides PRL2011 [10]. Such in silico analyses highlight the presence of a cytochrome bd oxidase-encoding complex in the genome of Bo. coagulans DSM 22924 (Additional file 1: Table S3), including the cydA and cydB, which code for structural subunit of the cytochrome, as well as cydC and cydD, encoding a transporter required for the cytochrome assembly [10]. Furthermore, these four genes were identified in 11 genomes of the Bifidobacteriaceae family, including six strains isolated from insects, i.e., B. actinocoloniiforme DSM 22766, Bifidobacterium bohemicum DSM 22767, Bifidobacterium bombi DSM 19703, Bifidobacterium commune R-52791, Bifidobacterium coryneforme LMG 18911 and B. indicum LMG 11587, highlighting a correlation between the presence of the cytochrome bd oxidase complex and the ecological niche of isolation (Additional file 1: Table S3).
Based on the Bifidobacteriaceae supertree reconstruction, five out of the six bifidobacterial species that were newly sequenced as part of this study are positioned within one of the previously identified bifidobacterial groups [24,56]. Interestingly, four strains appear to belong to the Bifidobacterium longum group, i.e., B. myosotis DSM 100196, B. reuteri DSM 23975, B. eulemuris DSM 100216 and B. lemurum DSM 28807 (Fig. 2). This finding is in line with the particular ecological origin of these strains. In fact, each one of these four species was isolated from feces of monkeys, similar to the other two members that had previously been assigned to the B. longum group, i.e., Bifidobacterium stellenboschense DSM 23968 and Bifidobacterium callitrichos DSM 23973. These findings therefore highlight that the B. longum group includes bifidobacterial (sub) species isolated from humans and other related primates (Fig. 2). Furthermore, B. hapali DSM 100202 exhibits a genetic relatedness with Bifidobacterium biavatii DSM 23969, which belongs to the B. bifidum group [24]. Notably, B. tissieri DSM 100201 was shown to occupy a unique position within the Bifidobacteriaceae supertree. This observation was confirmed through an average nucleotide identity (ANI) analysis showing the highest value of 88.02% with the most related Alloscardovia criceti DSM 17774 strain. Notably, an ANI value below to 95% is assumed to be sufficient to classify that taxon as a distinct species [43]. B. tissieri DSM 100201 is therefore an interesting strain for further investigation due to its genetic divergence from other bifidobacterial taxa. In order to validate the species assignment for the other 10 sequenced members of the Bifidobacteriaceae family (Table 2) from a genomic prospective, the decoded genomes were subjected to ANI comparisons with the other 56 (sub) species that had been sequenced. The sequenced genomes of the 10 Bifidobacteriaceae strains showed ANI values below 95%, confirming their status as distinct species, with a relatively high value of 93.8% observed between B. eulemuris DSM 100216 and B. lemurum DSM 28807 (Additional file 1: Table S4).
Furthermore, the Bifidobacteriaceae supertree confirmed the phylogenetic relatedness of the recently decoded genome sequences of B. commune R-52791 and Bifidobacterium aesculapii DSM 26737, which occupy positions within the same branches with B. bohemicum DSM 22767 and B. stellenboschense DSM 23968, respectively [49,57]. Such data further confirms the existence of a direct relatedness between the ecological origin of both strains, i.e., insects of the Bombus genus and monkeys of the Callitrichidae family.
Enzymatic profiling and evolutionary development of the Bifidobacteriaceae family A functional profiling analysis was performed to assess the presence of genes encompassing carbohydrate, amino acid and fatty acid degradation pathways in each (sub) species of the Bifidobacteriaceae family. Normalizing the obtained number of gene matches with the overall genetic arsenal (i.e. total amount of genes) of each strain, non-bifidobacterial strains retrieved percentages that were slightly lower compared to the bifidobacterial strains, i.e., 7%, 11% and 9%, respectively (Additional file 1: Table S5). Interestingly, all members of the Bifidobacteriaceae family possess the Bifidobacterium shunt pathway [6,58], including the gene xfp that encodes the enzyme D-xylulose 5-phosphate phosphoketolase/D-fructose 6-phosphate phosphoketolase, expanding the notion of this signature metabolic trait of bifidobacteria to the whole family.
To further investigate the carbohydrate utilization abilities encoded by the genomes of the 67 Bifidobacteriaceae (sub) species, an enzyme classification toward glycans was performed. This enzyme classification was based on the Carbohydrate Active Enzymes (CAZy) database [44], which encompasses all currently known genetic determinants involved in the breakdown and utilization of carbohydrates, and revealed that the pangenome of the Bifidobacteriaceae family includes 9742 genes predicted to encode carbohydrate-active enzymes, i.e., glycosyl hydrolases (GHs), glycosyl transferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs) and carbohydrate binding modules (CBMs), present at 43.4%, 43.8%, 0.2% 5.9% and 12.6%, respectively. This very substantial number of retrieved enzymes reflects the findings of a previous analysis conducted on the 47 type-strains of the Bifidobacterium genus, where the glycan-breakdown potential of bifidobacteria in the mammalian gut was subjected to an extensive scrutiny [54].
Focusing on GH identification, 3989 genes were predicted on the analyzed genomes of 55 bifidobacterial (sub) species, while the remaining 12 Bifidobacteriaceae genomes contain 573 such genes. The Bifidobacteriaceae genomes specify a large arsenal of GH families, were GH13, GH3, GH43, GH23, GH32 and GH25 outnumber the other identified families (Additional file 1: Table S6). Despite the higher average amount of predicted GHencoding genes for a given bifidobacterial genome, i.e., 72, as compared to that for a non-bifidobacterial member of the Bifidobacteriaceae family, i.e., 48, normalization of GH counts against the total amount of predicted genes provided similar GH indexes, i.e., 0.039 and 0.032, respectively (Additional file 1: Table S6). Nonetheless, genomes of the 20 strains that exhibit the highest GH index all belong to the Bifidobacterium genus, and are nearly all isolated from fecal samples of humans and monkeys, chief among them being Bifidobacterium catenulatum LMG 11043, B. eulemuris DSM 100216, B. hapali DSM 100202 and B. biavatii DSM 23969 (Table 1). These results suggest that bifidobacteria that reside in the primate/human gut have enjoyed a relatively high number of adaptive events related to carbohydrate metabolism to benefit from a wider source of different nutrients present in this particular environment. In contrast, Bifidobacteriaceae species that exhibit a relatively low GH index originate from a broad spectrum of environments, such as the gut of insects or other animals, and sewage (Table 1). Furthermore, genomes with a low GH index correspond with the smaller genomic complement of certain members of the Bifidobacteriaceae family, probably due to gene decay characteristic of those microbes that are considered harmful for human health, such as G. vaginalis, A. omnicolens, S. inopinata and S. wiggsiae [59].
While gene acquisition events that occur during evolution of microbial genomes support adaptation to new ecological niches, gene loss on the other hand contributes to genome simplification in order to preserve energy and biological compounds [59,60]. Prediction of the complete Bifidobacteriaceae glycobiome content allowed us to estimate the acquisition and loss rates of genes encoding carbohydrate-active enzymes within this saccharolytic family. In order to depict gain and loss events of genes with a predicted function in carbohydrate metabolism, we collected the BaeCOGs that include GH-encoding genes obtained from the pan-genome analysis. The resulting 846 BaeCOGs were allocated among the Bifidobacteriaceae supertree showing that evolution of the current bifidobacterial (sub) species involved only a limited number of ancestral gene loss events, yet a substantial number of GHencoding gene acquisitions (Fig. 3). Thus, genes encoding GHs appear to have been acquired early in the evolution of bifidobacteria, followed by a simplification of the GHassociated gene arsenal that has resulted in or followed specialization toward those ecological niches in which current bifidobacterial species have been identified. Interestingly, members of the B. bifidum group possess the highest number of GH-encoding gene acquisitions compared to those of other groups, probably in order to expand its metabolic ability towards different carbon sources present in the host, similar to what was mentioned above regarding GH index discrepancies ( Fig. 3 and Table 1). Taking into account only those BaeCOGs that include members of GH families known to be involved in hostglycan degradation, i.e., GH20, GH29, GH33, GH38, GH95, GH101, GH112, GH125 and GH129 [54], the B. bifidum group, once again, displays the highest number of GH-encoding gene acquisitions (up to seven BaeCOGs). In this context B. bifidum exhibits the highest host-glycan degradation BaeCOGs acquisition number of the Bifidobacteriaceae family, thereby highlighting the capability of this bifidobacterial species to feed on host-glycan [25,61,62].

Conclusions
A lot of effort has been invested in the dissection and characterization of the genomic content from different members of the Bifidobacterium genus [2,11,[24][25][26][27][28]. In contrast, very little is known about the genomics of other members of the Bifidobacteriaceae family, which include apart from the Bifidobacterium genus eight additional genera. In this study, genome sequencing allowed us to explore the genome content of known members across the Bifidobacteriaceae family, as represented by 67 (sub) species, and to scrutinize the phylogenetic relatedness between each taxon belonging to this family. Fig. 3 Analysis of evolutionary gain and loss events as based on predicted GHs within the Bifidobacteriaceae family. The image displays the core gene-based supertree of the Bifidobacteriaceae family, where each node reports the number of predicted GHs identified for each strain. Furthermore, gain and loss events are indicated by green and orange bars on the edge leading to each node, while the numbers placed at the supertree leaves represent the related bifidobacterial taxa presented in Table 1 Bifidobacteria exhibit a higher number of genes per genome compared to other members of the Bifidobacteriaceae family, perhaps reflecting an increased competitiveness based on a broad spectrum of ecological niches from which they were isolated. The more complex genome content of bifidobacteria is also reflected by the capability of these microorganisms to degrade multiple carbon sources [54,63]. Such findings were further validated in this study by the analyses of the pangenome of the members of the Bifidobacteriaceae family, which highlights the occurrence of a higher abundance of TUGs in bifidobacterial genomes dedicated to the carbohydrate metabolism and cell envelope biogenesis, when compared to such numbers for other members of the Bifidobacteriaceae family. Enzymatic gene profiling revealed that the 20 strains that showed the highest GH indexes belong to bifidobacterial strains that are nearly all isolated from humans or monkeys. Moreover, a gene gain/loss analysis shows that members of the Bifidobacterium genus isolated from such primates possess the highest number of GH-encoding gene acquisitions, probably in order to expand their ability to metabolize different carbon sources. These results highlight a relatively large number of adaptive events related to carbohydrate metabolism among members of the Bifidobacterium genus that reside in omnivorous organisms that consume a wide variety of nutrients.
Furthermore, our phylogenomic analysis revealed possible taxonomic inconsistencies in the classification of G. vaginalis and Bo. coagulans, which displayed a close phylogenetic relatedness with other bifidobacterial strains, i.e., B. subtile LMG 11597 and B. actinocoloniiforme DSM 22766. Such findings further corroborate the strengths of genome-based analyses as an essential approach to be incorporated in phylogenome-based taxonomic studies.