Demographic data on study participants are shown in Table 1. Two sets of subjects were randomly chosen for serum metabolomics analysis at 4 months and microbiome analysis at three timepoints: baseline, 4 months, and 12 months (with some subjects overlapping in the two sets). Approximately equal numbers of male and female subjects comprised the breast-fed reference cohort (BF) and the cohort provided the Lactobacillus paracasei spp. paracasei strain F-19 probiotic (F19); however, there were more female than male subjects in the standard formula group (SF) in our dataset. Additionally, a greater number of infants were delivered vaginally in the BF group compared to the F19 or SF groups, and the F19 group had a greater number of vaginally born infants compared to the SF group. Mean birthweight was not significantly different amongst the groups. We previously reported that both formula groups experienced similar infectious episodes during the intervention, but only the SF group had significantly more days and episodes of fever than the BF group (29). For these subsets of infants, we did not find any difference in the frequency of antibiotic use or diarrheal episodes among groups. We observed significantly higher serum IFN-γ in the SF group compared to the BF and F19 groups, which is what was published previously for this cohort (30). For the subset of infants used for metabolomics, we also observed higher serum IL-2 (general T cell stimulation) in the F19 group compared to the BF group, but this was not significant for the subset of infants used for microbiota analysis (data not shown).

Based on the weighted UniFrac distance matrix, principal coordinates analysis (PCoA) revealed that at baseline (3 wk ± 7 d), the infants assigned to the F19 group were significantly different from the BF group (pairwise PERMANOVA q = 0.033; n = 999 permutations) but were not different from the SF group (q = 0.174). At the end of the intervention, when infants were 4 months of age, the two formula groups were not different and together were significantly different from the BF group (pairwise PERMANOVA q = 0.015 (BF-SF), 0.015 (BF-F19); n = 999 permutations). By 12 months of age, fecal microbiota was indistinguishable between all three groups (Figure 1A).

Differential abundance analysis revealed that at 4 months, the feces of infants in the F19 group were enriched in members of the Lactobacillus genus compared to the SF group (ANCOM FDR < 0.05, tested on relative abundance), and were not significantly different from the BF group (Figure 1B). Aside from higher relative abundance of Lactobacillus, no differences at the genus level were observed in the F19 group relative to the SF group. Interestingly, the relative abundance of Bifidobacterium was not significantly different between the groups (Figure 1B). Aside from higher Lactobacillaceae in the F19 group compared to the SF group at 4 months, the relative abundances of bacterial families were comparable between the formula groups. Both formula groups were different from the BF reference group with respect to the relative abundance of several bacterial families including Lachnospiraceae, Veillonellaceae, Enterobacteriaceae, Erysipelotrichaceae, and Enterococcaceae (Figure 1C). Notably, group differences were no longer apparent at 12 months, as the relative abundance of bacterial families including Lactobacillaceae were similar in all three groups (ANCOM FDR>0.05). This showed that increased Lactobacillus through F19 supplementation did not persist 8 months after the end of intervention. Throughout the study, Bifidobacterium was a dominant genus in all groups. Analysis of the diversity of the microbiome revealed that BF infants had significantly lower Shannon diversity at all timepoints compared to the SF group, with lower evenness at all three timepoints, and lower growth of richness from baseline to 4 months (Figure 1D). The F19 group tended to have lower Shannon diversity than the SF group at 4 and 12 months, which was reflected as lower evenness at both timepoints, and was not significantly different from the BF group at either time point.

We have previously shown that in both the post-prandial and semi-fasting states, there are profound differences in the serum metabolome between BF and FF infants (4–7, 9, 11). In this study, as expected, the SF group differed from the BF reference group (Figures 2A,B), with metabolome differences similar to what was previously reported. Here, significantly higher levels of several amino acids (including valine, isoleucine, leucine, lysine, threonine, phenylalanine, asparagine, tyrosine, alanine, and histidine) and other metabolites [3-hydroxyisobutyrate, 2-hydroxybutyrate, 3-hydroxyisovalerate, creatine, carnitine, dimethyl sulfone (DMSO₂), galactonate, and urea], and lower levels of 1,2-propanediol, ketone bodies (3-hydroxybutyrate and acetone), glutamine, methanol, formate, fumarate, citrate, succinate, pyroglutamate, serine, dimethylamine, creatinine, myo-inositol, 2-oxoglutarate, and N,N-dimethylglycine were observed in the serum of SF infants compared to BF infants (p < 0.05). Comparison of F19 infants to BF infants revealed a similar pattern in metabolites but with fewer significantly different metabolites and lower magnitudes of effect sizes than comparison of SF infants to BF infants (Figure 2C). Metabolites higher in F19 infants compared to BF infants included 3-hydroxyisobutyrate, DMSO₂, creatine, urea, 2-hydroxybutyrate, and several amino acids (including valine, isoleucine, threonine, lysine, and leucine). Compared to BF infants, F19 infants exhibited lower levels of glutamine, 1,2-propanediol, 3-hydroxybutyrate, myo-inositol, serine, creatinine, formate, dimethylamine, N,N-dimethylglycine, betaine, proline, methanol, ornithine, and taurine.

Given that the F19 infants had a serum metabolome more similar to SF infants than to BF infants, we sought to determine the difference between the F19 and SF groups. Using a linear mixed-effect model, with hospital as a random effect, p-values and Cliff delta effect sizes were calculated and are summarized in a volcano plot (Figure 3A). Compared to the SF group, the F19 group had higher levels citric acid cycle intermediates (citrate, succinate, and fumarate) and methanol, and lower levels of several amino acids (tyrosine, lysine, phenylalanine, histidine, proline, leucine, and threonine), 3-hydroxyisovalerate, ornithine, betaine, urea, taurine, and mannose. The metabolites with at least a medium effect size and a p < 0.05 were summarized in Figure 3B. One of the branched chain amino acids, leucine, as well as taurine and urea were lower in F19 infants compared to SF infants with medium effect sizes, but these were not significant (p < 0.1).

Several SF (n = 6) and F19 (n = 12) infants had serum metabolomes consistent with a BF infant's metabolome, with substantial levels of 1,2-propanediol [a metabolite that is formed by fermentation of milk oligosaccharides by select strains of bifidobacteria that contain the genes to catabolize fucose (31–33)] and ketone bodies, with significantly lower levels of branched chain amino acids. Further, there were a few infants in the BF group that had levels of these metabolites that were consistent with FF infants (4–7, 9, 11). Removal of these samples from analysis resulted in no difference in our reported results suggesting that their presence in the dataset did not inadvertently skew the results.

This clinical trial was a multicenter double-blind, randomized controlled trial with infants from 21 ± 7 days until the end of the 4^th month of age. Infants were followed until the age of 1 year. The design and clinical results of the study were previously published (29). The trial was approved by the Institutional Review Board (IRB) at the University of California Davis and the regional Ethical Review Boards in Nanjing, Shanghai, and Beijing in China (ClinicalTrials.gov Identifier: NCT01755481). Inclusion criteria encompassed healthy infants born full term with a birthweight between 2,500 g and 4,000 g. For infants in the breast-fed reference group (BF), inclusion criteria also included exclusive breastfeeding from birth, with mothers intending to breast-feed for at least 5 months and providing at least 80% of the calories to their infants from breast milk. Additional inclusion criteria for the two formula groups (FF) included infants of mothers who either could not breast-feed, or voluntarily resigned from breast-feeding by 28 days.

Except for the addition of the probiotic Lactobacillus paracasei ssp. paracasei strain F19 to one of the formulas [at a dose of 10⁸ CFU/L; 8.33 × 10⁷ CFU/day based on an average formula consumption of 833 mL/day (29)] (F19), the composition of the study formulas was exactly the same. The formulas were manufactured in Hohot, China using a bovine milk powder provided by Arla Foods amba, Denmark, and the L. paracasei ssp. paracasei strain F19 provided by Chr. Hansen, Denmark. The nutrient composition of the formula was previously published (29).

Prior to the intervention, infants in the FF groups were provided standard formula (SF) if formula feeding had been initiated. Between the 5^th and 6^th month of life, infants in both FF groups were provided SF. For BF infants, if milk supply was insufficient, SF was used, but was not to exceed 20% of total calorie intake based on a monthly three-day intake record (29). Although feeding of other foods or formulas was not recommended during the intervention, data on consumption of other foods or formulas was not recorded. Parents were recommended to introduce complementary foods to the infants no later than 6.5 months of age. The group code was blinded from study staff and enrolled participants, and was not broken until sample extraction, lab analysis, and generation of concentration/relative abundance data had been completed.

To confirm consumption of L. paracasei ssp. paracasei strain F19 in the F19 group, the amount in the stool was determined using quantitative PCR using the same primers and conditions as described by Sieuwerts and Håkansson (54).

A total of 150 subjects (n = 50 per group) were randomly chosen from 179 BF, 167 SF and 167 F19 infants completing the study (29), and serum samples collected from those infants at 4 months of age were used for metabolomics analysis. Out of the potential 150 samples, 6 samples had too small (< 50 μl) or no volume, 15 samples were hemolyzed (identified by a bright red color), 3 samples from the F19 group were from infants with no detectable levels of L. paracasei ssp. paracasei strain F19 in their stool, and 3 samples had noisy NMR spectra or potential contamination. Exclusion of these samples left 123 samples (n = 42 BF, n = 40 SF, and n = 41 F19) for metabolomics analysis.

Serum samples were prepared by filtering through 3,000 MW cutoff Amicon filters (Merck Millipore, MA, USA) followed by the addition of potassium phosphate buffer (pH 6.7) and an internal standard containing 5 mM 3-(trimethylsilyl)-1-propanesulfonic acid–d6 (DSS-d6) and 0.2% NaN₃ in 99.8% D₂O as previously described (6). NMR spectra were acquired on a Bruker Avance 600 MHz NMR spectrometer at 25°C using the Bruker noesypr1d experiment as previously described (7). NMR spectra were processed using Chenomx NMR Suite v8.2 Processor (Chenomx Inc., Edmonton, Canada) (RRID:SCR_014682), and metabolites were identified and quantified using Chenomx NMR Suite Profiler v8.2 based on the internal standard DSS-d6 as described (55). Quantified metabolite concentrations were corrected for dilution and are expressed as absolute concentrations in μmol/L.

One hundred and fifty-one subjects (n = 51 BF, n = 50 SF and n = 50 F19) from the original cohort were randomly chosen for 16s rRNA gene sequencing, and a total of 453 fecal samples representing three time points [baseline (21±7 d), end of the 4^th month and 12 months] from 151 infants were used. Fecal samples were collected by a parent or guardian from the infant's diaper, and placed in a provided container, which was placed in a freezer bag and stored at −20°C in a home freezer until the hospital visit (29). At the hospital, samples were placed in a −20°C freezer and stored until transport on dry ice to the Netherlands Organization for Applied Scientific Research (TNO) for analysis.

DNA from fecal samples was isolated and V4 hypervariable region of the bacterial 16S rRNA gene was amplified with universal bacterial 515F−806R primer pair as described (6). Prepared libraries were sequenced on the paired-end 2 × 200 bases Illumina MiSeq platform with a MiSeq Reagent Kit v2 (MS-102-2003, Illumina).

Paired sequence reads from 453 selected samples (BF n = 51, SF n = 50, F19 n = 50) and 48 quality control samples were pre-processed in QIIME2 (version 2019.4; https://qiime2.org/) as described (6). The following subject's samples were excluded from sequencing analyses: (1) one subject in the F19 group identified as quantitative-PCR negative for the F19 strain at 4 months; (2) 13 non-exclusive breastfed infants in the BF group; and (3) infants treated with antibiotics during the intervention (3 BF, 9 SF, and 6 F19), leaving a total of 37 BF, 43 F19, and 41 SF. Raw sequence and processed files have been deposited through QIITA (study ID: 12874) in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB38295 (https://www.ebi.ac.uk/ena/browser/view/PRJEB38295).

To determine associations in the microbiota with diet, principal coordinate analysis (PCoA) was performed. Comparisons for beta diversity metric estimates were made by computing pairwise Permutational Multivariate Analysis of Variance (PERMANOVA) (56) based on 999 permutations. Kruskal-Wallis tests were used for comparing alpha diversity estimates between categorical variables. Longitudinal analyses were performed on changes in alpha diversity estimates over time using paired-difference testing with the q2-longitudinal plugin in QIIME2, followed by visualization using qiime2R (ver. 0.99.1) and ggplot2 (ver. 3.1.1) packages in R. Taxonomic assignments of representative sequences were conducted for 515F/806R primer pair through a naïve-Bayes classifier trained on the Greengenes ver. 13_8 database at 99% OTU similarities. Differential abundance analysis of each taxon was determined using ANCOM (Analysis of composition of microbiomes) (57). Taxa bar plots were generated using vegan (ver. 2.5-5) and ggplot2 packages in R.

For metabolomics, statistical computation was performed in R (ver. 3.5.3.) and visualized using the ggplot2 package unless stated otherwise. Principal component analysis (PCA) (prcomp function; mean-centered, non-scaled, and log-transformed data) was used to visualize the data. Groups were compared using pairwise PERMANOVA on Euclidean distances. Two-way Analysis of Variance (ANOVA, aov function; log-transformed data) were computed to test for significant effects of diet, sex, region, hospital, and mode of delivery as previously described (6). The magnitude of the diet effect on the metabolome was estimated by calculating Cliff's Delta effect size estimates (cliff.delta function).