Variation in the Conservation of Species-Specific Gene Sets for HMO Degradation and Its Effects on HMO Utilization in Bifidobacteria

Human milk provides essential nutrients for infants but also consists of human milk oligosaccharides (HMOs), which are resistant to digestion by the infant. Bifidobacteria are among the first colonizers, providing various health benefits for the host. This is largely facilitated by their ability to efficiently metabolize HMOs in a species-specific way. Nevertheless, these abilities can vary significantly by strain, and our understanding of the mechanisms applied by different strains from the same species remains incomplete. Therefore, we assessed the effects of strain-level genomic variation in HMO utilization genes on growth on HMOs in 130 strains from 10 species of human associated bifidobacteria. Our findings highlight the extent of genetic diversity between strains of the same species and demonstrate the effects on species-specific HMO utilization, which in most species is largely retained through the conservation of a core set of genes or the presence of redundant pathways. These data will help to refine our understanding of the genetic factors that contribute to the persistence of individual strains and will provide a better mechanistic rationale for the development and optimization of new early-life microbiota-modulating products to improve infant health.


Introduction
Aberrant human gut microbiome development during infancy is linked to various health conditions that can manifest themselves throughout adulthood.The colonization of pioneer microbes shortly after birth represents a key first step in this mutualistic relationship, shaping the developing microbiota and impacting numerous host physiological processes [1].This early-life developmental window represents a critical time for microbehost interactions, as the foundations for future health are established here [2].Breastfeeding is considered one of the most efficient strategies to shape a healthy gut microbiome during infancy [3]; however, its formation is also influenced by various other external factors such as the delivery mode, antibiotic use, and the environment [4][5][6][7][8][9].Bifidobacteria are among the first colonizers and commonly represent over 70% of the total microbial community during breastfeeding, providing various health benefits for the host [10].
The early colonization by bifidobacteria is largely facilitated by their ability to efficiently metabolize human milk oligosaccharides (HMOs).Despite being the third most abundant solid component in breastmilk, they provide no (direct) nutritional value to the infant [22].Instead, HMOs act as compounds that selectively promote the growth of mostly bifidobacteria; thus, the production of different HMOs has been industrialized to fortify formula milk [23].
The human milk glycans are diverse, and over 100 different HMO structures have been identified.This structural diversity is achieved with five building blocks-glucose, galactose, N-acetylglucosamine, fucose, and sialic acid.The core consists of lactose at the reducing end, often elongated by repeats of β1-3-linked lacto-N-biose (Galβ1-3GlcNAc; LNB), termed type 1 chains, or by β1-3/6-linked N-acetyllactosamine (Galβ1-4GlcNAc; LacNAc), termed type 2 chains.These core structures are frequently decorated with fucose and sialic acid residues via α1-2/3/4 and α2-3/6 linkages, respectively.Despite the structural diversity, a small number of HMO species can represent up to 70% of the total, and typically fucosylated HMOs make up the majority (>40%), unless the milk is derived from non-secretor or Lewis-negative mothers [24,25].
Although the gene repertoires and applied strategies are strongly species-specific, the growth rates, efficiency rates, preferences [27], and ranges of the utilized HMOs can differ among strains of the same species [15,19,20,26,27,31,33,49]. For instance, the fucosylated HMO (FHMO) consumption rates vary significantly between species, and generally all B. infantis, B. bifidum, and B. kashiwanohense strains utilize FHMOs.However, recent genomic analyses have revealed that the gene sets for FHMO degradation are also distributed in less than 10% of B. longum, B. breve, and B. pseudocatenulatum strains [17,21,50], providing strains with the capability to utilize FHMOs, likely through inter-species horizontal gene transfer [17].

Bacterial Strains and Growth Conditions
The 130 Bifidobacterium strains used in this study (Table S2) were previously purchased from the American Type Culture Collection (ATCC), isolated from a commercial probiotic product, or derived from the Chr Hansen Culture Collection (Chr Hansen, Denmark).The latter strains were isolated from feces (infants and adults) and are depicted with a Bif number.The bifidobacteria were routinely grown on Man-Rogosa-Sharpe (MRS; BD Difco, New York, NY, USA) medium, supplemented with 0.05% w/v L-cysteine-HCl (Merck KGaA, Darmstadt, Germany) (MRS-Cys) and incubated at 37 • C. All incubations were performed in airtight AnaeroPack boxes (Mitsubishi Gas Chemical Company, Inc., Tokyo, Japan) with AnaeroGen (Thermo Scientific, Oxoid, MA, USA) sachets added to create an anaerobic atmosphere.

Growth on Carbohydrate Sources
The six HMOs tested (2'-O-fucosyllactose (2 ′ FL), 3-O-fucosyllactose (3FL), lacto-Ntetraose (LNT), lacto-N-neotetraose (LNnT), 3'-O-sialyllactose sodium salt (3 ′ SL), and 6 ′ -Osialyllactose sodium salt (6 ′ SL) were kindly donated by Glycom A/S Denmark, through their scientific donation program.In addition, the 5 HMO building block monosaccharides (glucose, galactose, N-acetylneuraminic acid, L-fucose, N-acetyl-D-glucosamine) and the backbone di-saccharide lactose were all purchased from Merck.For the evaluation of growth on the 12 carbohydrates as the sole carbon source, the strains were cultured in modified MRS (mMRS) medium, which was manually prepared to obtain the following composition: tryptone (10.0 g/L), yeast extract (5 g/L), sodium acetate trihydrate (5 g/L), K 2 HPO 4 (2.6 g/L), di-ammonium hydrogen citrate (2.0 g/L), Tween-80 (1 g/L), L-cysteine-HCL (0.5 g/L), MgSO 4 •7H 2 O (0.2 g/L), MnSO 4 •H 2 O (0.05 g/L).Prior to autoclaving at 121 • C for 15 min, the medium was adjusted to pH 6.9.All carbohydrates were added directly to the mMRS medium (mMRS-HMO) at a final concentration of 2% (w/v), which was subsequently filter-sterilized (0.2 µm).All strains were cultured in biological triplicates.For the growth assays, the strains were first grown overnight at 37 • C in MRS-Cys.Ten-fold dilution series were prepared from the overnight cultures and incubated overnight under the same conditions with fresh medium.For each strain, late exponential-early stationaryphase cells were selected based on their optical density at 600 nm (OD 600nm ).These were washed twice in large volumes of mMRS medium and normalized to OD 600nm = 1.0 in mMRS, prior to inoculation.The growth experiments were performed in triplicate in 96-well microtiter plates with 180 µL of mMRS-HMO and 20 µL of bacterial cell suspension.Furthermore, each experiment contained a well with MRS-Cys as the positive control, mMRS-HMO without bacteria as the negative control, and mMRS medium as the background growth control.The microtiter plates were incubated at 37 • C under anaerobic conditions.The bacterial growth was measured at OD 600nm with a Synergy H1M (Biotek, Winooski, VT, USA) plate reader after 96 h of growth.

Genome Sequencing, Assembly, Annotation, and GH Identification
The bacteria were grown overnight in MRS-Cys and the pellets were dissolved in DNA/RNA Shield (Zymo research Europe, Freiburg, Germany) and subsequently shipped to BaseClear BV, Leiden, The Netherlands, where the paired-end sequencing was performed on Illumina platforms.
The draft genomes were de-novo-assembled with SPAdes v3.12.0 [51] using the '--careful' setting, and predictions of putative open reading frames (ORFs) and annotations based on BLASTP were performed with PROKKA v1.14.5 [52].As Illumina-based genome sequencing does not allow the complete closure of bacterial genomes, we assessed the genome completeness using CheckM2 v1.0.2 [53].The glycosyl hydrolases were predicted and annotated based on similarity to the Carbohydrate-Active enZYmes (CAZy) (http://www.cazy.org/)database [54].The counts of glucosyl hydrolase (GH) families were extracted using custom Python scripts and subsequently reduced to only the GH families associated with HMO degradation.

Homology Searches for HMO Utilization Genes
Queries of nucleotide sequences based on the literature were retrieved from full genome sequences of bifidobacteria deposited in the NCBI reference sequence database, accessed on 1 July 2023 (https://www.ncbi.nlm.nih.gov/), and the JGI repository accessed on 18 May 2022 (https://genome.jgi.doe.gov/portal/).The accession numbers are listed in Table S1.The query nucleotide sequences of the B. scardovii GHs were taken from the corresponding genome.The prevalence and similarity rates of the HMO utilization genes in the 130 draft Bifidobacterium genomes were examined via a BLAST analysis (BLAST + v2.9.0) using the Resfinder engine v3.2 [55], with the following criteria: identity >50%, query coverage > 50%.We ran the analysis with different thresholds to estimate the importance, and going lower or higher did not affect the main results.The heatmaps were generated with R version 4.2.0 [56] using the ComplexHeatmap v2.11.2 package [57].

Genome Sequencing
To identify the HMO-degrading machinery of the bifidobacteria, we performed genome sequencing and assembly and carried out a uniform open reading frame (ORF) prediction of 130 isolates.Illumina-based genome sequencing does not allow the complete closure of bacterial genomes; therefore, we assessed their genome quality using genome completeness (99.9 ± 0.2%), contamination (0.36 ± 0.68%, contig N50 (267k ± 270k nt), and L50 (5.3 ± 2.8) values, which all indicated that the genomes were of high quality and suitable for the subsequent analyses.S2).As an additional quality check based on biological genomic features, we calculated the previous genomic characteristics per species (Tables S3-S5).All metrics were in accordance with the published reference statistics for each species [19,26,[58][59][60][61].

GHs Involved in Degradation of Neutral Type I and II Chain HMOs
A heatmap summarizing all genes found in all genomes is supplied as Supplementary Figure S1.When reporting the statistics and prevalence of the common HMO-degrading genes in this chapter, we excluded B. animalis and B. angulatum, as they are not prototypical human species or do not degrade HMOs.However, their data are supplied as Supplementary Information.
The most prevalent β-N-acetylglucosaminidase (GH20), Blon_0732 from B. infantis [45], was 100% prevalent in B. infantis and B. longum.The version found in 100% of the B. breve strains, Bbr_1556 (NahA), was homologous to a different B. infantis GH20 enzyme, namely Blon_0459.In contrast, B. pseudocatenulatum (100%) and B. catenulatum (56%) shared a unique GH20 gene, BBPC_1688, not found in any other species.The prevalence of β-Nacetylglucosaminidases across all tested infant-derived bifidobacterial species supports their relevance in HMO degradation.However, their overall prevalence rate was lower than for LNT β-1,3-galactosidases, as in both B. pseudocatenulatum and B. catenulatum they were located inside LNB/LNT utilization clusters, which were sometimes entirely absent from the latter species (Figure S1).

Lnp Locus
The Lnp gene cluster is responsible for the import of mainly LNB by the ABC transporter (gltABC) and subsequent phosphorolysis by GH112 GNB/LNB phosphorylase (lnpA), releasing Gal-1-P,GlcNAc/GlcNAc-6-P through the modified Leloir pathway [31,38,49,63,64].The complete lnp locus was prevalent in 100% of the prototypical infant species B. longum, B. infantis, B. bifidum, and B. breve.Interestingly, in B. longum, we found two different versions of this gene cluster, on which we further elaborate in the species-specific chapter.
In contrast, all B. pseudocatenulatum and 4 out of 8 B. catenulatum strains encoded the gltABC transporter, although this was combined with LNB/LNT utilization machinery (nagA, nagB, and a GH20 GH) instead of the conventional lnp cluster layout that combines this transporter with the modified Leloir-like pathway, like the other infant bifidobacteria [65].B. adolescentis did not encode an lnp cluster or gltA (Figure S1).

Fucosyllactose-Fucose Utilization Loci
The Bifidobacterium fucose and FHMO utilization loci identified outside the typical FHMO utilizers B. infantis and B. bifidum are similarly organized.They consist of intracellular fucosidases (AfcA or AfcB), fucose catabolism genes, and an ABC transporter.The SBPs of the ABC transporter responsible for the intake of fucosylated (FL) HMOs are classified into 2 types, FL1-BP and FL2-BP.These are further classified into 4 clades based on their substrate specificity, comprising cluster 1-IV (FL1-BP) and clusters 2-I, 2-II, and 2-III (FL2-BP).These transporters are exclusively distributed in bifidobacterial genomes.Transporters belonging to cluster 1-IV can import 2 ′ FL and 3FL, whereas clusters 2-I and 2-II can import 2 ′ FL, 3FL, lactodifucotetraose (LDFT), and lacto-N-fucopentaose I (LNFP I).Cluster 2-III can take up the widest range of FHMOs-2 ′ FL, 3FL, LDFT, LNFP I, LNFP II, and lacto-N-difucohexaose I/II (LNDFH I/II).The substrate specificity of the FL transporter type in the FHMO cluster corresponds to the presence of corresponding fucosidase genes, i.e., transporters with a narrow substrate specificity are paired with a GH95 AfcA homolog 1,2-α-L-fucosidase, while transporters with wider substrate specificity rates are paired with a GH95 and complementary GH29 AfcB homolog 1,3/4-α-l-fucosidase to cover the larger range of imported FHMOs [19].The additional GH29 comes paired with a L-fucose mutarotase FucU.Genomic analyses have revealed that these gene sets for fucosylated HMO degradation are distributed in almost all strains of B. kashiwanohense, although in less than 10% of strains of B. longum, B. breve, and B. pseudocatenulatum [17,21,50].In this dataset, we identified versions of this complete FHMO cluster in one B. pseudocatenulatum strain but none in 67 B. longum strains.Interestingly, in all 8 B. breve strains we only found partial, non-functional gene sets.

Species-Specific HMO Degradation Patterns and Genetic Machinery
Despite these shared genetic characteristics, most gene sets were species-specific in their homology, content, and organization, exemplified by the B. infantis 43,000 kb H1 HMO cluster, the extracellular B. bifidum GHs, and the B. breve lnt and nah loci.However, we also observed within-species strain-level genetic variation, of which the extent differed between species.The HMO degradation profiles were mostly strain-specific due to the similarity in the utilized HMOs, except for B. infantis and B. bifidum, strains whose degradation patterns were similar within species and highly divergent from other species.Due to the clear species-specific genotype and strain-specific phenotype, we performed species-specific gene trait matching to determine the impact of genetic variation on HMO utilization.

B. breve Results
All B. breve strains grew well on LNT, while 3 out of the 8 strains did not grow on LNnT.All 8 strains contained at least one member of the required GH families for LN(n)T/G degradation (GH2 (lacZ2/Z6), GH42 (nahA), and GH20), as well as one GH112 LNB-GNB phosphorylase, GH33 α-sialidase, and GH95 1,2-α-fucosidase.As the HMO GH profile was conserved in all strains, the reasons for these phenotypical differences must be in the auxiliary HMO degradation machinery.
Both the lnt and lnp−glt loci were extremely well conserved in all strains.However, we observed small variations in the LN(n)T SBP (nahS) and transcriptional response regulator nahR from the nah locus.One group encoded a nahS gene with 97% identity versus 99% identity, as compared to the gene from strain UCC2003.This different SBP correlated with the lower overall growth on LNnT.Additionally, we observed some small changes in the nahR transcriptional regulator, although this did not seem to have a large effect on the final OD 600 value; however, it might affect the response time and growth rate in vivo.
Despite the presence of fucosidases, no strains grew on 2 ′ FL/3FL.All tested strains encoded the complete fucose catabolism machinery and fucosidases yet interestingly lacked the complete FL ABC transporter.To better evaluate the prevalence of this genotype, we additionally analyzed 7 B. breve genomes from our database that were not part of the HMO screening presented here.These strains displayed the same genetic configuration.Lastly, all strains encoded a GH33 α-sialidase but did not grow on 3/6 ′ SL (Figure 1).
cluster in one B. pseudocatenulatum strain but none in 67 B. longum strains.Interestingly, in all 8 B. breve strains we only found partial, non-functional gene sets.

Species-Specific HMO Degradation Patterns and Genetic Machinery
Despite these shared genetic characteristics, most gene sets were species-specific in their homology, content, and organization, exemplified by the B. infantis 43,000 kb H1 HMO cluster, the extracellular B. bifidum GHs, and the B. breve lnt and nah loci.However, we also observed within-species strain-level genetic variation, of which the extent differed between species.The HMO degradation profiles were mostly strain-specific due to the similarity in the utilized HMOs, except for B. infantis and B. bifidum, strains whose degradation patterns were similar within species and highly divergent from other species.Due to the clear species-specific genotype and strain-specific phenotype, we performed species-specific gene trait matching to determine the impact of genetic variation on HMO utilization.

Results
All B. breve strains grew well on LNT, while 3 out of the 8 strains did not grow on LNnT.All 8 strains contained at least one member of the required GH families for LN(n)T/G degradation (GH2 (lacZ2/Z6), GH42 (nahA), and GH20), as well as one GH112 LNB-GNB phosphorylase, GH33 α-sialidase, and GH95 1,2-α-fucosidase.As the HMO GH profile was conserved in all strains, the reasons for these phenotypical differences must be in the auxiliary HMO degradation machinery.
Both the lnt and lnp−glt loci were extremely well conserved in all strains.However, we observed small variations in the LN(n)T SBP (nahS) and transcriptional response regulator nahR from the nah locus.One group encoded a nahS gene with 97% identity versus 99% identity, as compared to the gene from strain UCC2003.This different SBP correlated with the lower overall growth on LNnT.Additionally, we observed some small changes in the nahR transcriptional regulator, although this did not seem to have a large effect on the final OD600 value; however, it might affect the response time and growth rate in vivo.
Despite the presence of fucosidases, no strains grew on 2′FL/3FL.All tested strains encoded the complete fucose catabolism machinery and fucosidases yet interestingly lacked the complete FL ABC transporter.To better evaluate the prevalence of this genotype, we additionally analyzed 7 B. breve genomes from our database that were not part of the HMO screening presented here.These strains displayed the same genetic configuration.Lastly, all strains encoded a GH33 α-sialidase but did not grow on 3/6′SL (Figure 1).
The lack of growth on LN(n)T by 3 of the included 8 strains was discordant with a group of 24 strains that all showed good growth on both substrates [20].However, inconsistencies in the growth phenotypes of strains have been observed previously, such as growth [20] and lack of growth [17] on LNnT by the type strain DSM20213 T .Nevertheless, all tested strains contained at least one member of the required GH families for LN(n)T degradation and one GH112 LNB-GNB phosphorylase, GH33 α-sialidase, and GH95 1,2-α-fucosidase, in concordance with these being reported as part of the B. breve core glycobiome [58].B. breve is known to possess broad but variable saccharolytic capacity, although the GHs specific for HMO degradation seem to be well conserved.
One strain encoded an extra GH29 α-1-3/4-fucosidase.This gene was also identified in 4 out of 24 strains, of which 15 were isolated from 3-4 month-old exclusively breastfed term infants from the USA [20], and 0 out of 13 isolates from Japanese infants [21], making this a fairly rare genotype.Despite the presence of fucosidase, no strains grew on 2 ′ FL/3FL.Nevertheless, this phenotype has been reported for B. breve strains encoding the FHMO locus with an FL2-II/AfcA and FL1-IV/AfcA configuration [21].All tested strains (and 7 of the additionally analyzed genomes) encoded the complete fucose catabolism machinery and fucosidases yet interestingly lacked the complete FL ABC transporter.The presence of a strain with the additional GH29 AfcB combined with a FucU gene hints at B. breve strains with wider HMO assimilation phenotypes that might still be discovered.Interestingly, the presence of this partial cluster in our dataset is in contrast to other strains where either the whole cluster was present or absent [17,19,50].
As the B. breve strains with FL-SBs were all previously isolated from Japanese infants from 2 studies [21,66], we wondered whether a 2 ′ FL-consuming strain isolated from a North American infant (SC95) would also display this genomic feature.However, we found the same partial transporter-deficient genotype.This discrepancy still remains unexplained, although the lack of growth on 3FL for this strain, which should be degraded by the additional GH29 α-1-3/4-fucosidase in its genome, might hint at the use of a novel transporter that only imports 2 ′ FL.
Although all strains encoded a GH33 α-sialidase, none grew on 3/6 ′ SL, which is a common feature, as B. breve strains have been shown to consume acidic HMOs, although they prefer more complex structures such as LSTb and monosialyllacto-N-hexaose (S-LNH) [20].Similarly, one strain that did not grow on 2 ′ FL did utilize larger fucosylated HMOs from a complex human milk mixture [20].This also points to the previously mentioned presence of novel undiscovered transporters or yet unknown SBP specificity for B. breve strains.In conclusion, the pan genome based on 73 B. breve strains has been shown to be open and approaching saturation.The observed genetic diversity included genes involving various dietary-or host-derived carbohydrate utilization capabilities [58].This allows B. breve to utilize dietary substrates associated with infancy up to weaning (HMOs, lactose, and plant-derived carbohydrate sources).Additionally, the tightly controlled system for the transcriptional regulation of HMO metabolism allows B. breve strains to rapidly switch their metabolic processing to and from HMOs, lactose, and plant-derived carbohydrate sources, which is a regular occurrence during weaning.This variability in dietary carbohydrates usage seems to extend to the HMO consumption profiles in the case of B. breve.The cores of the transporters and GHs are sufficient for the degradation of LN(n)T; therefore, the majority of strains grow well on both.However, the apparent variation in the degradation of more complex structures points to the prevalence of yet unknown acidic and fucosylated HMO transporters in combination with a GH29 α-1-3/4 fucosidase.

B. longum subsp. infantis Results
All B. infantis strains grew well on all substrates, although the OD 600 values were generally lower on 3 ′ SL and 6 ′ SL.In our set of 8 strains, two out of three GH20 β-Nacetylglucosaminidases were conserved in all strains, namely Blon_2355 located in H1 and Blon_0732 located outside of H1.Blon_0459 was missing from 2 strains.Both the GH42 LNT β-1,3-galactosidase Bga42A/Blon_2016 and GH2 Bga2A/Blon_2334, located in H1, were conserved in all strains.In contrast, the prevalence of other GH2/GH42 βgalactosidases was strain-specific, although these are active on β-1,4-gal in GOS rather than LN(n)T [62,67].Four out of the five previously reported fucosidases in the type strain ATCC15697 T were conserved in all strains, whereas the final fucosidase belonging to GH151 (Blon_0346) was missing in three out of six strains.All strains encoded GH33 α-sialidases, NanH2/Blon_2348 located in H1, and NahH1/Blon_0646 located in H4.All strains also contained one GH112 GNB/LNB phosphorylase located in the lnp (H5) cluster.
All strains contained at least one GH for the internal sequential degradation of both type I and type II cores, FHMOs, and acidic HMOs; therefore, the specialization towards specific HMOs in our set is not based on GH profile alone but rather the complement of SPBs can predict the spectrum of available substrates and is responsible for phenotypical differences in vitro (and likely in vivo) (Figure 2).
controlled system for the transcriptional regulation of HMO metabolism allows B. breve strains to rapidly switch their metabolic processing to and from HMOs, lactose, and plantderived carbohydrate sources, which is a regular occurrence during weaning.This variability in dietary carbohydrates usage seems to extend to the HMO consumption profiles in the case of B. breve.The cores of the transporters and GHs are sufficient for the degradation of LN(n)T; therefore, the majority of strains grow well on both.However, the apparent variation in the degradation of more complex structures points to the prevalence of yet unknown acidic and fucosylated HMO transporters in combination with a GH29 α-1-3/4 fucosidase.

Results
All B. infantis strains grew well on all substrates, although the OD600 values were generally lower on 3′SL and 6′SL.In our set of 8 strains, two out of three GH20 β-Nacetylglucosaminidases were conserved in all strains, namely Blon_2355 located in H1 and Blon_0732 located outside of H1.Blon_0459 was missing from 2 strains.Both the GH42 LNT β-1,3-galactosidase Bga42A/Blon_2016 and GH2 Bga2A/Blon_2334, located in H1, were conserved in all strains.In contrast, the prevalence of other GH2/GH42 βgalactosidases was strain-specific, although these are active on β-1,4-gal in GOS rather than LN(n)T [62,67].Four out of the five previously reported fucosidases in the type strain ATCC15697 T were conserved in all strains, whereas the final fucosidase belonging to GH151 (Blon_0346) was missing in three out of six strains.All strains encoded GH33 αsialidases, NanH2/Blon_2348 located in H1, and NahH1/Blon_0646 located in H4.All strains also contained one GH112 GNB/LNB phosphorylase located in the lnp (H5) cluster.
All strains contained at least one GH for the internal sequential degradation of both type I and type II cores, FHMOs, and acidic HMOs; therefore, the specialization towards specific HMOs in our set is not based on GH profile alone but rather the complement of SPBs can predict the spectrum of available substrates and is responsible for phenotypical differences in vitro (and likely in vivo) (Figure 2).Considering the complete HMO degradation machinery, we could broadly divide our strains into four genetic configurations based on the presence or absence of (partial) genetic loci.Group I had retained the most complete utilization machinery (M-63, EVC001, ATCC15697 T ).The three other groups missed the complete lnp/H5 transporter (gltABC).Group II additionally missed nahS (Bif175, Bif181) and group III missed both Considering the complete HMO degradation machinery, we could broadly divide our strains into four genetic configurations based on the presence or absence of (partial) genetic loci.Group I had retained the most complete utilization machinery (M-63, EVC001, ATCC15697 T ).The three other groups missed the complete lnp/H5 transporter (gltABC).Group II additionally missed nahS (Bif175, Bif181) and group III missed both nahS and fucosyllactose cluster 3 (Bi-26).Lastly, group IV missed fucosyllactose cluster 3 and partially H2 (NLS superstrain, Bifin02).
The growth profiles of our strains confirmed the previous work, showing that the degradation capabilities of individual HMOs were comparable (i.e., all strains grow well on all tested HMOs), although small differences have been reported, such as one out of 21 strains displaying poor growth on 6 ′ SL [30].Nevertheless, at a deeper level, strains have been shown to differ in both their substrate efficiency (max OD 600 ) and growth rates on individual HMOs [27,33,66].
The aforementioned different configurations of HMO degradation gene sets suggest specialization towards specific or more broad HMO categories.Some strains have lost the unique gltA from lnp/H5 (Blon_2177) for the import of larger type 1 HMOs [27,69], whose absence likely renders these strains incapable of internalizing these HMOs whole.In contrast, none of our tested strains (or other strains) have been shown to miss both FL1-II BP and FL1-IV BP, suggesting that the competitive advantage of FHMO utilization is indispensable for B. infantis and is not exchanged for the broader range of complex HMOs that gltA offers.Therefore, strains within the subspecies B. infantis may be broadly divided into generalists or those focused more on the utilization of smaller (fucosylated) HMOs.
The observed configuration (SBPs missing from H1 and the absence of partial or complete genetic loci) is concordant with research from 2010, which used B. infantis strains isolated before 1990, before this subspecies was commercially produced as a probiotic.Furthermore, the absence of the gltABC transporter was observed in <25% of 21 publicly available B. infantis genomes [49], as well as in a recently isolated native strain from Bangladesh [74].This strongly suggests that this feature is prevalent in natural B. infantis populations, rather than the result of gene reduction in the production of probiotics, as previously suggested [26].
Interestingly, these 4 genotypical configurations coincide with 4 groups that were based on the total B. infantis genome size.The group that included the type strain had an average genome size of 2796 kb, while another group including Bi-26 displayed a 7.3% smaller average genome size (2606 kb).Part of that genome size difference was attributed to selected gene loss in Bi-26, largely of various transporters and glycosidases, which likely reflects adaptation to the utilization of specific carbohydrates [27].However, Bifin02, with a similar configuration of its HMO clusters as Bi-26, has a genome size of 2758 kb (closed circular genome), indicating that other genes are also responsible for these genome size differences.
Recently, regulatory mechanisms for HMO degradation in the generalist (gltA+) B. infantis ATCC15697 T were identified.GlcNAc and its phosphorylated derivatives, GlcNAc-6P and GlcNAc-1P, were identified as potential nagR transcriptional effectors, suggesting that the release of GlcNAc during the internal degradation of any GlcNAc-containing glycans by GHs (e.g., LN(n)T, and in particular fucosylated HMOs (e.g., LNFP I) and milk N-glycans), results in the upregulation of nagR-controlled genes.These de-repress the nag and H1 cluster, including genes encoding LNT and LNnT transporters, as well as the lnp/H5 cluster, explaining the similarity of the transcriptomic responses of B. infantis to LNT and LNnT in vitro [30,33].In contrast, in response to fucosyllactoses (2 ′ FL, 3FL, and DFL), all top upregulated genes occurred in alternative operons outside H1 (fucosyllactose cluster 1-3) [30,70], indicating that dedicated pathways for fucosyllactose utilization have evolved in B. infantis.
Therefore, some strains seem to have evolved to use certain pathways more efficiently, both via the gene deletion of gltABC, as observed in our dataset, and by differing in their global regulatory responses to HMO structures (suggesting changes in global regulatory response networks) [27,69].This has previously been shown in the group III strain Bi-26 (missing gltABC, nahS, and fucosyllactose cluster 3).Compared to the generalist type strain, small fucosylated HMOs exerted far-reaching regulatory roles comparable in effect size to that of lactose-the preferred carbon source of B. longum [75] and potentially one of its global regulators [76,77].In contrast, the type strain perturbed the expression of a much smaller complement of genes.This led to higher growth rates of Bi-26, potentially providing a competitive advantage over other B. infantis strains in infants breastfed with milk high in 2 ′ FL or 3FL.Another way to gain a competitive advantage would be to employ different sequential preference utilization profiles.
Taken together, these findings demonstrate major strain-specific adaptations to the efficient utilization of FLs via two routes: a different configuration of the HMO degradation machinery caused by the gene deletion of the transporters GltABC and nahS to restrict the uptake range of (neutral) HMOs and deletion of the redundant FHMO transporters (FL1-IV) with lower substrate specificity; a change in global regulatory response networks to different HMOs.Other configurations might lead to varying efficiencies (more biomass from the same substrate), response times, growth rates, and sequential preferences for other single HMOs.Conversely, the group I strains, including the type strain, have maintained all genes needed for the metabolism of all HMOs, as well as genes more likely related to mucin usage.These differences also result in altered metabolic outputs or shifts in the ratios of metabolic end-products, including metabolites used in cross-feeding with other microbes, such as 1,2-propanediol and fucose [78], which eventually would impact ecosystem dynamics.
While most of the functionality is retained through genetic redundancy, the genetic specialization of some B. infantis strains (and the complex NagR regulon structure in the generalist B. infantis ATCC15697 T ) likely reflects the evolutionary adaptation to the simultaneous utilization of multiple HMOs.This suggests that the use of rationally formulated HMO mixtures rather than individual oligosaccharides as prebiotics may be a more efficient solution for the selective stimulation of B. infantis growth in the neonatal gut, as it considers the nuanced regulatory mechanisms and physiology of the target organism [65].

B. longum subsp. longum Results
Almost all of the tested 67 B. longum strains grew moderately to well on LNT only, whereas a few strains showed minor growth on LNnT.We found copies of GH42 LNT-β-1,3galactosidase (BLNG_01753 [31], a homolog of Bga42a/lntA) and GH20 β-N-acetylglucosaminidase [79] (a homolog of Blon_0732) in all strains.However, BLNG_00015 [31], a homolog of GH2 β-1,4-galactosidase Bga2A, was only present in 51 out of 67 strains (76%).Therefore, all strains encoded the full complement of genes to internally degrade LN(n)T, as Bga42a/lntA is also effective on β-1,4-linked-gal.Despite all strains encoding one GH112 GNB/LNB phosphorylase [38], we identified two different versions of the gene-one related to BLLJ_1623 [73] from the adult isolate and type strain JCM1217 T , while the other was more similar to BLNG_00163 [31], from the infant isolate SC596.Lastly, the B. longumspecific extracellular membrane-bound GH136 lacto-N-biosidase (lnbX) was present in 43 out of 67 (64%) strains, and the co-occurrence rate with the chaperone needed for the proper folding and function of lnbX (lnbY) [80] was 100%.
The prevalence rate of the complete lnp locus was 100%.However, the gltABC transporter was also present in two versions; here, 47 out of 67 (70%) strains encoded a homolog to the adult isolate JCM1217 T (BLLJ_1624-1626), while the rest encoded a homolog to the infant strain SC596 (BLNG_00160-00162).All genomes also contained a well-conserved group of genes homologous to the B. infantis nag cluster, although the additional gltFGH (GNB) transporter was missing in all B. longum strains (Figure 3).
Most of our tested strains grew moderately to well on LNT, while only a few displayed minor growth on LNnT, in line with previous observations [31].None of the tested strains grew on 2 ′ FL/3FL, and as expected we found no homologs to any of the genes for FHMO uptake, hydrolysis, or L-fucose catabolism encoded by the FHMO locus [17,19,21].The prevalence rate of this cluster in B. longum was only 3% in 151 genomes in public repositories [49], so its absence in our set of 67 strains was to be expected.
Most of our tested strains grew moderately to well on LNT, while only a few displayed minor growth on LNnT, in line with previous observations [31].None of the tested strains grew on 2′FL/3FL, and as expected we found no homologs to any of the Although the prevalence rate of the complete lnp locus was 100%, interestingly the gltABC transporter and GH112 GNB/LNB phosphorylase were present in two versions.Here, 47 out of 67 (70%) strains encoded homologs to the adult isolate JCM1217 T (BLLJ_1623-1626), while the rest encoded homologs to the infant strain SC596 (BLNG_00160-00163).Although gltA from JCM1217 T can bind LNT, its binding affinity is approximately 1000x and 100x lower than for GNB and LNB, respectively [73], and it is not considered to import LNT, as the disruption of lnbX in another adult strain with the same gltA (100% nt identity), JCM31944, eliminated the growth on LNT [82].In contrast, gltA from the infant strain SC596 (BLNG_00160) displayed similar affinities for type 1 HMOs as its homolog in B. infantis (high affinity for LNB, GNB, LNT, and other structurally related glycans, such as sialyl LNT, and even several complex bi-and tri-antennary N-glycans, such as those found in milk and intestinal glycoproteins [31]).
All genomes also contained genes homologous to the B. infantis nag cluster, encoding the GlcNAc catabolic pathway (BLNG00457-4060), without the additional gltFGH (GNB) transporter.When present, both the lnp and nag (nagK, nagB, and nagA) clusters are organized in the nagR regulon, as previously mentioned for other Bifidobacterium species.
In the single B. longum strain lnbX+ that was tested, this gene was additionally under the control of nagR [65].
Therefore, the main differences in the LNT-utilizing machinery were in the presence and absence of lnbX and the 2 versions of the lnp cluster.The prevalence of lnbX was 64% in all strains and 55% and 68% in combination with the infant (BLNG_00160) and adult (BLLJ1626) gltA types, respectively.In all cases, the growth of the lnbX+ strains was better than for lnbX-(OD 600 values in all strains: 1.16 vs. 0.89, p = 0.0003; combined with BLNG_00160: 1.22 vs. 0.91, p < 0.001; combined with BLLJ1626: 1.14 vs. 0.88, p = 0.02).Despite this apparent growth advantage in vitro, the actual proportion of secretory LnbXpositive B. longum in an infant population was shown to be only 0.2% of the B. longum total [82,83].This corresponds with the lower percentage of infant-type gltA combined with lnbX.Furthermore, B. longum prioritizes many other sugars such as Gal, GlcNAc, and lactose over GNB or LNB utilization, which suggests that most infant gut B. longum strains prefer transporters for assimilating mainly LNT and prefer simple sugars over LNB.In this context, it is noteworthy that all strains encoded a complementary ABC transporter specific for GNB but not LNB (BLNG_00933-936) [31].The higher prevalence rate of lnbX+ in combination with BLLJ_1626 might exist to compensate for the poor binding affinity for LNT of the latter.Nevertheless, BLLJ_1626+/lnbX-strains still grew to comparable OD 600 values as other lnbX-strains, in contrast to the observation made by Yamada et al. [82], who reported the elimination of growth on LNT after the disruption of lnbX.We have no explanation for this discrepancy, although in the majority of strains, lnbX does not seem indispensable for growth on LNT.Our results, focusing on the HMO-degrading machinery, are in line with previous observations of a higher genetic diversity among B. longum strains compared to other Bifidobacterium species, except for B. adolescentis [15,[58][59][60].In the particular case of B. longum strains, their diversity and capabilities to metabolize a large range of carbohydrates have been suggested to arise from the influence of the intra-individual environment on epigenetic mechanisms, resulting in differential growth rates on carbohydrate substrates as adaptations to dietary changes during the early life developmental window [14], making them particularly powerful microbial competitors in the dynamic and complex human gut environment.We have shown that this diversity extends to the HMO-degrading machinery, although the phenotype was relatively stable due to the presence of complementary or redundant pathways for the degradation of LNT.Despite the finding that the presence of extracellular lnbX resulted in higher growth rates in vitro, thereby indicating a potential competitive feature, the percentage of lnbX-positive strains studied earlier seems relatively low at only 38% [49] compared to 64% in this study.However, it is likely that the "selfish" transporter-dependent intracellular digestion strategy enables the bifidobacteria to efficiently capture preferred carbon sources in the competitive ecosystem, providing an advantage over other gut microbes.
All 14 tested strains encoded the aforementioned essential genes at very high homology rates.In accordance with the genetic conservation, the degradation of all tested HMOs was relatively homogeneous.The growth on all individual HMOs was good, with better growth on 3′SL than 6′SL, in concordance with the preference of 2,3 over 3,6linkages by SiaBb2 [46].Nevertheless, a notable exception was the lack of growth on 3/6′SL by JCM1255 T , although an absence of exo-α-sialidase activity has been reported for this strain previously [46].This strains also contains a 17nt frameshift mutation in lnpA [84], which has been suggested to account for the limited growth of the strain on HMOs [85].Indeed, we did observe significantly poorer growth on all other individual HMOs by this strain.Additionally, two other strains did not grow on 2′FL.
We observed the same distinctly high genetic conservation rate of all other genes encoded by the lnp cluster.However, one permease subunit of the ABC transporter (gltB, BBPR_1056) has been shown to be frequently absent (6 out of 10 strains) [48,86], and this genotype was assumed to be associated with reduced growth on mucin [48].However, all 14 strains tested here, including JCM1255 T , encoded BBPR_1056.Moreover, in these 2 studies, BBPR_1056 was shown to be both absent [48] and present [86] in JCM1255 T .We
All 14 tested strains encoded the aforementioned essential genes at very high homology rates.In accordance with the genetic conservation, the degradation of all tested HMOs was relatively homogeneous.The growth on all individual HMOs was good, with better growth on 3 ′ SL than 6 ′ SL, in concordance with the preference of 2,3 over 3,6-linkages by SiaBb2 [46].Nevertheless, a notable exception was the lack of growth on 3/6 ′ SL by JCM1255 T , although an absence of exo-α-sialidase activity has been reported for this strain previously [46].This strains also contains a 17nt frameshift mutation in lnpA [84], which has been suggested to account for the limited growth of the strain on HMOs [85].Indeed, we did observe significantly poorer growth on all other individual HMOs by this strain.Additionally, two other strains did not grow on 2 ′ FL.
We observed the same distinctly high genetic conservation rate of all other genes encoded by the lnp cluster.However, one permease subunit of the ABC transporter (gltB, BBPR_1056) has been shown to be frequently absent (6 out of 10 strains) [48,86], and this genotype was assumed to be associated with reduced growth on mucin [48].However, all 14 strains tested here, including JCM1255 T , encoded BBPR_1056.Moreover, in these 2 studies, BBPR_1056 was shown to be both absent [48] and present [86] in JCM1255 T .We currently do not have an explanation for the significant discrepancies in the conservation rates between these datasets.Notably, it seems unlikely this gene is sensitive to deletion in different working stocks of the type strain, considering the high conservation rates of all other genes involved in host glycan degradation.
We found homologs for all genes from the nagR regulon, which has recently been identified in most human bifidobacteria in the lnp and nag (nagK, nagB and nagA) clusters, in concordance with the high conservation rates of other HMO-degrading machinery.However, due to the specific extracellular strategy of B. bifidum, nagR also regulates multiple genes encoding GHs participating in extracellular HMO (LnbB, LnbX, BbgIII) and mucin O-glycan degradation in strain PRL2010.This transcription factor might, therefore, function as a global regulator of both HMO and mucin-O-glycan degradation in this species [65].
Interesting is the presence of lnbX in this regulon.Although all strains contained a homolog to B. longum GH136 lacto-N-biosidase LnbX and its chaperone LnbY (<50% coverage, 50% identity threshold) and both proteins are expressed in the cell, the gene is thought to be non-functional in B. bifidum, due to missing signal peptides in either LnbY or LnbX [80].
Interestingly, strain Bif176 apparently encoded a homolog to the gltFGH transporter from the B. infantis nag cluster that bound to LNB/GNB/Lewis a and type 1 H-trisaccharides (present in glycolipids and glycoproteins in several cell types, including the intestinal epithelium) [87], as well as a homolog to the GH29 α-1-3/4-fucosidase that was unique to B. scardovii.This genotypical feature may be interesting but apparently it did not lead to discernable changes in HMO utilization within our setup, although it might confer benefits in the competitive environment of the infant gut.
In conclusion, all surveyed genes showed very high prevalence and conservation rates in this species, which were mirrored by a relatively homogeneous HMO degradation phenotype.However, in contrast to other species, very small mutations have been shown to lead to large phenotypical differences, especially for the type strain JCM1255 T , which has a 17nt frameshift mutation in lnpA that (potentially) interferes with LNB utilization.We also observed the significantly lower growth rates of this strain on substrates that produce LNB as intermediate steps during their catabolism.There might also be mutations in its siaBb2, which completely lacks exo-α-sialidase activity [46], although all genes in our set were highly alike, similarly to the GH95 1,2-α-L-fucosidase (AfcA) in strains Bif072 and Bif011, which showed no growth on 2'FL.The main distinction between B. bifidum and other bifidobacteria is its capacity to grow on both HMOs and mucin.This seemingly unique capability, facilitated by its large repertoire of extracellular GHs, may have led to the lack of loss or acquisition of other carbohydrate-metabolizing abilities, suggesting this capability provides a strong selective advantage over other (bifido)bacteria [86].

B. pseudocatenulatum Results
All B. pseudocatenulatum strains grew on LNT only, with one strain growing well on 2 ′ FL and 3FL.We found homologs for both the prototypical β-1,3/1,4-galactosidases Bga42A/lntA (Blon_2016/Bbr_0529) and Bga2A (blon_2334) in all strains and also identified a GH20 β-N-acetylhexoaminidase.Therefore, all five tested strains contained at least one member of the required GH families for LN(n)T degradation (GH2, GH42, GH20).The GH20 β-N-acetylhexoaminidase displayed no significant homology to any of the 'main' Bifidobacterium GH20 β-N-acetylhexoaminidases, and its presence was restricted to the species B. pseudocatenulatum and B. catenulatum.
With regards to the transporters, we identified homologs to structural units of the GNB/galactosyl/lactosamine transporter from B. longum (BLNG_0933-0936, missing the ATP-binding module) and of both B. breve nah and lnt transporters, missing both SBPs for the uptake of LN(n)T and LNT, respectively.In contrast, B. pseudocatenulatum contained the gltABC transporter, which was combined with LNB/LNT utilization machinery (nagA, nagB, and the GH20 β-N-acetylhexoaminidase).Instead of the conventional lnp cluster layout that combines this transporter with the modified Leloir-like pathway for processing LNB/GNB, type I HMOs and glycolipids, and intestinal glycans, like the other infant bifidobacteria [38,73,88], we additionally found a homolog to the nagK from B. breve outside this cluster.We confirmed the presence of a complete FHMO cluster in the type strain DSM20438 T , which contains a GH95 α-1,2 fucosidase, which was more similar to genes from other FHMO clusters [17,21,50] than Blon_2335/AfcA in H1 from B. infantis (Figure 5).With regards to the transporters, we identified homologs to structural units of the GNB/galactosyl/lactosamine transporter from B. longum (BLNG_0933-0936, missing the ATP-binding module) and of both B. breve nah and lnt transporters, missing both SBPs for the uptake of LN(n)T and LNT, respectively.In contrast, B. pseudocatenulatum contained the gltABC transporter, which was combined with LNB/LNT utilization machinery (nagA, nagB, and the GH20 β-N-acetylhexoaminidase).Instead of the conventional lnp cluster layout that combines this transporter with the modified Leloir-like pathway for processing LNB/GNB, type I HMOs and glycolipids, and intestinal glycans, like the other infant bifidobacteria [38,73,88], we additionally found a homolog to the nagK from B. breve outside this cluster.We confirmed the presence of a complete FHMO cluster in the type strain DSM20438 T , which contains a GH95 α-1,2 fucosidase, which was more similar to genes from other FHMO clusters [17,21,50] than Blon_2335/AfcA in H1 from B. infantis (Figure 5).

Discussion
Unlike other Bifidobacterium species more common to breastfed infants such as B. infantis, B. longum, B. breve, or B. bifidum, B. pseudocatenulatum is underexamined, despite its frequent presence in breastfed infant feces [12,[18][19][20][21], although it is also prevalent in adults [12,16,19,89].All 6 tested strains contained at least one member of the required GH families for LN(n)T degradation (GH2, GH42, GH20), although the 100% prevalence of Bga42A/lntA in this strain set is in contrast to the ~30% previously reported for 64 publicly available genomes [49].The B. pseudocatenulatum-and B. catenulatum-specific GH20 β-Nacetylhexoaminidase displayed no significant homology to any of the 'main' Bifidobacterium GH20 β-N-acetylhexoaminidases, although we can assume its function due to the growth of all strains of LNT.
For the uptake of LNT, strains of this species use the B. pseudocatenulatum-and B. catenulatum-specific LNB/LNT utilization cluster (gltABC, nagA, nagB, and the β-Nacetylhexoaminidase, combined with nagK and bga42A), which is under influence of the nagR regulon, resulting in the de-repression of nagR-controlled genes due to the release of GlcNac by the activity of the GH20 β-N-acetylhexoaminidase [46], thereby explaining the upregulation of the LNB/LNT utilization cluster by LNFPI but not 2′FL [19].This machinery explains the growth of all strains on LNT, which is in concordance with previous work [19].Nevertheless, one study reported minimal growth on LNnT for one out of 10 infant isolates [18], which might have been non-specific growth.

Discussion
Unlike other Bifidobacterium species more common to breastfed infants such as B. infantis, B. longum, B. breve, or B. bifidum, B. pseudocatenulatum is underexamined, despite its frequent presence in breastfed infant feces [12,[18][19][20][21], although it is also prevalent in adults [12,16,19,89].All 6 tested strains contained at least one member of the required GH families for LN(n)T degradation (GH2, GH42, GH20), although the 100% prevalence of Bga42A/lntA in this strain set is in contrast to the ~30% previously reported for 64 publicly available genomes [49].The B. pseudocatenulatumand B. catenulatum-specific GH20 β-Nacetylhexoaminidase displayed no significant homology to any of the 'main' Bifidobacterium GH20 β-N-acetylhexoaminidases, although we can assume its function due to the growth of all strains of LNT.
For the uptake of LNT, strains of this species use the B. pseudocatenulatumand B. catenulatum-specific LNB/LNT utilization cluster (gltABC, nagA, nagB, and the β-Nacetylhexoaminidase, combined with nagK and bga42A), which is under influence of the nagR regulon, resulting in the de-repression of nagR-controlled genes due to the release of GlcNac by the activity of the GH20 β-N-acetylhexoaminidase [46], thereby explaining the upregulation of the LNB/LNT utilization cluster by LNFPI but not 2 ′ FL [19].This machinery explains the growth of all strains on LNT, which is in concordance with previous work [19].Nevertheless, one study reported minimal growth on LNnT for one out of 10 infant isolates [18], which might have been non-specific growth.
Additionally, we confirmed the presence of a complete FHMO cluster in the type strain DSM20438 T [19].The catalytic specificity rates of GH95 α-fucosidases differ among Bifidobacterium species, and this 1,2-α-L-fucosidase has shown cross-reactivity on 1,3-linkages in 3FL, as strains of B. pseudocatenulatum and B. breve grew well on both 2 ′ FL and 3FL without a dedicated GH29 α-1,3/4 fucosidase [19].This explains the growth of this strain on both 2 ′ FL and 3FL, although it is only equipped with a GH95 α-fucosidase.Nevertheless, the presence of a GH29 α-1,3/4-fucosidase in this cluster [19] has been shown to promote the cleavage of an additional range of HMO moieties that are poorly cleaved by the GH95 enzyme, suggesting that the addition of the second α-fucosidase expands the pool of fucosylated HMOs that are catabolized.Therefore, similar to other species, although growth experiments with a limited range of HMOs show comparable growth characteristics, the substrate competition in vivo may still show a growth advantage to strains possessing the complementary GH29 α-1,3/4-fucosidase [19].
The prevalence of these complete FHMO utilization loci in B. pseudocatenulatum is estimated to be low (<10%).However, this might be strongly influenced by their origin, as no α-fucosidase-positive strains were identified in 45 strains from Vietnam (isolated from children and adults), while in two 3-5-month-old exclusively breastfed, vaginally born UK infants, 1 infant contained a GH95+ strain.Similarly, 3 out of 6 strains isolated from Japanese vaginally born breastfed infants [21] contained the complete FHMO locus [21], suggesting a strong association of FHMO+ strains with infant ecosystems.Nevertheless, surprisingly for B. longum, strains capable of metabolizing FHMOs were also isolated from a formula-fed baby that only received standard non-supplemented (i.e., no prebiotics or synthetic HMOs) formula [14].
The fact that B. pseudocatenulatum seems to have an open pan-genome, approaching saturation [60], combined with overrepresentation of the FHMO cluster in infants and the presence of B. pseudocatenulatum in adults, suggests genomic flexibility and is indicative of (rapid) adaptation to resource availability.The pan-genome of B. pseudocatenulatum may expand mainly with genes responsible for carbohydrate transport and metabolism through horizontal gene transfer (HGT) with other Bifidobacterium species, thereby expediting the diversification of clonal bacteria.B. pseudocatenulatum populations were shown to undergo a prolonged period of within-host evolution and expansion, leading to a strain-level differential response to a dietary intervention due to differences between strains in gene copy numbers of carbohydrate-active enzymes targeting plant polysaccharides [90].

B. catenulatum Results
Four out of 8 B. catenulatum strains grew moderately on LNT.In all strains, we found homologs for both prototypical β-1,3/1,4-galactosidases, Bga42A/lntA (Blon_2016/Bbr_0529) and Bga2A (Blon_2334).Four strains additionally encoded a GH20 β-N-acetylhexoaminidase related to the B. pseudocatenulatum gene BBPC_1688.In these strains, we also identified homologs to the same complement of transporters as B. pseudocatenulatum, except for BLNG_0933-0935.We also found the LNB/LNT utilization cluster and other genes under the NagR regulon influence [65], as described previously for B. pseudocatenulatum.However, the genes in this cluster are slightly different, while the gltABC transporter has slightly higher similarity rates.Strains with and without the LNB/LNT utilization cluster grew to OD 600 0.70 ± 0.15 vs. 0.32 ± 0.14 (p < 0.001) on LNT, respectively (Figure 6).

Discussion
The presence of a full complement of GHs (GH2, GH42, GH20) for the internal sequential degradation of LN(n)T, combined with the LNB/LNT utilization cluster and other genes under the NagR regulon influence [65], as described above for B. pseudocatenulatum, explains the growth on LNT by strains containing the LNB/LNT utilization cluster, similar to B. pseudocatenulatum.B. catenulatum consists of the subspecies kashiwanohense and catenulatum, of which the latter is mainly found in adults [61,89,91] but also infants [61] and was shown to be shared between an infant-mother pair [92].Previous comparative genomic analyses between the two subspecies kashiwanohense and catenulatum identified the FHMO cluster in all (5) kashiwanohense strains but not (yet) in B. catenulatum (11 strains, including 2 infant isolates).Additionally, the GH reservoir of B. catenulatum was more directed to plant-derived glycans present in the adult intestine.Only 1 infant and 1 adult-derived strain encoded a GH20 gene, giving them the full complement of GHs for the internal sequential degradation of LN(n)T [61], which is in contrast to our findings, where 4 out of 8 strains encoded a GH20 gene.As far as we know, this is the first set of B. catenulatum strains to be tested for growth on HMO, and our results suggest that the infant type B. catenulatum strains encode and express genes closely related to the B. pseudocatenulatum LNB/LNT utilization cluster, enabling the internalization and degradation of LNT, while the adult type has lost this ability, which provides a rationale for the presence and prevalence of this species in both infants and adults.Additionally, although the number of strains was low, the pan-genome of B. catenulatum tended to be more open than B. kashiwanohense, suggesting more flexible environmental adaptability for the former.In contrast, kashiwanohense is specialized to infants, with 100% prevalence of the FHMO cluster so far [61], while this has not been identified yet in B. catenulatum.

Discussion
The presence of a full complement of GHs (GH2, GH42, GH20) for the internal sequential degradation of LN(n)T, combined with the LNB/LNT utilization cluster and other genes under the NagR regulon influence [65], as described above for B. pseudocatenulatum, explains the growth on LNT by strains containing the LNB/LNT utilization cluster, similar to B. pseudocatenulatum.B. catenulatum consists of the subspecies kashiwanohense and catenulatum, of which the latter is mainly found in adults [61,89,91] but also infants [61] and was shown to be shared between an infant-mother pair [92].Previous comparative genomic analyses between the two subspecies kashiwanohense and catenulatum identified the FHMO cluster in all (5) kashiwanohense strains but not (yet) in B. catenulatum (11 strains, including 2 infant isolates).Additionally, the GH reservoir of B. catenulatum was more directed to plant-derived glycans present in the adult intestine.Only 1 infant and 1 adult-derived strain encoded a GH20 gene, giving them the full complement of GHs for the internal sequential degradation of LN(n)T [61], which is in contrast to our findings, where 4 out of 8 strains encoded a GH20 gene.As far as we know, this is the first set of B. catenulatum strains to be tested for growth on HMO, and our results suggest that the infant type B. catenulatum strains encode and express genes closely related to the B. pseudocatenulatum LNB/LNT utilization cluster, enabling the internalization and degradation of LNT, while the adult type has lost this ability, which provides a rationale for the presence and prevalence of this species in both infants and adults.Additionally, although the number of strains was low, the pan-genome of B. catenulatum tended to be more open than B. kashiwanohense, suggesting more flexible environmental adaptability for the former.In contrast, kashiwanohense is specialized to infants, with 100% prevalence of the FHMO cluster so far [61], while this has not been identified yet in B. catenulatum.

Discussion
This explains why none of the strain grew on any of the tested HMOs, except for Bif034 and Bif129 on LNnT.One strain (Bif106) showed potential towards host-derived glycans due to the presence of (a combination of) the B. longum transporter (BLNG_0933-0936) and a homolog of Bga42A/lntA (Blon_2016/Bbr_0529).However, this was not a strain that grew on LNnT, and we currently have no explanation for the latter.These observations largely support a pan-genomic analysis of several B. adolescentis genomes, showing that this species has specialized itself towards the utilization of plant-derived glycans, which are normally present with high abundance in adult diets [59].

Discussion
This explains why none of the strain grew on any of the tested HMOs, except for Bif034 and Bif129 on LNnT.One strain (Bif106) showed potential towards host-derived glycans due to the presence of (a combination of) the B. longum transporter (BLNG_0933-0936) and a homolog of Bga42A/lntA (Blon_2016/Bbr_0529).However, this was not a strain that grew on LNnT, and we currently have no explanation for the latter.These observations largely support a pan-genomic analysis of several B. adolescentis genomes, showing that this species has specialized itself towards the utilization of plant-derived glycans, which are normally present with high abundance in adult diets [59].

B. scardovii Results
The B. scardovii strain grew well on LN(n)T and 3FL.We found homologs to the GH42 β-1,3-galactosidase Bga42A/lntA (Blon_2016/Bbr_0529) and GH20 β-N-acetylhexoaminidase (Blon_0732) but not to the prototypical GH2 β-1,4-galactosidase Bga2A/lacZ6 (Blon_2334/ Bbr_1552).Nevertheless, we predicted 6 other GH2 enzymes that were unique for this species, as well as 2 unique GH29 genes.With regards to the transporters, we identified homologs to a partial B. breve lnt cluster (missing lntR), the complete GNB/galactosyl/lactosamine transporter from B. longum (BLNG_0933-0936), and a partial lnp cluster (missing nahK/lnpB).Interestingly, it contains the type II HMO/mucin glycan transporter, similar to the B. infantis H1 cluster hmoA2B2C2 (Blon_2342-4346) [65], which was not found in the B. scardovii type strain [65].Strikingly, it is the only HMO-degrading strain without a homolog to the B. breve nah transporter.The genome also contained a homolog to the partial nag cluster of nagA, nagB, and nagK, giving it a similar nagR regulon structure as B. longum, with hmoA2B2C2 potentially included as well.It also contained homologs to several fucose catabolism genes from B. infantis H1, namely fucD, one ABC transporter permease subunit, and an L-fucose dehydrogenase (Figure 8).
Nutrients 2024, 16, x FOR PEER REVIEW 21 of 27 glycan transporter, similar to the B. infantis H1 cluster hmoA2B2C2 (Blon_2342-4346) [65], which was not found in the B. scardovii type strain [65].Strikingly, it is the only HMOdegrading strain without a homolog to the B. breve nah transporter.The genome also contained a homolog to the partial nag cluster of nagA, nagB, and nagK, giving it a similar nagR regulon structure as B. longum, with hmoA2B2C2 potentially included as well.It also contained homologs to several fucose catabolism genes from B. infantis H1, namely fucD, one ABC transporter permease subunit, and an L-fucose dehydrogenase (Figure 8).

Discussion
This genetic machinery explains the good growth of this strain on LN(n)T because Bga42A/lntA has cross-reactivity towards β-1,4-linkages in LNnT.Most bifidobacteria have a dedicated GH2 β-1,4-galactosidase, and one of the predicted GH2 enzymes unique for this species had no known substrate, meaning it might still encode a specialist enzyme with activity towards the terminal Gal in type II chains.However, we cannot explain the growth on 3FL.Despite the genome encoding the necessary GH29 fucosidases and homologs to several fucose catabolism genes, we did not identify a FL transporter.Although Bifidobacterium species are frequently isolated from the human intestine, the type strain JCM 12489 T was isolated from human blood.B. scardovii can use LNB as a supplementary carbon source [93], although as far as we know, this is the first strain for which growth on HMOs has been tested.

Conclusions
In conclusion, our findings highlight the extent of the genetic diversity between strains of the same species and demonstrate the effects of strain-level genetic diversity on species-specific HMO degradation pathways and utilization.Most strains retained their species-specific HMO degradation profile through the conservation of a core set of genes or through the presence of overlapping or redundant genes or pathways.In particular, the B. infantis strains showed large differences in their configuration of complete HMO utilization loci to potentially alter their HMO utilization preference via the gene deletion of transporters to restrict the range of internalized HMOs, while retaining most of the functionality through genetic redundancy.These different types of genetic diversity within species and the specialization of species, although not visible with growth on the individual HMOs only measured at one timepoint as in this study, likely influence competitive behavior for HMO foraging in situ between different Bifidobacterium strains and species.
Additionally, the presence of a partial FHMO cluster in all analyzed B. breve strains, in contrast to its complete presence or absence in other datasets, combined with one strain containing an additional GH29 α-1-3/4-fucosidase, hints at the potential of B. breve strains

Discussion
This genetic machinery explains the good growth of this strain on LN(n)T because Bga42A/lntA has cross-reactivity towards β-1,4-linkages in LNnT.Most bifidobacteria have a dedicated GH2 β-1,4-galactosidase, and one of the predicted GH2 enzymes unique for this species had no known substrate, meaning it might still encode a specialist enzyme with activity towards the terminal Gal in type II chains.However, we cannot explain the growth on 3FL.Despite the genome encoding the necessary GH29 fucosidases and homologs to several fucose catabolism genes, we did not identify a FL transporter.Although Bifidobacterium species are frequently isolated from the human intestine, the type strain JCM 12489 T was isolated from human blood.B. scardovii can use LNB as a supplementary carbon source [93], although as far as we know, this is the first strain for which growth on HMOs has been tested.

Conclusions
In conclusion, our findings highlight the extent of the genetic diversity between strains of the same species and demonstrate the effects of strain-level genetic diversity on species-specific HMO degradation pathways and utilization.Most strains retained their species-specific HMO degradation profile through the conservation of a core set of genes or through the presence of overlapping or redundant genes or pathways.In particular, the B. infantis strains showed large differences in their configuration of complete HMO utilization loci to potentially alter their HMO utilization preference via the gene deletion of transporters to restrict the range of internalized HMOs, while retaining most of the functionality through genetic redundancy.These different types of genetic diversity within species and the specialization of species, although not visible with growth on the individual HMOs only measured at one timepoint as in this study, likely influence competitive behavior for HMO foraging in situ between different Bifidobacterium strains and species.
Additionally, the presence of a partial FHMO cluster in all analyzed B. breve strains, in contrast to its complete presence or absence in other datasets, combined with one strain containing an additional GH29 α-1-3/4-fucosidase, hints at the potential of B. breve strains with wider HMO assimilation phenotypes that might still be discovered.Further, we complement the little data available for B. catenulatum and B. scardovii by showing that some strains of the former grow on LNT and some of the latter on LN(n)T and 3FL by identifying the utilized gene sets.We also found that B. longum strains encoded two different versions of the lnp cluster in combination with lnbX and determined the effect of this genotype on the growth on LNT.
Exploring the genomic diversity of strains across different species will help to refine our understanding of the genetic factors that contribute to the persistence of individual strains in different ecological contexts, such as the ability to utilize a specific set of HMOs and the formation of cross-feeding networks with other (Bifidobacterium) species.The current methods of microbiome research generally lack the resolution to discern strainspecific differences that shape the complex network of host-microbe interactions in the human gut.Even within the simplest example, the infant gut, where the microbiome's complexity is a fraction of that in an adult, the impact of strain-specific metabolic variations on this ecosystem is not well understood [17,27,33,69].Defining individual strain's roles within the complex system will be essential for understanding the numerous interactions affecting host health throughout life, and could be used as a guide for health-promoting supplementation.Therefore, detailed mechanistic knowledge of the relationships between substrates and specific strains within a species will inform effective and rationalized prebiotic and probiotic strategies.Single or multi-strain synbiotic supplements could be designed to pair specific HMOs with specific strains that will consume them via known strategies to achieve predefined functional outcomes and for selective stimulation in a competitive environment.

Figure 1 .
Figure 1.B. breve, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 2 .
Figure 2. B. infantis, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 2 .
Figure 2. B. infantis, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 3 .
Figure 3. B. longum, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 3 .
Figure 3. B. longum, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 4 .
Figure 4. B. bifidum, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 4 .
Figure 4. B. bifidum, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 5 .
Figure 5. B. pseudocatenulatum, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 5 .
Figure 5. B. pseudocatenulatum, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 6 .
Figure 6.B. catenulatum, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 6 .
Figure 6.B. catenulatum, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Nutrients 2024 , 27 (
16,  x FOR PEER REVIEW 20 of Blon_2334) and 5/13 Bga42A/lntA (Blon_2016/Bbr_0529)), contrary to all other species tested.None of the strains encoded α-sialidases, while four of the strains contained a GH95 α-fucosidase.With regards to the transporters, we found homologs to the complete GNB/galactosyl/lactosamine transporter from B. longum (BLNG_0933-0936) and structural units of the B. breve nah transporter, although nahS (the LN(n)T SBP) was missing.No homologs to any lnp clusters or the B. pseudocatenulatum-B. catenulatum LNB/GNB/LNT catabolic cluster was found, nor were genes under the influence of the NagR regulon, as found in HMO-utilizing bifidobacteria[65] (Figure7).

Figure 7 .
Figure 7. B. adolescentis, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 7 .
Figure 7. B. adolescentis, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 8 .
Figure 8. B. scardovii, from left to right: heatmap of growth on HMOs and building blocks (OD600); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.

Figure 8 .
Figure 8. B. scardovii, from left to right: heatmap of growth on HMOs and building blocks (OD 600 ); number of predicted GH genes and prevalence and sequence conservation rates of HMO utilization genes.The query gene clusters and predicted global function are color-coded on top, while the locus tag and protein name are shown at the bottom.
The numbers of predicted ORFs ranged from 1550 for a B. animalis to 2607 for a B. infantis.The average G + C content was 59.8 ± 1.6%, and the values ranged from 55.9% for a B. catenulatum to 64.8% for B. scardovii.The average genome size was 2.32 ± 0.19 Mbp, whereby B. scardovii had the largest genome (3.13 Mbp) and B. animalis the smallest (1.92 Mbp) (Table