Human milk oligosaccharide metabolism and antibiotic resistance in early gut colonizers: insights from bifidobacteria and lactobacilli in the maternal-infant microbiome (2026)

KEYWORDS:

• Human milk
• oligosaccharides
• antibiotic resistance
• bifidobacteria
• infant
• mother

Bacterial growth on a pool of human milk oligosaccharides

The ability of the isolated bacterial strains to utilize a pool of HMOs was assessed as previously described.^Citation27 Briefly, pooled HMOs (pHMOs) were extracted and purified following the method outlined by Jantscher-Krenn et al.^Citation28 Human milk samples, donated by healthy volunteers to a milk bank, were pooled and centrifuged to remove the lipid layer. Cold ethanol was added to precipitate proteins, which were then removed by centrifugation. Ethanol was eliminated via rotary evaporation. Lactose (Lac) and salts were removed using size exclusion and fast protein liquid chromatography. The resulting oligosaccharide mixture was lyophilized and stored at −20 °C until use.

Strains were cultured overnight in modified peptone-yeast-glucose broth (mPYG, Supplementary S1) under anaerobic conditions (5% H₂, 10% CO₂, balanced with N₂) in a chamber (Whitley A45, Don Whitley Scientific, West Yorkshire, UK) to obtain the seed cultures for the assays. Assays were performed in diluted complex media, rather than minimal media, so that both positive and negative growth effects could be compared to a baseline level of growth. For anaerobic conditions, the media were degassed prior to inoculation.

Strains were grown under 4 different conditions: 1) mPYG dextrose-free medium, as the negative control (mPYG); 2) mPYG medium (mPYG_gluc); 3) mPYG dextrose-free medium supplemented with 15 g/L of pHMOs (pHMOs); 4) mPYG medium supplemented with 15 g/L of pHMOs (pHMOs_gluc). Each medium condition was inoculated with 5% (v/v) seed culture in biological triplicate in 96-well plates and overlaid with 50 μL sterile mineral oil to prevent evaporation. Plates were transferred into a BioTek Epoch 2 microplate spectrophotometer (Agilent, Santa Clara, CA, USA) for overnight growth measurements, which was placed under anaerobic conditions (5% H₂, 10% CO₂, balanced with N₂). Plates were incubated for 24–48 h at 37 °C while optical density (OD₆₀₀) readings were recorded every 30 min with 10s double orbital shaking immediately prior to each reading.

Individual HMOs and by-products utilization assays of bacterial strains

The growth of selected strains on individual HMOs or their degradation by-products was established as previously described,^{Citation29,Citation30} using the mPYG medium. For Lactobacillaceae strains, N-acetyllactosamine (LacNAc; Elicityl, Crolles, France) and lacto-N-neotetraose (LNnT; Glycom, Esbjerg, Denmark) were added individually to the mPYG dextrose-free medium at a concentration of 4 mm. For Bifidobacterium strains, lacto-N-tetraose (LNT; Glycom, Esbjerg, Denmark) and LNnT were added to the mPYG dextrose-free medium at a concentration of 2 mm. Glucose (Sigma-Aldrich, St. Louis, MO, USA), Lac (PanReac AppliChem, Barcelona, Spain) and galactose (PanReac AppliChem, Barcelona, Spain) were used as positive controls of growth, and mPYG dextrose-free medium was used as a negative control. Briefly, the strains were cultured overnight in mPYG medium and diluted to an OD₆₀₀ of 0.1 in 200 μl of mPYG dextrose-free medium supplemented with each carbohydrate in 96-well plates. Cultures were grown at 37 ºC under the anaerobic conditions described above. Bacterial growth was determined by measuring the OD₆₀₀ with Cerillo StratusTM Microplate Reader (Cerillo, Charlottesville, VI, USA) every 30 min with a constant shaking of 80 rpm for 48 h. At least two independent biological replicates and three technical replicates were performed for each growth assay.

HMO utilization by each bacterial strain

The utilization of HMOs was estimated by measuring the HMO profile of the bacterial growth supernatants. Bacterial cultures were centrifuged 10 min at 5,000 rpm and supernatants were collected and stored at −80 ºC. Each supernatant (20 μL) was spiked with maltose (20 mg/mL) as an internal standard at the beginning of sample preparation to correct for sample losses during sample processing and allow for absolute oligosaccharide quantification. Oligosaccharides were fluorescently labeled with 2‐aminobenzamide (2AB, Sigma) at 65 °C for 2 h. Labeled oligosaccharides were analyzed by high performance liquid chromatography (Dionex Ultimate 3000, Dionex, now Thermo Fisher) on an amide‐80 column (15 cm length, 2 mm inner diameter, 3 µm particle size; Tosoh Bioscience)] with a 50‐mmol L − 1 ammonium formate–acetonitrile buffer system. Separation was performed at 25°C and monitored with a fluorescence detector at 360 nm excitation and 425 nm emissions. HMO peaks were annotated based on standard retention times and quantified in reference to the internal standard.

The measured HMO structures included 2’-fucosyllactose (2′FL), 3-fucosyllactose (3FL), 3’-sialyllactose (3′SL), 6′-sialyllactose (6′SL), difucosyllactose (DFLac), LNT, LNnT, lacto-N-fucopentaose (LNFP) I, LNFP II, LNFP III, sialyl-lacto-N-tetraose b (LSTb), sialyl-lacto-N-tetraose c (LSTc), difucosyllacto-N-tetrose (DFLNT), disialyllacto-N-tetraose (DSLNT), lacto-N-hexaose (LNH), fucosyllacto-N-hexaose (FLNH), difucosyllacto-N-hexaose (DFLNH), fucodisialyllacto-N-hexaose (FDSLNH) and disialyllacto-N-hexaose (DSLNH).

Genome annotation and phylogenetic analysis

Quality of the reads was assessed using FastQC (version 0.11.5).^Citation33 Specifically, reads were filtered for length, while those with a quality score less than Q30 were discarded using Trimmomatic (version 0.39).^Citation34 De-novo contigs were generated using SPAdes v4.0.0 (parameters: –careful -k 21,33,55,77,99,127).^Citation35 Assembly statistics computed using QUAST v5.1.0rc1.^Citation36 Genome completeness and contamination were estimated with CheckM v1.1.2.^Citation37 Mean coverage was computed using CMseq v1.0.4 (https://github.com/SegataLab/cmseq). Genomes with completeness >95% and contamination <1.5% were included in the downstream analyses. Taxonomic assignment was performed using GTDB-Tk v.2.4.0^Citation38 against bacterial and archaeal genomes based on the Genome Database Taxonomy (GTDB; Release 09–RS220 (24th April 2024)). Unless specified, all the computational tools used in this work were applied with default parameters. Average Nucleotide Identity (ANI) was calculated using FastANI (with default settings (K-mer size = 16, threads count for parallel execution = 1, fragment length default = 3,000).^Citation39 We analyzed multiple genomes using a genome list with command options of “—ql” and “–al.” For visualization, we used the “−matrix” command. For subsequent analysis, genomes with over 80% similarity were selected using FastANI. The assemblies and plasmid were annotated using Prokka (v.1.14.5).^Citation40 A phylogenetic tree was constructed using all species from Bifidobacterium and Lactobacillaceae, following the methodology with GToTree (v1.8.2),^Citation41 utilizing the protein sequences of 138 curated single-copy genes specific to Actinomycetota (i.e., Bifidobacterium), and 119 curated single-copy genes for Bacillota (i.e., Lactobacillaceae). Homologous sequences were aligned and concatenated using MUSCLE v3.8.1551.^Citation42 The resulting phylogenetic tree was constructed with FastTree v2.1.10.^Citation43 The tree was rooted using Bifidobacterium animalis subsp. lactis (B. lactis) for Bifidobacterium and Lactobacillus paragasseri for Lactobacillaceae. The iTOL web tool version 6^Citation44 was employed for the tree visualization and annotation.

Genomics-based reconstruction of HMO metabolism

We used a subsystems-based approach implemented in the SEED platform^Citation45 for the reconstruction of HMO utilization pathways across the set of 39 Bifidobacterium and 14 Lactobacillaceae genomes. The genomes were annotated via RAST^Citation46 and uploaded to mcSEED (microbial community SEED), a clone of the SEED annotation environment. Homology and genomic context-based methods in mcSEED were used to analyze the representation of groups of orthologous genes implementing a curated set of 67 and 57 functional roles in bifidobacteria and lactobacilli, respectively, that contribute to HMO metabolism. These functional roles included components of glycan transporters, extracellular and cytoplasmic glycoside hydrolases (GHs), monosaccharide catabolic enzymes, and transcriptional regulators.^{Citation23,Citation47} The genomic distribution of functional roles was used to reconstruct glycan utilization pathways in two metabolic subsystems containing the genomes of bifidobacteria and lactobacilli isolates, respectively. For bifidobacteria, we reconstructed pathways for Lac, lacto-N-biose (LNB), LNT, LNnT, 2’FL, 3FL, DFLac, LNFP I, 3’SL, 6’SL, and lcHMO catabolism (Supplementary Figure S2 and Supplementary Table S2). The lcHMOs are HMOs with a degree of polymerization ≥5 (e.g., LNFP II/III, LNDFH I/II) for which separate utilization pathways could not be unequivocally predicted due to the insufficient information on the substrate specificity of corresponding HMO transporters. For lactobacilli, pathways for individual HMO (LNT and LNnT), and the products of HMO degradation (Lac, LacNAc, LNB, and lacto-N-triose [LNTriose], L-fucose [Fuc], N-acetylneuraminic acid [Neu5Ac], N-acetyl-D-glucosamine [GlcNAc], D-galactose [Gal], and D-glucose [Glc]) were reconstructed ( and Supplementary Table S3).

The reconstructed pathways contained alternative biochemical routes implemented by different subsets of catabolic enzymes and diverse glycan transporters corresponding to different ecological strategies of HMO metabolism (extracellular vs. intracellular). To enable comparison across different strains, the detailed pathway variants were simplified to binary utilization phenotypes (1 and 0) corresponding to the predicted ability to metabolize a specific HMO species or its degradation by-products.

Identification of antibiotic resistance genes, virulence factors and plasmids

The antimicrobial resistance genes (ARGs) were detected from the assembly output of Resistance Gene Identifier (RGI) v6.0.3 with the CARD database v3.3.0,^Citation48 and Strict significance based on CARD curated bitscore cutoffs. The virulence factors were identified by ABRicate v1.0.1 (https://github.com/tseemann/abricate.) using the Virulence Factor Database (VFDB)^Citation49 with a cutoff coverage of 80 and 70% identification. Plasmids were detected by the Platon v1.7,^Citation50 which is based on a taxon-independent approach to extract both plasmid and genomic sequences. The pathogenicity toward humans was predicted by PathogenFinder v1.1 (https://cge.food.dtu.dk/services/PathogenFinder/).

Growth curve analysis

Growth curves were analyzed using the R package gcplyr^Citation53 which carries out model-free and non-parametric analyses to directly quantify attributes of growth dynamics. The underlying OD₆₀₀ of the base medium was subtracted from each condition at each timepoint (blank subtraction). The area under the curve (AUC), the maximum growth rate (r) and carrying capacity (k) were calculated to study the impact of HMOs on the growth properties of the strains.

Log fold changes in growth curve metrics were calculated by dividing the values of the growth curve metric from condition 1 by condition 2 and applying logarithmic transformation. Growth curve metrics between the different conditions of each strain were compared by Mann-Whitney U tests with p-values adjusted using the BH method.

HMO utilization analysis

HMO utilization at 24–48 h was calculated in reference to the initial HMO concentrations in the medium at the beginning of each experiment (t = 0) and represented as the percentage of HMOs remaining in the medium. The ComplexHeatmap R package^Citation54 was used to draw the HMO utilization profile of each group of species.

Resistance and HMO consumption profiles association analyses

Bray-Curtis distance-based Multivariate Redundancy Analysis (db-RDA) was conducted to test the correlation between HMO consumption capacity and antibiotic resistance profiles, by terms of MIC values. Both datasets where log transformed and the mean of overall HMO consumption and MIC values of each species were used for db-RDA. Antibiotic susceptibility was indicated as a MIC value of 0. The ‘adonis2’ function from vegan R package^Citation55 was used to perform permutations to evaluate overall differences in HMO consumption between species and to assess the association with the resistance profile (taken as the variables). Mann-Whitney U tests with p-values corrected using the BH method were performed to study the differences between the log transformed mean HMO consumption values and the log transformed mean MIC values of each species. The significance value was determined based on 999 permutations.

Genotype-to-phenotype predictions of HMO utilization

We analyzed the representation of reconstructed HMO utilization pathways encoded in the genomes of Bifidobacterium strains and associated them with their HMO utilization profiles, expressed as the percent of remaining HMOs in the culture supernatant after 48 h of growth ( and Supplementary Figure S2). A detailed description of genotype–phenotype comparisons can be found in Supplementary Note S2. Overall, the in silico predictions were in good agreement with in vitro HMO utilization profiles for most strains; however, several discrepancies were observed.

B. bifidum strains degraded most measured HMO structures, consistent with the presence of multiple HMO-acting extracellular GHs encoded in the genomes of these strains (Supplementary Table S2). However, some strain-level variability of HMO utilization was noted. For example, B. bifidum IATA039 could not degrade 2’FL, likely due to a 5’-end truncation of the gene encoding the alpha-1,2-fucosidase (BbAfcA; GH95) required for 2’FL metabolism^Citation56 (Supplementary Figure S3). We also observed an interplay between concentrations of 3FL and DFLac in the supernatants of all strains except B. bifidum IATA049 and IATA039 (). The increased 3FL concentration could be explained by the unequal (separated in time) removal of fucose residues from DFLac by the alpha-fucosidases of B. bifidum. The concentration of LNnT in the supernatants of B. bifidum strains varied considerably, decreasing in some strains (IATA049, IATA102, IATA139, IATA148) and increasing in others (IATA001, IATA016). Since LNnT is an intermediate product of the degradation of multiple more complex HMOs (e.g., LNFP III), the observed results may reflect strain-level variability in the kinetics of long-chain HMO degradation. B. bifidum strains IATA102, IATA005, IATA001, and IATA049 were able to grow in the medium supplemented with LNnT as a sole carbon source, confirming their capacity to utilize this HMO (Supplementary Figure S4a).

B. infantis strains efficiently utilized multiple HMO species, including 2’FL, 3FL, LNFP I, LNT, LNFP II, and DFLNT, consistent with the in silico predictions (). B. infantis IATA105, however, did not utilize LNnT, potentially due to mutations in genes encoding the respective utilization pathway. B. longum strains displayed varying capacities to utilize HMOs. B. longum strains IATA034, IATA062, and IATA075 degraded LNT, LNFP I, LNH, and FLNH via a mechanism involving an extracellular lacto-N-biosidase (LnbX; GH136)^Citation57 (Supplementary Figure S2). Other B. longum strains utilized only LNT, likely via an intracellular mechanism, although the variant of the LNT transporter (GltABC) found in most of these strains was previously shown to have low affinity to LNT.^Citation58

All B. pseudocatenulatum strains were predicted to utilize LNT, which was confirmed by their ability to consume over 95% of the LNT present in the pHMO mixture (). In agreement with in silico predictions, B. breve IATA027 and IATA153 efficiently consumed both LNT (>90%) and LNnT (>75%) (). The remaining B. breve strains depleted only LNT. Consistent with glycoprofiling data, B. breve IATA048 and IATA131 grew in a medium supplemented with LNT but not LNnT (). The discrepancy between genomic predictions and experimental observations of LNnT utilization could be attributed to mutations in the genes involved in the metabolism of this oligosaccharide. For example, in B. breve IATA131, the nahS gene was disrupted by a premature stop codon, suggesting that the encoded substrate-binding component of the LNnT-specific ABC transporter might be nonfunctional (Supplementary Figure S3).

B. adolescentis and B. lactis strains were predicted to lack the capability to utilize any HMOs due to the absence of required transport systems and GHs ( and Supplementary Table S2). Glycoprofiling results supported this prediction for most strains. Unexpectedly, B. lactis strains IATA010, IATA020, and IATA029 depleted LNT from the pHMO mixture and grew in a medium where LNT served as the sole carbon source (Supplementary Figure S4a).

We also analyzed the representation of HMO utilization pathways in 14 Lactobacillaceae and compared it with HMO glycoprofiling data ( and Supplementary Table S3). Unlike the bifidobacteria tested, most lactobacilli except L. mucosae IATA082, L. rhamnosus IATA115, and L. sakei IATA088 were predicted not to utilize any HMO structures. Despite lacking complete HMO utilization pathways, all analyzed strains were predicted to metabolize mono- and disaccharide products of HMO degradation, namely, lactose, N-acetylglucosamine, galactose, and glucose, and some strains were predicted to also utilize N-acetyllactosamine, lacto-N-biose, and lacto-N-triose II ().

L. mucosae IATA082, L. rhamnosus IATA115, and L. sakei IATA088 were predicted to partially degrade LNnT via a mechanism involving an extracellular beta-galactosidase (Supplementary Table S3). Consistent with this prediction, L. mucosae IATA082 and L. sakei IATA088 depleted LNnT from the pHMOs mixture and grew in the medium with LNnT as the sole carbon source ( and Supplementary Figure S4b). Notably, these were the only strains lacking the N-acetyllactosamine utilization pathway and failed to grow when N-acetyllactosamine was the only carbon source in the medium (Supplementary Figure S4b). L. rhamnosus IATA115, however, did not degrade LNnT present in the pHMOs mixture. Other discrepancies between genomic predictions and experimental data include the unexpected partial degradation of LNFP III by L. paracasei strains IATA083 and IATA109 (). Additionally, L. mucosae IATA082 depleted LNT, LNH as well as fucosylated HMOs (LNFP III, FLNH, DFLNH) despite lacking genes encoding alpha-fucosidases ( and Supplementary Table S3).