Effect of Bovine Milk Fat Globule Membrane and Lactoferrin in Infant Formula on Gut Microbiome and Metabolome at 4 Months of Age

Abstract Background Milk fat globule membrane (MFGM) and lactoferrin (LF) are human-milk bioactive components demonstrated to support gastrointestinal and immune development. Significantly fewer diarrhea and respiratory-associated adverse events through 18 mo of age were previously reported in healthy term infants fed a cow-milk–based infant formula with an added source of bovine MFGM and bovine LF through 12 mo of age. Objectives The aim was to compare microbiota and metabolite profiles in a subset of study participants. Methods Stool samples were collected at baseline (10–14 d of age) and day 120. Bacterial community profiling was performed via 16S rRNA gene sequencing and alpha and beta diversity were analyzed (QIIME 2). Differentially abundant taxa were determined using linear discriminant analysis effect size (LefSE) and visualized (Metacoder). Untargeted stool metabolites were analyzed (HPLC/MS) and expressed as the fold-change between group means (control to MFGM+LF ratio). Results Alpha diversity increased significantly in both groups from baseline to 4 mo. Subtle group differences in beta diversity were demonstrated at 4 mo (Jaccard distance; R2 = 0.01, P = 0.042). Specifically, Bacteroides uniformis and Bacteroides plebeius were more abundant in the MFGM+LF group at 4 mo. Metabolite profile differences for MFGM+LF versus control included lower fecal medium-chain fatty acids, deoxycarnitine, and glycochenodeoxycholate, and some higher fecal carbohydrates and steroids (P < 0.05). After applying multiple test correction, the differences in stool metabolomics were not significant. Conclusions Addition of bovine MFGM and LF in infant formula was associated with subtle differences in stool microbiome and metabolome by 4 mo of age, including increased prevalence of Bacteroides species. Stool metabolite profiles may be consistent with altered microbial metabolism. This trial was registered at https://clinicaltrials.gov as NCT02274883.


Introduction
Maturation of the infant immune system and gastrointestinal (GI) tract, including dynamic development of the gut microbiome, establishes the foundation for long-term metabolic and immune health (1). Diet is central to infant growth, health, and intestinal development, and human milk (or its substitution) is recognized as the most significant factor shaping the microbiome (2). Ongoing investigation into hu-man milk has provided new information on functional components, such as the milk fat globule membrane (MFGM). MFGM envelopes milk-fat globules, facilitates their secretion from the mammary gland, and is rich in both nonglycosylated and glycosylated bioactive proteins and lipids (3) and includes components with antimicrobial effects including gangliosides (4), oligosaccharides (5), and the glycoproteins butyrophilin, lactadherin, and mucins (6). Previous research suggests that human MFGM components may affect bacterial attachment, colonization, clearance, and metabolism in the oral cavity and the gut (7).
Due to homology between human and bovine MFGM, the incorporation of bovine MFGM (bMFGM) in infant formula (IF) may better approximate the composition of complex human-milk lipids (8,9). Animal studies have demonstrated the impact of bMFGM on microbiome outcomes. In piglet models, animals who received formula with added bMFGM displayed significant group differences in ascending colon and stool microbiota (10), as well as higher stool Proteobacteria and reduced Firmicutes (11). In rat pups receiving formula with added bMFGM, richness and evenness of gut microbial species increased compared with a control formula with added vegetable fat, and the abundance of Firmicutes and Proteobacteria were more similar to a group receiving dam's milk (3). In the same study, Lactobacilli were most abundant in damreared pups, present in the bMFGM group, and not detected in the control group. In stress-induced juvenile rats receiving a diet with bMFGM, lactoferrin (LF), and prebiotics, improvements in sleep outcomes corresponded to increased growth of select beneficial taxa and protection against the stress-induced decrease in alpha diversity (12). In infants fed formula with added bMFGM, moderate effects on overall microbiome structure (oral and stool) were demonstrated (13) and the common middle-ear pathogen Moraxella catarrhalis was less prevalent in oral swabs at 4 mo of age (vs. infants receiving formula with no added bMFGM) (14).
In healthy term infants randomly assigned to a cow-milk-based IF (MFGM+LF) with added bMFGM (Lacprodan ® MFGM-10; Arla Foods Ingredients) and bovine lactoferrin (bLF; 0.6 g/L; Friesland-Campina DMV) compared with a similar control formula (no added bMFGM) through 12 mo, an accelerated neurodevelopmental profile by 12 mo of age and effects on language by 18 mo of age was reported (15). Added bMFGM and bLF in formula was also associated with significantly fewer episodes of respiratory adverse events and diarrhea through 18 mo of age, in line with a prior report of decreased diarrhea in infants receiving MFGM (16). Because the effects of early nutrition on the GI microbiome may potentially contribute to overall health and developmental outcomes, stool samples were collected in a subset of study participants at 4 mo of age, with the goal of comparing GI microbiome outcomes and metabolite profiles.

Participants and study design
In this double-blind, controlled, prospective study (trial registration: https://clinicaltrials.gov/ct2/show/NCT02274883), healthy term infants (born at 37-42 wk of gestation) were randomly assigned to receive the following: a cow-milk-based IF or a similar IF that had added bMFGM and LF from baseline (10-14 d of age) through 365 d of age. Mothers who had decided to exclusively provide IF were screened for study eligibility and eligible infants had received IF feeding for ≥3 d prior to randomization. As described previously, the primary outcome was the Bayley Scales of Infant Development, 3rd edition (Bayley-III), composite score at day 365 (15). Stool samples were collected in a subset of participants to evaluate stool microbiota at baseline and day 120 (MFGM+LF: n = 27; control: n = 35) and metabolite profiles at day 120 (MFGM+LF: n = 26; control: n = 33). Additional criteria for in-clusion in the subset included the following: vaginal birth, no antibiotic/corticosteroid treatment or prebiotic/probiotic supplements, and access to a home freezer. For each analysis, a stool sample of approximately 5 mL was collected at home and stored in the home freezer. Samples were de-identified at the time of collection. After delivery to the study site, all samples were maintained at −80 • C until processing and analysis. The research protocol and informed consent forms observing the Declaration of Helsinki (including October 1996 amendment) were approved by the Shanghai Nutrition Academy Medical Ethical Committee.
Alpha diversity (within-sample diversity) analyses were performed using the QIIME 2 q2-diversity plugin to measure richness (as observed ASVs). To quantify the dissimilarity between microbial samples by study time point and group, we examined multiple beta diversity metrics in QIIME 2: unweighted UniFrac (a phylogenetic, qualitative metric), weighted Unifrac (phylogenetic, abundance-weighted), Bray-Curtis dissimilarity (nonphylogenetic, abundance-weighted), and Jaccard distance (nonphylogenetic, qualitative) (23). Feature tables were evenly subsampled at 27,000 sequences per sample prior to alpha and beta diversity analyses. Two-factor ANOVA tests and Wilcoxon signedrank tests were performed using the q2-longitudinal plugin (24) to test whether alpha diversity estimates differed between study time point and group. Pairwise distances were used to construct principal coordinates analysis plots and perform 2-way permutational multivariate ANOVA (PERMANOVA) tests (25) [as implemented in the Adonis method in the vegan R package (26), wrapped via the q2-diversity plugin] to test whether beta diversity estimates partitioned by study time point and group.
Differential heatmaps of taxonomic compositions in each sample (by study time point and group) were visualized using Metacoder (27). Analysis of composition of microbiomes (ANCOM) (28), linear discriminant analysis effect size (LefSE) (29), and Wilcoxon rank-sum tests (implemented in the Metacoder package) were also used to test differential ASV abundance by study time point and group.

Stool metabolite analysis
A minimum of 100 mg of stool sample was needed to prepare and complete metabolite analysis. Briefly, samples were prepared using a methanol extraction to remove the protein fraction while allowing maximum recovery of small molecules and subsequently prepared for either LC/MS or GC/MS analysis. Metabolomic profiling was performed by Metabolon, Inc. (Research Triangle Park, NC), as previously described (30,31). Differential abundance of metabolites between groups was identified using Kruskal-Wallis H-tests as implemented in Scipy (32) followed by Benjamini-Hochberg correction for multiple testing as implemented in Statsmodels (33).

Correlations between bacteria and metabolites
The QIIME 2 plugin version of mmvec (https://github.com/biocore /mmvec) was used to predict associations between stool bacterial ASVs and stool metabolites using the microbe-metabolite vectors (mmvec) method (34). This method utilizes neural networks to predict microbe-metabolite log-conditional probabilities (i.e., the probability of  observing a metabolite at a certain frequency as a function of microbial composition). Log-conditional probabilities between each microbemetabolite pair were calculated using mmvec with default settings, ranked by score, and top scores were visualized with a heatmap.

Age and feeding effects on stool bacterial community diversity
Observed species (richness, representing the count of unique ASV observed in each sample) were significantly different between baseline and day 120 (Kruskal Wallis, P < 0.001; data not shown). No group differences in bacterial composition or species richness (ASVs) from base-line to day 120 were detected (paired Wilcoxon tests; Supplemental  Figures 1 and 2). Age exerted the strongest effect on bacterial community composition by all measures of beta diversity ( Table 1), although age explained only a small degree of variance in beta diversity (max R 2 = 0.06 with unweighted UniFrac). Using Jaccard distance (R 2 = 0.01, P = 0.042), low-abundance, phylogenetically related species differed between groups (Figure 1). No significant differences between study feeding groups at baseline or day 120 were detected by Bray-Curtis and Unifrac (weighted or unweighted) analysis.

Relative microbial abundance
Analysis of core microbes for all samples revealed a high prevalence of select bacterial taxa at the species level by heatmap visualization (Supplemental Figure 3). Both study time points and groups were  (Figure 2). Differential stool bacterial species abundance between study feeding groups was tested using ANCOM, LEfSe, and MetaCoder (Wilcoxon rank-sum tests). These 3 tests were used to evaluate results using tests with different degrees of sensitivity and precision; the most conservative test, ANCOM, detected no significant features at either time point. LefSE and Wilcoxon tests were subsequently compared to assess what features may differ between feeding groups at baseline and day 120. No significant group differences in species were detected by ANCOM (data not shown), whereas LEfSe has shown several significant differences (Figure 3). At baseline, Trabulsiella sp. was more abundant in the MFGM+LF group and Clostridium tertium was more abundant in the control group. At day 120, Bacteroides plebeius and Bacteroides uniformis were more abundant in the MFGM+LF group and Bifidobacterium adolescentis and Ruminococcus sp. were more abundant in the control group.

Effects of study feeding on stool metabolite abundance
Untargeted metabolite analysis was performed on day 120 samples. Of the 569 total biochemicals detected (418 named and 151 unnamed), 20 were differentially abundant between feeding groups (Kruskal-Wallis, P < 0.05; Figure 4). Metabolites at higher abundance in the MFGM+LF group included carbohydrates (fructose, 2deoxyibose), steroids (lanosterol, pregn steroid monosulfate), lipids (lignocerate, 1-myristoylglycerol), N-acetylhistamine, and biliverdin. Metabolites at higher abundance in the control group samples included the medium-chain fatty acids (caprylate, caprate, and laurate), the bile acid glycochenodeoxycholate, and deoxycarnitine. Otherwise, the fecal metabolic profiles were generally similar between feeding groups. However, after applying Benjamini-Hochberg multiple test correction, the differences were not significant.

Discussion
The aim of the present analysis was to examine the impact of bMFGM and bLF in IF on stool microbiome and metabolome. We previously demonstrated improvements in neurodevelopmental testing outcomes and reductions in respiratory adverse events and diarrhea in infants receiving added bMFGM and bLF in formula (compared with a control) through 12 mo of age (15). In a subset of these infants (n = 62), we utilized V4-16S rRNA gene sequencing to characterize fecal bacterial community development during the first 4 mo of life.
Bacterial microbiota diversity can be quantified as the amount of variation in the microbial community within a sample (alpha diversity) or between samples (beta diversity). In many circumstances, diversity is considered an indicator of community stability and resiliencethat is, higher diversity is associated with more stable community and lower diversity relates to less stable community. However, the microbiota of infants receiving human milk is more stable over time and is characterized by lower alpha diversity compared with infants receiving IF (35). Results in the current study demonstrated no significant effects of the diet on alpha diversity at the bacterial family level, suggesting that species richness increased with infant age but study feeding had less impact on the rate of change. Our results agree with a previous clinical study that reported a minimal effect on stool bacterial microbiota diversity in infants receiving bMFGM in IF (36). Timby et al. (37) showed that Moraxella catarrhalis was significantly lowered in oral microbiota in infants fed formula with bMFGM. M. catarrhalis is a common pathogen in otitis media, which was significantly less abundant in infants receiving added bMFGM in formula (14). Our analysis has shown that only low-abundance, phylogenetically related species differed between groups, explaining a very low degree of variation in beta diversity.

FIGURE 5
Microbe-metabolite vector analysis identifies co-occurrence probabilities between microbial ASVs and metabolites. Each row represents an individual bacterial ASV, labeled by its species ID, and each column represents an individual metabolite. The intersection of each row/column indicates their co-occurrence probability as an LCP. High LCPs (red) indicate strong co-occurrence probabilities between those features, and low LCPs (blue) indicate negative or null correlations. Margin vectors are colored according to bacterial family affiliation of taxa displayed along the y-axis and metabolite super family classification of metabolites displayed along the x-axis. Only the most significant correlations (features with at least 1 LCP >1.5) are shown. ASV, amplicon sequence variant; LCP, log conditional probability score.
In the current study, LefSE analysis demonstrated an increase in 2 species of Bacteroides in the MFGM+LF group at day 120: B. uniformis and B. plebeius. Previously, stool B. uniformis has been reported as being higher in breastfed infants (38) and the authors further speculated that B. uniformis uses human-milk oligosaccharides (HMOs) to thrive in the infant gut from the early postnatal stage. This is supported by a separate report of a glycan utilization system identified in infant gut-associated B. uniformis that specifically facilitates the metabolism of mother's milk (39). Evidence for the effects of B. uniformis includes the finding that the strain B. uniformis CECT 7771, present in human milk, ameliorates metabolic and immune dysfunctions associated with obesity in a mouse model (40,41). Since both bovine and human MFGM contain glycolipids and glycoproteins, these may provide glycan components utiliz-able by B. uniformis. Furthermore, it has been observed that the genus Bacteroides as a whole are one of the taxa most affected by delivery mode; lower levels were reported in infants born by cesarean delivery (42) or prematurely (43), whereas infants born vaginally have higher amounts of Bacteroides throughout the first year of life (42,(44)(45)(46)(47). This suggests that Bacteroides, among other taxa, could be considered as a marker of relative eubiosis in infancy.
Other microbiome changes have been associated with components of bMFGM. Infants receiving a formula with added bovine gangliosides, a component of bMFGM, displayed an increase in Bifidobacteria and reduction in fecal E. coli compared with infants receiving a standard formula (48). In the present study, differences between the feeding groups were also observed at baseline. Because infants in the MFGM+LF group began receiving added bMFGM, bLF, and a blend of polydextrose and galacto-oligosaccharides in formula within the first few weeks of life, it is possible that an influence of even a few days of feeding early in life had an impact on the bacterial microbiota results.
Metabolomics analysis in the current study detected no significant differences between MFGM+LF and control groups. However, several important trends that are likely related to microbial metabolism were observed. For example, in the current study, stool lactate increased by day 120 in infants receiving added bMFGM and bLF in IF. Lactate is produced by Bifidobacteria, Bacteroides, Enterococci, and Streptococci (49), which are also primary colonizers of the infant gut and often associated with higher stool lactate in breastfed infants (50).
In a previous study of infants receiving added bMFGM in IF, lower fecal concentrations of amino acids (e.g., isoleucine, glutamate, phenylalanine, tyrosine, valine, and glycine) and 2-hydroxyisovalerate (branched-chain amino acid degradation product), cadaverine (lysine breakdown product), and 4-hydroxyphenyllactate (tyrosine breakdown product) were observed (vs. infants receiving control formula) (36). The present study also demonstrated a modest impact on infant stool metabolome, with added bMFGM and bLF in formula associated with lower medium-chain fatty acids and deoxycarnitine, although another amino acid metabolite, N-acetylhistamine, was elevated relative to the control group. Following Benjamini-Hochberg multiple test correction, none of the stool metabolites were significant despite differential abundance in MFGM+LF versus control groups.
The correlation analysis has shown several co-occurrence clusters (Supplemental Figure 3). None of these metabolites are among the differentially abundant metabolites, suggesting that the co-occurrence patterns here likely reflect the "core microbiome" and "core metabolome" in these subjects, rather than processes impacted by study feeding. We also acknowledge that the low sample to observation ratio limits our statistical power, and the uncorrected P values may contain many false positives. Given these limitations, we used this analysis as an exploratory approach for assessing how added bMFGM and bLF in IF may impact the stool metabolome. In addition, the analysis would have benefited from stool samples collected from a larger number of participants as well as samples from infants delivered by cesarean section.
In summary, our study has shown modest effects of supplementing IF with bMFGM and bLF on stool microbiome and metabolites, including (subtle) group differences in beta diversity (Jaccard distance) at 4 mo of age and higher abundance of select Bacteroides species. While such changes may contribute to the beneficial GI outcomes previously shown for formula with added bMFGM and bLF, this modest alteration of the gut microbiome is unlikely to represent the sole or primary driving factor, and additional mechanisms remain to be explored.

Acknowledgments
The authors' responsibilities were as follows-MC: developed the research plan for microbiome and metabolome outcomes and drafted the initial manuscript; NB: performed bioinformatics and statistical analysis; CLH: participated in study design and performed statistical analyses; JLW: participated in the study design; FL: conducted research and provided oversite at local study sites; CLB and CR: conceptualized and designed the study; SSW: supervised data collection, provided study oversight, and has primary responsibility for final content; and all au-thors: participated in interpretation of data, drafting or revising of the manuscript, and read and approved the final manuscript.

Data Availability
Data described in the manuscript will be made available upon request.