The Integrated Probiotic Database: a genomic compendium of bifidobacterial health-promoting strains

Genome sequences of bifidobacterial commercial strains

In accordance with scientific literature surveys, publicly available chromosomal sequences of 34 bifidobacterial strains used in commercial dietary probiotic products were retrieved from the National Center for Biotechnology Information (NCBI) public database [Table 1]. To ensure a consistent genomic analysis, open reading frames (ORFs) from each bifidobacterial genome sequence were re-predicted and annotated using the most recent release of the MEGAnnotator pipeline^[9]. In detail, contigs greater than 1000 bp were employed to predict protein-encoding ORFs through Prodigal v2.0 (Linux command line “./prodigal -f gff -a [protein_translation_to_selected_file] -i [input_filename.fasta] -o [output_filename]”)^[10]. 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^[11] together with hidden Markov model profile (HMM) searches (http://hmmer.org/) against the manually curated Pfam-A database (cutoff e-value of 1 × 10^-10). Then, tRNA genes were detected through tRNAscan-SE v1.4^[12], while rRNA genes were identified using RNAmmer v1.2^[13].

Comparative genomic analysis

All 34 genome sequences of Bifidobacterium members were employed for a pan-genome analysis using the Pangenome Analysis Pipeline (PGAP) v1.1 (http://pgap.sf.net)^[14]. The predicted proteome of each bifidobacterial genome was classified into functional gene clusters through the gene family (GF) method, consisting of pairwise protein-similarity search using blast software v2.2.28+ (cutoff e-value of 1 × 10^-10 and exhibiting at least 50 % identity across at least 80 % of both protein sequences). The obtained data were used to assign proteins to so-called clusters of orthologous groups (COGs) employing MCL (graph theory-based Markov clustering algorithm)^[15]. A pan-genome profile was then built using a presence/absence matrix encompassing all COGs identified in the analyzed genomes (Linux command line “./PGAP.pl --strains [input_strain_list] --input input_path/ --output output_path/ --thread 20 --identity 0.5 --coverage 0.8 --cluster --method GF --evolution --pangenome”). Subsequently, the core genome of commercial bifidobacterial strains was obtained by selecting the protein families shared between all genomes, while truly unique genes (TUGs) of a given genome were identified based on those protein families that are not present in other bifidobacterial chromosomes. Functional annotation of each TUG arsenal was accomplished employing the eggNOG database^[16]. Each pairwise average nucleotide identity (ANI) was calculated using the program fastANI^[17].

Phylogenomic analysis

To disentangle the phylogenetic relationships between the 34 collected bifidobacterial probiotic strains, the concatenated sequence of amino acids belonging to the core genome of each bifidobacterial strain was aligned using the MAFFT software^[18]. The resulting phylogenetic tree was built using the neighbor-joining method through the ClustalW v2.1 program^[19], and the graphical viewer of phylogenetic trees FigTree v1.4 (http://tree.bio.ed.ac.uk/software/figtree/) was used to its visual representation.

Glycobiome prediction and identification of genes conferring antimicrobial resistance

The genome sequences of the publicly available 34 bifidobacterial probiotic strains were subjected to assessment of genes encoding for glycosyl hydrolase (GH), glycosyl transferase (GT), and polysaccharides lyase (PL) enzymes through sequence similarity search in the carbohydrate-active enzyme (CAZy) database^[20] using HMMER v3.3^[21] (cutoff e-value of 1 × 10^-15) and BLASTP analysis^[22] (cutoff e-value of 1 × 10^-10).

The proteome of each bifidobacterial probiotics genome was also screened for the presence of bacterial antimicrobial resistance based on sequence similarity to genes classified in the Comprehensive Antibiotic Resistance Database^[23] (BLASTP cutoff e-value of 1 × 10^-5). Outcomes were then manually validated to eliminate possible false positives. Moreover, the Transporter Classification DataBase (TCDB)^[24] was employed to assess the putative transporter specificity.

Identification of sortase-dependent pilus-encoding loci, and bacteriocins-encoding genes

Sortase-dependent (SD) pilus-encoding loci (type I and II pili) were identified through homology search tool RAPsearch (cutoff E value of 1 × 10^-5 with minimal alignment length 20)^[25] exploiting the custom sortase-dependent pilus genes database previously built^[26]. Then, a detailed manual inspection was performed to identify complete pilus gene clusters.

Likewise, bacteriocin-encoding genes were detected using RAPsearch analysis (cutoff E value of 1 × 10^-5 with minimal alignment length 20) employing the BAGEL4 database^[27].

Assessment the genetic background for exopolysaccharides, virulence, and bile salt hydrolases production

To identify the loci encoding exopolysaccharides (EPSs), the protein sequences of well-known priming glycosyltransferases (pGTFs) were retrieved from NCBI database and were used to inspect the 34 bifidobacterial genome sequences. Subsequently, for each bifidobacterial chromosome, the genomic regions flanking the putative pGTF were investigated to identify EPS-encoding key genes (such as glycosyltransferases, flippases, ABC transporters, and carbohydrate precursor biosynthesis/modification enzymes). The presence of putative virulence genes and bile salt hydrolases were identified through sequence similarity (homology) search in the Virulence Factor Database^[28] and in the protein sequence NCBI database, respectively (cutoff E value of 1 × 10^-5). Thus, the resulting hits were manually inspected to remove false positives.

Statistical analyses

All statistical analyses were computed using SPSS software (www.ibm.com/software/it/analytics/spss/).

The IPDB

Bifidobacterial strain names from labels of commercially available probiotic products were identified based on comprehensive scientific literature research, and all associated publicly available genomes (complete and draft) were retrieved from NCBI [Table 1]. As a result, we collected a total of 34 bifidobacterial probiotics, including 4 B. longum subsp. infantis, 10 B. longum subsp. longum, 4 B. bifidum, 3 B. breve, and 13 B. animalis subsp. lactis chromosomes sequences constituting the IPDB in its first iteration. Notably, to ensure consistency in the gene prediction, all bifidobacterial genomes used in this study were re-annotated using the MEGAnnotator pipeline as described in the Methods Section^[9]. Subsequently, the 34 bifidobacterial commercial probiotic genomes were employed to perform a comparative genome analysis to identify peculiar genetic traits possibly involved in intestinal colonization and host-microbe interaction.

All the re-annotated genome sequences, along with strain-specific functional details and information concerning the comparative analysis results, are included in the newly developed IPDB available at http://probiogenomics.unipr.it/cmu/ (direct download at http://probiogenomics.unipr.it/files/Pro-biotic_Bifidobacteria_DataBase.zip). Note that IPDB will be expanded to include the genomes of non-bifidobacterial commercialized probiotic strains in the near future.

General genome features of the bifidobacterial strains encompassed in the IPDB

According to the genome prediction and annotation processes, we identified a number of predicted ORFs ranging from 2567 for B. longum subsp. infantis EVC001 to 1565 for B. animalis subsp. lactis BLC1 [Table 1]. As previously reported, B. longum subsp. infantis showed the largest genomes among the probiotic collection, ranging between 2.83 and 2.61 Mb^[29], while B. animalis subsp. lactis resulted in the taxon with the smallest genome sizes (average of 1.95 Mb).

Notably, the ANI investigation highlighted a higher degree of genome identity among the 13 strains belonging to the B. animalis subsp. lactis species used as probiotics (average of 99.8%), compared to all the other considered (sub)species (average of 98.1%) [Supplementary Table 1]. Although the high degrees of synteny and sequence homology between members of this taxon is well-known (ANI ~99.7%)^[30], 53% of the B. animalis subsp. lactis strains showed ANI ≥ 99.99%, indicating that, presumably, identical strains have been effectively deposited and commercialized with different strain names. Moreover, according to the ANI analysis, the strain B. longum subsp. longum 35624, previously misclassified as a member of the B. longum subsp. infantis taxon, is still promoted commercially with an incorrect classification [Supplementary Table 1].

Overview of the commercial probiotics pan-genome

The genome sequences of the 34 bifidobacterial probiotic strains were used to predict five (sub)species-specific pan-genome profiles by classifying each strain-specific proteome into protein families named COGs [Figure 1A]. Combining the obtained five (sub)species-specific pangenomes, we identified the core genome of the bifidobacterial probiotics (BPBs-CG) by taking into account a total of 657 COGs shared by all collected bifidobacterial (sub)species [Supplementary Table 2]. Similarly, five (sub)species-specific core genomes were obtained considering the COGs shared by all the strains belonging to a given sub(species) while being absent in the others. Accordingly, these latter were characterized by 150 B. longum subsp. infantis (Binf-CG), 90 B. longum subsp. longum (Blon-CG), 343 B. bifidum (Bbif-CG), 169 B. breve (Bbre-CG), and 445 B. animalis subsp. lactis (Blac-CG) (sub)species-specific core genes (SSCore genes) [Figure 1A]. Notably, B. longum subsp. longum showed the fewest SSCore genes (ANOVA P-value < 0.05), suggesting that the evolutionary dynamics of this taxon have not led to achieving substantially unique genetic traits, while, in contrast, the phylogenetically correlated subspecies B. longum subsp. infantis showed a marked SSCore comparable to the B. breve species [Figure 1A and B]. Conversely, the relatively high number of SSCore genes of Blac-CG and Bbif-CG could reflect the evolutionary distance between B. animalis subsp. lactis and B. bifidum with respect to the other taxa included in this study, as pointed out by the phylogenetic reconstruction based on BPBs-CG [Figure 1A and B].

Figure 1. Pangenome of the five bifidobacterial (sub)species and phylogenetic relationships reconstruction. (A) The five (sub)species-specific pangenomes profiles. The core gene pools characterizing each bifidobacterial (sub)species, i.e., Binf-CG, Blon-CG, Bbif-CG, Bbre-CG, and Blac-CG, are highlighted as part of each core genome. (B) The phylogenomic tree based on the BPBs-CG describing the phylogenetic relationships among the 34 collected bifidobacterial probiotics. Each (sub)species-based cluster is highlighted with a different color.

The pan-genome analysis also revealed the strain-specific genes repertoire, i.e., TUGs, highlighting a variable number of TUGs ranging from 403 to 12 (average of 122.6 TUGs per genome) [Figure 1A and Supplementary Table 2]. Notably, based on eggNOG analysis, bifidobacterial TUG arsenals included an average of 38.3% of genes with general or unknown function (R/S); an average of 16.5%, 13%, and 9% of genes predicted to be involved in DNA replication (M), carbohydrate (G), and amino acid (E) metabolisms, respectively; and the remaining 23.2% were related to cell wall/membrane biogenesis (M), defense mechanism (V), translation (J), transcription (K), and inorganic ion transport (P) [Supplementary Table 2]. Interestingly, B. longum subsp. infantis showed the highest number of TUGs (average of 322) [Supplementary Table 2]. This observation indicated peculiar features that may characterize the B. longum subsp. infantis strains employed as commercial probiotics^[29,31]. Indeed, the relatively high degree of B. longum subsp. infantis genotype variation could be associated with the high rate of horizontal gene transfer events previously observed within this taxon^[29]. In contrast, B. animalis subsp. lactis exhibited the fewest TUGs (average of 62.7) [Supplementary Table 2], corroborating the limited genetic variability among members of this taxon^[30], as evidenced by the abovementioned ANI analysis.

Distribution of host-derived glycans metabolizing capabilities providing probiotic properties

Probiotic strains can metabolize the complex dietary carbohydrates that cannot be processed by host enzymes through the production of specific GHs, enhancing digestion and conferring health benefits to the host by releasing health-promoting compounds, such as Short-Chain Fatty Acids^[32]. With the aim to investigate the differences in carbohydrate metabolizing capabilities of bifidobacterial probiotics, we explored the metabolic enzyme arsenal for complex carbohydrates, i.e., the glycobiome, catalyzing the breakdown of both dietary and host-derived carbohydrates. For each bifidobacterial probiotic strain, the complete glycobiome profile, including GHs, GTs, and PLs, is reported in Supplementary Table 3.

Based on the CAZy database^[20], we identified about 120 GHs per genome, corresponding to an average of 40.2 different GH families. In particular, 22 of the latter, including enzymes deputed to plant-derived carbohydrates metabolism as well as GH families active on glycosidic linkages of lactose, resulted to be included in the BPBs-CG, thus shared by all bifidobacterial probiotics [Supplementary Figure 1 and Supplementary Table 3].

Remarkably, additional enzymes degrading host-derived glycan structures (HMOs and intestinal mucin) such as GH101 (endo-α-N-acetylgalactosaminidase), GH20 (β-hexosaminidase), GH33 (sialidase), and GH129 (α-N-acetylgalactosaminidase) (www.cazy.org/) were detected in all bifidobacterial probiotics, except B. animalis subsp. lactis, and 27% of the B. longum subsp. longum strains. Consequently, these data highlight strain-dependent abilities of B. longum subsp. longum to digest HMOs-derived structures and, thus, to promote the absorption of nutrients during infant breastfeeding [Supplementary Figure 1 and Supplementary Table 3]. Furthermore, the GH29 family (α-L-fucosidases) was observed to be highly represented in B. bifidum and B. longum subsp. infantis chromosomes [Supplementary Figure 1 and Supplementary Table 3], while GH84 (exo-/endo-β-N-acetylglucosaminidases) and GH89 (extracellular soluble α-N-acetylglucosaminidases) were found exclusively within Bbif-CG, reflecting expanded metabolic capabilities toward host-derived glycan utilization of the abovementioned taxa, compared to the other Bifidobacterium probiotic (sub)species^[33-35].

Interestingly, members of the recently discovered GH136 family, which exert the role of extracellular lacto-N-biosidase^[36], beyond being shared by all B. bifidum probiotics, were found in 63% of those belonging to B. longum subsp. longum [Supplementary Figure 1 and Supplementary Table 3]. This observation might reveal a crucial survival strategy adopted by specific B. longum subsp. longum strains to increase their competitiveness in the infant gut ecosystem, although this subspecies is also adapted to utilize plant-derived oligosaccharides present in the adult diet^[37].

Overall, the genomes of B. longum subsp. longum showed the highest number of accessory GH genes [Figure 2 and Supplementary Table 3]. Indeed, within the chromosomes of this taxon, 24 GH families were found in 9%-90% of the strains, in comparison of only 2-8 GH families constituting the accessory GH arsenal of the other considered probiotic (sub)species [Figure 2 and Supplementary Table 3]. In particular, GH families involved in the degradation of HMOs and host glycan structures, i.e., GH129, GH136, GH85 (endo-β-N-acetylglucosaminidase), and GH29, were found in, respectively, 72.7%, 63.6%, 45.5%, and 9.1% of the probiotic strains belonging to B. longum subsp. longum.

Figure 2. Accessory glycosyl hydrolase (GH) profiles of the bifidobacterial probiotics. For each bifidobacterial probiotic strain, the occurrence of the accessory GHs, i.e., GHs shared by a subset of the considered probiotic strains, are depicted through a bar-plot graph. The accessory GH families active on the host-derived glycans mentioned in the text are highlighted with different colors.

Although carbohydrate utilization capabilities are often associated with (sub)species-specific features, IPDB analyses reported differences in carbohydrate-metabolizing enzymes between commercialized probiotics of the same (sub)species. Such differences can have functional and ecological implications worthy of consideration for probiotic formulation and consumption.

Extracellular structures involved in microbe-host interactions

Bacterial extracellular appendages, such as pili or fimbriae, are long and non-flagellar structures strategically localized to the cell surface to promote bacterial adhesion in the gut, simultaneously impacting microbe-host dialogue^[38,39]. In the Bifidobacterium genus, SD pili (types I and II), collectively representing the SD fimbriome, as well as type IV pili, have been previously described^[26]. While these latter are highly conserved among bifidobacterial genomes, the SD pili showed a considerable variability^[40]. Thus, we explored the SD pili-encoding genes arsenal of the 34 collected bifidobacterial probiotic strains exploiting a custom database built in the contest of a previous study^[26].

Overall, SD pilus gene clusters, composed of a sortase-encoding gene for assembling pilus subunits and two pilin subunit-encoding genes, were found in 91% of the inspecting genomes [Figure 3 and Supplementary Table 4]. Interestingly, while genome sequences of B. longum susp. infantis appear unable to encode this type of pili, probiotic strains belonging to B. bifidum possessed the highest number of SD pili-encoding clusters [Figure 3 and Supplementary Table 4]. In particular, the genomes of B. bifidum strains BGN4 and PRL2010 showed three SD pili loci, thus suggesting putative improved adherence and persistence features. Furthermore, a diverse array of genes required for the production of SD pili was observed between probiotic strains belonging to the same (sub)species. In particular, B. longum subsp. longum showed a variable number of SD pili, ranging 0-2 [Figure 3 and Supplementary Table 4], highlighting possible different abilities to colonize and persist in the human gastrointestinal tract.

Figure 3. Occurrence of genetic probiotic features in the bifidobacterial strains. For each of the 34 considered bifidobacterial probiotic strain, the heat map shows the predicted number of glycosyl hydrolase enzymes involved in host glycan metabolisms, antimicrobial resistance determinants, pili- and bacteriocins-encoding genes, bile salt hydrolases, and exopolysaccharides (EPSs)-encoding loci.

Overall, these data collected in the IPDB, in addition to (sub)species-specific features, highlight considerable strain-level variabilities in the environment interaction structures that could therefore determine different individual probiotic properties.

Production of bacteriocins by commercial bifidobacterial probiotics

In addition to external structures, bifidobacteria exploit molecule-based systems to compete for intestinal colonization directly. Although the inhibitory activity of bifidobacteria could partially derive from the production of organic acids, it is hypothesized that some members of the Bifidobacterium genus can produce antimicrobial molecules such as bacteriocins^[41,42]. These latter are ribosomally synthesized peptides with antimicrobial activities against other bacteria, either belonging to the same species or even across genera^[43,44]. Consequently, these compounds are regarded as a probiotic trait contributing to higher niche competitiveness and inhibition of intestinal pathogens^[45]. For this reason, we investigated the occurrence of bacteriocins-encoding genes among the 34 bifidobacterial probiotics using the BAGEL4 database^[27].

As a result, five potential bacteriocin genes were predicted to be codified by B. longum subsp. infantis and B. longum subsp. longum probiotic strains. In particular, a Class I lantibiotic (BLD_1648) was found in B. longum subsp. infantis EVC001 and in four members of B. longum subsp. longum taxon, i.e., strains CECT7347, CECT7894, DLBL07, and DLBL09 [Figure 3 and Supplementary Table 4].

Based on these in silico analyses results, only a limited number of Bifidobacterium species encode for bacteriocins, and intra-(sub)species variabilities have been found when comparing different strains. In particular, only (certain) strains belonging to B. longum subsp. longum and B. longum subsp. infantis showed strain-specific abilities in producing antimicrobial compounds, which may facilitate the introduction of the (probiotic) producer into an established niche by directly inhibiting competing strains or pathogens. Thus, these findings evidenced by analysis of the IPDB reinforce the need for a precise assessment of desirable probiotic properties, such as bacteriocins production, at a strain-specific level.

Antibiotics resistance prediction and their distribution among commercial bifidobacterial probiotics

Probiotics are specifically selected not to carry AMR, with particular attention to AMR determinants located in the proximity of transposable elements or falling inside (integrated) bacterial plasmids, which could contribute to the spread of AMR^[46,47]. Notably, AMR determinant surveys across the Bifidobacterium genus revealed that, except for tetracyclines resistance (tet genes) in specific cases, the resistance phenotypes are independent of the presence of particular genes, or they do not fall in genomic regions involved in horizontal gene transfer events. Hence, they rarely represent a risk for transfer to unrelated pathogenic or potentially pathogenic bacteria^[48,49]. Conversely, AMR can enhance the survival of the probiotics in the presence of antimicrobial compounds due to medical treatments, thus constituting a beneficial feature^[50].

In this contest, the collected 34 bifidobacterial probiotic strains were inspected for putative antibiotic resistance determinants. Even if our in-silico analysis remains only predictive, such an approach can provide indications for further in vitro validations. As a result, an average of 7.6 AMR genes per chromosome were identified [Supplementary Table 5]. Among these, a putative ATP-binding cassette (ABC) transporter that exports macrolides (ARO:3000535) and putative rifampicin (ARO:3004480), fosfomycin (ARO:3003785), and mupirocin (ARO:3003730) resistances were found shared by all probiotics, while a gene conferring resistance to cationic antimicrobial peptides (ARO:3003577) was shared by 76% of bifidobacterial probiotics strains [Figure 3 and Supplementary Table 5]. Notably, prediction of the putative transporter specificity was assessed with manual validation employing the TCDB^[24].

In addition, (sub)species-specific AMR genes were also observed, including genes putatively involved in resistance against several different classes of antibiotics, i.e., multidrug efflux transporter, within Blon-CG (ARO:3000816) and Bbre-CG (ARO:3002813), and genes conferring putative resistance to mycinamicin (ARO:3001301) and tetracycline antibiotics (ARO:3003980, ARO:3000194) in Blac-CG [Figure 3 and Supplementary Table 5].

Focusing on the unique gene pools, which could be horizontally acquired, strain-specific AMR genes were found in 4 out of the 34 screened probiotic strains [Figure 3 and Supplementary Table 5]. In particular, B. animalis subsp. lactis BS05, showing the lowest ANI value, i.e., 99.3%, when compared to the other genomes of the same species, was predicted to encode a mosaic tetracycline resistance gene (tetW/N/W, ARO:3004442) and a gene possibly involved in resistances to carbapenems, rifamycin, and peptide antibiotics (ARO:3005059) [Figure 3 and Supplementary Table 5]. Moreover, two different macrolide resistance systems (ARO:3000616 and ARO:3004626) were noticed in the B. longum subsp. longum CECT7894 and JDM301 strains, respectively, whereas a generic antibiotic efflux pump (ARO:3000838) was observed in the strain DLBL09 [Figure 3 and Supplementary Table 5].

Based on the data collected in the IPDB, bifidobacterial probiotics appear to possess a relatively low acquired resistance compared to members of the Enterococcus and Lactobacillus genus used as probiotics in humans and farm animals^[51], for which resistance to a wide range of antimicrobials carried on plasmids or in the proximity of conjugative transposons has been identified^[52-54]. Nevertheless, strain-specific AMR determinants have been observed, highlighting the need for case-by-case assessments.

Screening of additional genetic features involved in colonization and persistence

To further characterize the 34 bifidobacterial probiotic strains included in the database, attributes including bile salt tolerance mediated by bile salt hydrolases, production of exopolysaccharides (EPSs), and the presence of putative virulence factors were investigated. Notably, the ability to hydrolyze bile salts is often regarded as a desirable feature for probiotic strain selection since it can promote probiotic fitness and colonization by detoxifying bile^[55]. According to in silico analyses, bile salt hydrolase activity has been predicted for all 34 bifidobacterial strains [Supplementary Table 6], resulting to be a widespread trait among the bifidobacterial probiotics. Furthermore, screening of potential virulence-related genes revealed the presence of homologous genes associated to surface carbohydrates polymers and response regulator proteins which typically mediate the interaction with the surrounding environment [Supplementary Table 7]. However, such structures are not recognized as harmful. Instead, they are well-known to participate in the host-microbe dialogue underlying and supporting the claimed probiotic effects^[56]. Consistently, the analysis revealed the absence of genes associated with clear detrimental effects, remarking the safe use of bifidobacterial strains as probiotics.

Among the interesting and attractive characteristics of probiotic strains, the production of EPSs has grasped the attention because of its important role in maintaining commensalism between human host and (bifido)bacteria as well as for their putative health-promoting properties^[57,58]. EPSs are extracellular carbohydrate polymers, and, for their biosynthesis, a gene cluster including a pGTF and additional genes, such as ABC transporters, subunit polymerization enzymes, and carbohydrate precursor biosynthesis/modification enzymes, are required^[59]. In particular, pGTF is an essential enzyme that catalyzes the first step of the EPS biosynthetic pathway^[59].

In this context, the 34 bifidobacterial probiotics were explored for EPS loci employing well-known pGTF gene sequences as molecular indicators, as previously performed^[60,61].

Accordingly, the production of EPSs was predicted in all bifidobacterial chromosomes except for those belonging to B. bifidum species. Specifically, the presence of two highly conserved EPS loci were observed in each B. animalis subsp. lactis probiotic [Supplementary Table 6], while a single EPS-producing locus with a significant intra-(sub)species variability was detected among probiotic strains belonging to B. longum subsp. longum and B. longum subsp. infantis taxa [Supplementary Table 6]. The precise location of the pGTFs predicted in each bifidobacterial genome is reported in Supplementary Table 6.