Decoding the Genomic Variability among Members of the Bifidobacterium dentium Species

Members of the Bifidobacterium dentium species are usually identified in the oral cavity of humans and associated with the development of plaque and dental caries. Nevertheless, they have also been detected from fecal samples, highlighting a widespread distribution among mammals. To explore the genetic variability of this species, we isolated and sequenced the genomes of 18 different B. dentium strains collected from fecal samples of several primate species and an Ursus arctos. Thus, we investigated the genomic variability and metabolic abilities of the new B. dentium isolates together with 20 public genome sequences. Comparative genomic analyses provided insights into the vast metabolic repertoire of the species, highlighting 19 glycosyl hydrolases families shared between each analyzed strain. Phylogenetic analysis of the B. dentium taxon, involving 1140 conserved genes, revealed a very close phylogenetic relatedness among members of this species. Furthermore, low genomic variability between strains was also confirmed by an average nucleotide identity analysis showing values higher than 98.2%. Investigating the genetic features of each strain, few putative functional mobile elements were identified. Besides, a consistent occurrence of defense mechanisms such as CRISPR–Cas and restriction–modification systems may be responsible for the high genome synteny identified among members of this taxon.


Introduction
For millennia, animals and bacteria have coevolved and coadapted to establish interspecies relationships [1]. In this context, the gut microbiota composition of humans has been studied to identify which species grew together with their host and their role in human health [2,3]. Among hundreds of bacterial species that inhabit humans, members of the Bifidobacterium genus were identified as one of

Isolation of B. dentium Strains
One gram of each fecal sample was mixed with 9 mL of phosphate-buffered saline (PBS), pH 6.5. Serial dilutions and successive plating were performed using de Man-Rogosa-Sharpe (MRS) agar (Scharlau Chemie, Barcelona, Spain), supplemented with 0.05% (w/v) L-cysteine hydrochloride and 50 µg/mL mupirocin (Delchimica, Italy). The addition of the latter chemicals was used to enhance the selective growth of members of the Bifidobacterium genus. Agar plates were incubated in a chamber (Concept 400; Ruskin, Bristol, UK) with an anaerobic atmosphere (2.99% H 2 , 17.01% CO 2 , and 80% N 2 ) at 37 • C for 48 h. Among fecal samples reported in this study, about 1000 colonies were selected and subcultivated in MRS broth supplemented with 0.05% (w/v) L-cysteine hydrochloride in an anaerobic chamber at 37 • C for 16 h. Morphologically different colonies grown on MRS plates were randomly picked and restreaked to isolate purified bacterial strains. Then, DNA was extracted from all isolates and characterized as previously described [36]. Specifically, DNA was submitted to a B. dentium species-specific identification using a PCR methodology with primers Bdent1 5 -ATCCCGGGGGTTCGCCT-3 and Bdent2 5 -GAAGGGCTTGCTCCCGA-3 , which had been designed on the 16S rRNA gene sequence of this taxon. PCR amplification was performed according to the following protocol: one cycle of 94 • C for 5 min, followed by 30 cycles of 94 • C for 30 s, 54 • C for 30 s and 72 • C for 50 s, and finally one cycle of 72 • C for 5 min. This procedure allowed the identification of 32 strains belonging to the B. dentium species estimating a charge of B. dentium on the selective media of 3.2%.

Microbial DNA Extraction
B. dentium strains were inoculated in de Man-Rogosa-Sharpe (MRS) (Scharlau Chemie, Barcelona, Spain) medium to perform chromosomal DNA extraction. The latter medium was supplemented with 0.05% (w/v) L-cysteine hydrochloride and incubated at 37 • C in an anaerobic atmosphere (2.99% (v/v) H 2 , 17.01% (v/v) CO 2 , and 80% (v/v) N 2 ) employing an anaerobic chamber (Concept 400, Ruskin). Cells from 10 mL of an overnight culture were harvested by centrifugation at 6000 rpm for 8 min. The obtained cell pellet was used for DNA extraction using the GenElute TM Bacterial Genomic DNA kit (Sigma-Aldrich, Darmstadt, Germany) following the manufacturer's guidelines. Then, internal transcribed spacer (ITS) sequences were amplified from extracted DNA using primer pair Probio-bif_Uni/Probio-bif_Rev [37] and sequenced to avoid the decoding of clonal strains. Finally, sequenced ITS sequences were compared to a public database composed of bifidobacterial ITS sequences (http://probiogenomics.unipr.it/pbi/) involving the basic local alignment search tool (BLAST). This procedure allowed for discarding clonal strains within the same sample and to reduce the number of different B. dentium strains from 32 to 18. Table 1 lists the isolated strains employed in this study.

B. dentium Chromosomal Sequencing and Assemblies
The chromosomal DNA of 18 B. dentium strains was decoded by GenProbio Srl (http://genprobio. com) using a MiSeq platform (Illumina, San Diego, CA, USA) according to the supplier's protocol as previously described [30] by using the Nextera XT DNA Library Prep Kit (Illumina). The library samples obtained were then pooled into a Flow Cell V3 600 cycle (Illumina). Fastq files of paired-end reads generated from each genome sequencing effort were used as input for the genome assembly through the MEGAnnotator pipeline [38]. The SPAdes program v3.14.0 was used for de novo assembly of each bifidobacterial genome sequence with the pipeline option "-careful" and a list of k-mer sizes 21, 33, 55, 77, 99, 127 [39]. MEGAnnotator then employed contigs greater than 1000 bp to predict protein-encoding open reading frames (ORFs) using Prodigal [40]. Predicted ORFs were then functionally annotated using RAPSearch2 (reduced alphabet based protein similarity search) (cutoff e-value of 1 × 10 −5 and minimum alignment length 20) employing the NCBI reference sequences (RefSeq) database [41] together with hidden Markov model profile (HMM) searches (http://hmmer.org/) performed against the manually curated Pfam-A database (cutoff e-value of 1 × 10 −10 ). Then, tRNA genes were detected through tRNAscan-SE v1.4 [42], while rRNA genes were identified by means of RNAmmer v1.2 [43]. To guarantee a consistent genomic analysis, ORFs from each B. dentium genome sequence retrieved from NCBI public databases were repredicted and reannotated employing the identical bioinformatics pipeline used for the 18 B. dentium strains isolated in this study. The presence of extrachromosomal DNA was investigated using a curated database of plasmid sequences downloaded from NCBI [44] and using the tool mlplasmids v1.0 [45].

Comparative Genomic Analysis
The genome content of the 38 reconstructed B. dentium strains (Table 1) was submitted to a pangenome calculation by means of the pan-genomes analysis pipeline (PGAP) [46]. Predicted proteome of a specific B. dentium strain was screened for orthologues against the proteome of the other collected B. dentium strain employing BLAST analysis [47] (cutoff e-value of <1 × 10 −5 and 50% identity over at least 80% of both protein sequences). The data obtained were then clustered into protein families, also named as clusters of orthologous groups (COGs), employing MCL (graph theory-based Markov clustering algorithm) [48] by means of the method gene family (GF). Then, using the PGAP software, a pangenome profile was built based on a presence/absence matrix including all identified COGs of the 38 analyzed strains. The core genome of the B. dentium species was then defined based on protein families shared among each analyzed genome. Furthermore, truly unique genes (TUGs) encoded by a single analyzed genome that are not present in other B. dentium chromosomes were also identified.

Phylogenetic and Phylogenomic Analyses
The concatenated sequence of amino acids belonging to the core genome of each B. dentium was aligned by means of the MAFFT software [49], and the resulting phylogenetic trees was built using the neighbor-joining method through the ClustalW v2.1 program [50]. Then, the graphical viewer of phylogenetic trees FigTree v1.4 (http://tree.bio.ed.ac.uk/software/figtree/) was used to build a B. dentium phylogenetic tree. A value of average nucleotide identity (ANI) was calculated using the program fastANI using each genome pair as input [51].

Genomic Analyses
The proteome of each B. dentium strain was screened for genes predicted to encode carbohydrate-active enzymes based on sequence similarity to genes classified in the carbohydrate-active enzyme (CAZy) database [52]. Each predicted proteome of a given bifidobacterial strain was screened for orthologues against data of 17,538 bacterial genomes using HMMER v3.3 [53] (cutoff e-value of 1 × 10 −15 ) and BLASTp analysis [47] (cutoff e-value of 1 × 10 −30 ). Thus, a screening was performed employing the dbCAN2 meta server [54], followed by a BLASTp validation of the obtained results. The values of the glycosyl hydrolases (GH) index were attributed to dividing the number of GH counts of each bifidobacterial type strain against the total number of its predicted genes ( Table 1). The presence of genes encoding restriction enzymes was performed using the REBASE database [55]. The prediction of clustered regularly interspaced short palindromic repeats (CRISPR) and associated Cas-encoding genes was achieved through CRISPRCasFinder v4.2.2 [56]. The presence of genes that are predicted to be acquired by horizontal gene transfer (HGT) events was performed using COLOMBO v4.0 [57], with a sensitivity value of 0.4 as previously applied to type strains of the genus Bifidobacterium [58]. Identification of prophagelike elements was assessed using a custom database composed of bifidobacterial genes previously classified as prophage genes from each type strain [59]. Then, predicted prophages were further profiled by means of the PHAge Search Tool Enhanced Release (PHASTER) [60]. Each predicted proteome was also screened for the presence of insertion-sequence (IS) elements through the IS Finder online tool (https://isfinder.biotoul.fr/). Predicted genes were manually evaluated to remove false positives from each genomic analysis. B. dentium proteome was screened for proteins that can act as antibiotic resistance proteins by employing the MEGARes database through a BLASTp analysis (cutoff e-value of 1 × 10 −30 ) [61]. Putative bacteriocin-encoding genes were predicted using BAGEL3 with default parameters [62]. Toxin-antitoxin systems were predicted using the TADB v2.0 database through a BLASTp analysis (cutoff e-value of 1 × 10 −5 ) [63].

Carbohydrate Growth Assays
B. dentium strains were cultivated on semisynthetic MRS medium supplemented with a 0.5% (w/v) concentration of a specific sugar, and the optical densities (measured at a wavelength of 600 nm) were recorded using a plate reader (BioTek, Winooski, VT, USA). The plate was read in intermittent mode, with absorbance readings performed at 3 min intervals for three times after 48 h of growth, where each reading was ahead of 30 s of shaking at medium speed. Cultures were grown in biologically independent triplicates, and the reported growth data were expressed as the means of these replicates. Carbohydrates used in this study were acquired from Sigma and Carbosynth (Berkshire, UK).

Genomic Overview of the B. dentium Taxon
Many published studies discussed the peculiar ecological niche of B. dentium, highlighting the occurrence of this taxon with dental caries [20][21][22][23][24]. However, a genomic comparison between the sequenced B. dentium strains is somewhat limited, and comparative genome analyses characterizing the genomic diversity and genetic features of members of this species have not been described so far. Thus, in the framework of a previous ecologic survey of bifidobacterial communities in mammals, using a culture-dependent approach, we isolated 17 different B. dentium strains from fecal samples of primates and one single strain from a fecal sample of Ursus arctos (Table 1) [26]. More specifically, we isolated B. dentium strains from seven primate species, i.e., Colobus guereza, Eulemur macaco, Lemur catta, Macaca silenus, Pan troglodytes, and Pongo borneo, as well as Homo sapiens (Table 1).
Accordingly, a next-generation sequencing (NGS) approach was used to decode the genomes of the 18 newly isolated B. dentium strains. The sequences obtained were then analyzed together with 20 additional publicly available B. dentium genomes listed in Table 1. Genomic sequences of newly isolated B. dentium strains were sequenced to a coverage depth that ranged from 40-fold to 255-fold, which upon assembly resulted in 10 to 37 contigs (Table 1). By using the complete genome of the type strain B. dentium DSM 20436 as the reference sequence (NCBI accession number AP012326.1), we predicted the contig orientation for each draft genome. After the assembly, the genome length of each B. dentium strain was retrieved, resulting in genomes whose size was shown to range from 2,431,695 to 2,697,950 bp (Table 1). Overall, the number of predicted open reading frames (ORFs) ranged from 1936 for B. dentium VB29 to 2371 for B. dentium BIOML-A2 (Table 1). Thus, collected general genome features may suggest a noticeable genome variability at first glance, with a variation of more than 200 Kbps and 400 ORFs between genomes. Notably, the genome variability identified between strains was not correlated with the presence of extrachromosomal DNA, since no plasmid sequences were identified among B. dentium chromosomes.

Pangenome and Core Genome of the B. dentium Species
Comparative genome analysis of the B. dentium species was performed employing identified genes of each collected strain. Accordingly, the pangenome analysis of the B. dentium taxon was undertaken using a previously described method based on the generation of clusters of orthologous groups (COGs) [65]. The latter analysis allowed the identification of 5754 COGs, representing the pangenome of the B. dentium species. Accordingly, 1287 COGs of the B. dentium pangenome were shared among the 38 genomes of this taxon, thus representing the so-called core genome of the B. dentium species corresponding to 22.4% of its pangenome (Figure 1). shared among the 38 genomes of this taxon, thus representing the so-called core genome of the B. dentium species corresponding to 22.4% of its pangenome (Figure 1). Figure 1. B. dentium pangenome and core genome: (a) displays the number of core genes (green), unique genes (yellow), and dispensable genes (blue) identified in the pangenome analysis; (b) shows the pangenome size based on the sequential addition of the 38 B. dentium genomes; (c) illustrates the core genome size based on the sequential addition of the same genomes. Moreover, genes shared between two or more strains, i.e., dispensable genes, and genes present in just one of the analyzed strains, i.e., truly unique genes (TUGs), were also identified. This analysis allowed us to detect a number of TUGs that ranged from 19 for B. dentium BIOML-A3 to 219 for B. dentium 793B, with an average of 60 TUGs per genome (Table 1). Thus, the average TUG number highlighted that the variability of this taxon was sufficiently explored since the B. dentium ORFeome contains a small number of TUGs when compared to the previously analyzed bifidobacterial pangenomes of B. adolescentis and B. bifidum [27,28].
The pangenome size, when plotted versus the 38 B. dentium genomes, shows that the power trendline tends to reach a plateau. As reported in Figure 1, the genomic data from the last strain added to the pangenome analysis slightly expand the total gene pool. Thus, in accordance to the reported data, the resulting pangenome curve suggests a "closed" B. dentium pangenome. More specifically, after the addition of the 38th B. dentium genome, any additional genome included in the analysis will result in the increases of the pangenome size of about 60 COGs (Figure 1).
The extended core genome coupled with the relatively small number of TUGs and the prediction of a "closed" pangenome, suggested that this taxon does not possess a high genome variability in respect to other bifidobacterial species previously scrutinized, such as Bifidobacterium asteroides, B. longum, and B. pseudolongum [6,30,66]. Instead, the variability between B. dentium strains resulted in being comparable with that previously inspected in members of the B. bifidum and Bifidobacterium breve species [27,67].

Phylogenetic Analyses and Evolutionary Development of the B. dentium Taxon
The assessment of the B. dentium phylogeny was performed applying a previously described methodology [58,68,69] by using data collected from the comparative genomic analysis reported above. Thus, gene sequences of the 38 B. dentium strains were employed as well as those of Bifidobacterium moukalabense DSM 27321, which served as an outgroup. Then, paralogs were excluded from the 1287 B. dentium core genes identified in the core genome analysis using the PGAP pipeline (see Materials and Methods), resulting in 1140 genes of which the concatenated amino acid sequences were employed to build a phylogenetic tree (Figure 2). The resulting B. dentium-based phylogenetic tree showed the presence of three clusters named BdA, BdB, and BdC, composed by 14, 8, and 16 strains, respectively, which seem to have a very close phylogenetic relatedness (Figure 2).
B. dentium strains that were isolated from human beings did not fit into a specific cluster, as well as strains isolated from the other species of primates. Conversely, B. dentium strains isolated from different mammals were found evenly distributed among all clusters. Thus, the reconstructed B. dentium-based phylogenetic tree highlighted that none of the 38 strains considered in this study presented divergent genetic signatures among its core genome sequence. Unfortunately, few genome sequences belonging to strains isolated from the oral cavity were publicly available. So, all four B. dentium retrieved from this ecological niche included in the phylogenetic tree represented the type strain of this species (Figure 2).
To further explore the overall genomic variability among B. dentium, we employed the average nucleotide identity (ANI) approach between all collected genome sequence pairs. Accordingly, we highlighted the genome synteny among members of this taxon, with associated ANI values ranging from 98.2% to 99.9% (Supplementary Materials Table S1). Notably, two strains that display an ANI value <95% may be considered to belong to two distinct species [70]. Based on this notion, in our previous Bifidobacterium based phylogenomic studies, we identified a suitable ANI threshold of about 94% to discriminate between bifidobacterial species [6,68,69]. Thus, ANI values above 98.2% obtained from this phylogenomic analysis highlighted a low variability between genomes of members of the B. dentium taxon. B. dentium strains that were isolated from human beings did not fit into a specific cluster, as well as strains isolated from the other species of primates. Conversely, B. dentium strains isolated from different mammals were found evenly distributed among all clusters. Thus, the reconstructed B. dentium-based phylogenetic tree highlighted that none of the 38 strains considered in this study

Carbohydrate-Active Enzymes among Members of the B. dentium Species
The complete carbohydrate-active enzyme repertoire of the B. dentium taxon was explored to unveil predicted genes encoding glycosyl hydrolases (GHs) and glycosyltransferases (GTs) of each analyzed strain. The screening allowed us to identify members of 29 GH families and 12 GT families, highlighting a predominance of genes predicted to encode GHs belonging to the GH13, GH43, GH3, and GH2 families, and GTs of the GT2 and GT4 families (Figure 3). Together with GH5, GH36, GH31, GH26, and other minority families, which were also identified in the glycobiome of each B. dentium, these GH families represent the wide core glycobiome of the B. dentium taxon involving 19 different GH families. Notably, the most abundant GH family identified in B. dentium, represented by GH13, is also widespread in bacteria and it is characterized in the degradation of a wide range of carbohydrates, including plant-derived polysaccharides, such as amylose, amylopectin, maltodextrins, and pullulan. The degradative abilities of this taxon related to plant-derived polysaccharides were also confirmed by the other most abundant GHs, i.e., GH43 and GH3, involved in the degradation of complex plant polysaccharides, such as (arabino)xylan, (arabino)galactan, and other arabinans. Furthermore, the high number of GH2 family enzymes were represented by β-galactosidases and β-mannosidase, highlighting an advanced ability into the degradation of simple carbohydrates as well. Nonetheless, genes encoding additional 15 GH families are shared among each analyzed B. dentium strain, highlighting a massive degradative ability towards a wide range of carbohydrates ( Figure 3). Additionally, each analyzed genome revealed a conserved distribution of genes encoding for carbohydrate esterases (CEs), while no polysaccharide lyases (PLs) were identified. More specifically, members of four CE families were identified in each B. dentium genome, i.e., CE1, CE2, CE4, and CE9, showing a distribution of one to three family genes per genome.
To identify differences between GH profiles of each B. dentium strain, hierarchical clustering was performed using the complete repertoire of GHs and GTs identified. As disclosed in Figure 3, the resulting clusters did not show any appreciable differences in the GH and GT family profiles. Furthermore, no correspondence was found between phylogenetic clustering (BdA, BdB, and BdC clusters of Figure 2) and clustering based on GHs or GTs families ( Figure 3). Remarkably, B. dentium BRDF23 strain was reported as the outlier of both carbohydrate-active enzymes clustering, but any peculiar gene encoding carbohydrate-active enzymes arise from the analysis.
In order to compare the amount of genetic material involved in the carbohydrate metabolism, the GH count of each B. dentium strain was normalized with the number of predicted genes. The latter process generated GH indexes ranging from 3.59% in B. dentium 831B to 4.24% in B. dentium VB29, highlighting a conserved GH distribution among strains. As previously observed, B. dentium DSM 20436 possesses the highest GH index of the genus Bifidobacterium together with other species isolated from the gut of different primates [68,71]. Thus, the vast repertoire of identified GHs, corroborate with the ability of members of this taxon to establish a relationship with the host in multiple compartments, such as the oral cavity rich in simple sugars and the gut, where the host has not processed complex carbohydrates.
The carbohydrate fermentation capabilities of the B. dentium taxon were then investigated through growth experiments involving newly B. dentium isolated strains together with the type strain of this species. To obtain a complete as possible overview into the B. dentium metabolic abilities, we included both simple and complex plant-derived glycans as the unique carbon source on a semisynthetic medium. As displayed in Figure S1, each B. dentium strain grew on common simple sugars, such as glucose and galactose. Furthermore, all strains could metabolize a broader array of complex sugars, such as maltodextrin, starch, and xylan. However, analyzed strains displayed little if any growth on arabinogalactan, arabinose, and N-acetyl-D-glucosamine ( Figure S1). As previously observed through the in silico carbohydrate-active enzyme profiling, each B. dentium strains' growth performance resulted in a conserved pattern between glycans, except for growths on mannose ( Figure S1). To identify differences between GH profiles of each B. dentium strain, hierarchical clustering was performed using the complete repertoire of GHs and GTs identified. As disclosed in Figure 3, the

B. dentium Mobilome
To investigate the amount of putative horizontal gene transfer (HGT) events that occurred among the B. dentium taxon, the genome of the 38 strains were screened using the software Colombo [57,58].
The analysis resulted in identifying an average of 8.4% putative genes involved in HGT events with respect to the total amount of B. dentium predicted genes (Table 1). Thus, members of this species were found to be less suitable to acquire alien genes with respect to the same analysis recently applied to members of the B. pseudolongum subsp. globosum taxon [30]. This characteristic was also reflected by the low number of TUGs previously identified ( Table 1).
As previously shown through comparative genomic analysis of numerous species belonging to the Bifidobacterium genus, prophage sequences and other elements involved in the transposition of genes, such as transposases, represent a large part of the genetic material implicated in HGT events [58]. Thus, we performed a prophage profiling taking into account the pangenome of the B. dentium taxon, by using a manually curated database composed by genes of prophages identified among the Bifidobacterium genus [59]. Thus, 75 prophage sequences were identified, ranging from zero to three between B. dentium chromosomes (Table 1). In this context, retrieved prophage sequences were predicted to belong to the Siphoviridae family, encompassing modules that putatively encode functions involved in lysogeny, DNA replication, DNA packaging, head and tail morphogenesis, and host lysis [72]. Nonetheless, an additional classification of the latter sequences using the PHASTER tool, classified them as incomplete prophages. Nevertheless, each reported phage sequences displayed several representative genes of the modules listed above, such as integrases, capsid, portal and terminase proteins, tape measure proteins, and lysin and holing proteins (Table S2).
A further investigation into the B. dentium mobilome involved a screening of insertion sequence (IS) elements employing the IS Finder online tool. Overall, IS21 and IS3 were those IS families retrieved from the proteome of each strain, highlighting a few occurrences of one or two IS elements per genome, while in half of these genomes, no IS elements were identified (Table 1). Nevertheless, the B. dentium DE-29 strain was the only exception, displaying four IS elements in its genome (Table 1). Furthermore, across the analyzed genomes, many additional genes were predicted to encode transposases, but we observed a high frequency of partial IS elements that are no longer functional, showing an evolutionary development of this taxon aimed to avoid genome rearrangements.

B. dentium Defense Mechanisms
To enhance their fitness, (bifido)bacteria evolved to contrast HGT events through the ability to defend themselves against the invasion of foreign DNA [73]. To explore the defense mechanism of each B. dentium strain, we investigated the presence of CRISPR-Cas systems, as previously reported for each type-strains of the genus Bifidobacterium [74]. Among the 38 genomes analyzed, one or two CRISPR-Cas systems were identified in each B. dentium strain, except for the 2091B strain. Altogether, 46 CRISPR-Cas systems were shared among the 38 B. dentium chromosomes, highlighting a high occurrence of type I systems (I-C and I-U) characterized by the cas3 gene, and type II systems (II-C) by the cas9 gene (Table 1 and Table S3). Based on our screening, type III systems characterized by the cas10 gene appear to be absent in our assessed strain collection. Interestingly, type II systems were found in 11 B. dentium strains in possession by another type I system. Thus, among the B. dentium taxon, we observed a high degree of conservation of type I systems employed by this taxon to defend themselves from phage invasions.
Another defense mechanism complex that prevents the acquisition of foreign DNA is represented by restriction-modification (RM) systems [75]. The screening revealed that the type IV RM system represented the predominant RM system among the B. dentium taxon, being present in 31 out of 38 genomes, followed by the type II RM system identified in the genomes of 22 strains (Table 1 and  Table S4). Furthermore, type I and III RM systems were identified at a lower frequency, i.e., in 16 and 4 strains, respectively (Table S4). Additional RM genes were retrieved from this analysis but were not identified in specific gene clusters and were excluded from the classification (Table S4). Thus, members of the B. dentium species encoded an average of two putative functional RM systems aiming to prevent the invasion by foreign DNA sequences. In this context, only B. dentium 70T2 lacked in RM systems (Table 1). Altogether, these findings revealed conserved defense mechanisms in almost the entire B. dentium collection, proving a direct correlation between a reduced mobilome and an efficient defense against the invasion of foreign DNA. Correlation between phylogenetic grouping, as reported in Figure 2, and mobilome/defense mechanisms profiles of each strain were not supported by any statistical significance. Thus, note that members of the B. dentium taxon share a close phylogenetic relatedness with low variability between genomes.
Predicted protein-coding genes were then investigated for the presence of antibiotic resistance (AR) elements. However, putative AR genes encoding transporters were omitted from this analysis due to the inaccuracy of their bioinformatic prediction [76,77]. Thus, false positives were excluded through manual editing of the results based on the protein domains identified within predicted AR genes. The overall number of AR genes was of two putative AR enzymes in each of the 38 B. dentium genomes acting against beta-lactams and tetracyclines. More specifically, each B. dentium possesses in its core genome a transglycosylase characterized by a penicillin-binding domain, and a translation elongation factor with a TetM-like domain. In the same fashion, predicted protein-coding genes were screened for the presence of bacteriocins able to inhibit the growth of other bacteria. However, the analyzed proteome of the B. dentium taxon lacked in these peptidic toxins. Furthermore, the presence of toxin-antitoxin systems was also investigated. The screening highlighted 177 putative toxin-antitoxin (TA) protein-encoding genes shared among the 38 B. dentium chromosomes (Table S5). The large majority of the putative TA genes were identified to be widespread across the B. dentium genomes, while 60 were closely linked, representing actual type II TA systems. The most conserved system was represented by the HipAB system, identified in 25 out of 38 B. dentium strains (Table S5). Additional TA systems were represented by three RelBE systems identified in B. dentium BIOML-A1 and B. dentium BIOML-A2, and two HicAB systems identified in B. dentium 793B and B. dentium 2078B (Table S5).

Conclusions
Since their discovery, members of the B. dentium species have been associated with the development of plaque and dental caries in humans [20][21][22][23]. In the last few years, however, a wider distribution of this taxon in the fecal samples of several mammals has been shown [26]. Based on this notion, through culture dependent isolation attempts, we reconstructed the genome sequences of 17 B. dentium strains identified in the gut microbiota of primates and one single strain from a fecal sample of Ursus arctos. Comparative genomic analyses and phylogenetic reconstruction of the taxon revealed that it does not possess a high genome variability regarding other bifidobacterial species, highlighting a very close phylogenetic relatedness between strains. In order to better understand this notion, we explored the mobilome and the defense mechanisms shared between members of the B. dentium species, resulting in identification of a few transposable elements within their genome coupled with consistent machinery to defend against the invasion of foreign DNA. Furthermore, the profiling of the carbohydrate-active enzyme repertoire of each analyzed strain highlighted a taxon that is evolved with a consistent wide core glycobiome aimed to degrade an equally wide range of carbohydrates.