This study used a subset of samples from the Child Health and Residence Microbes (CHaRM) study (30, 31) conducted in Brisbane, Australia, selected based on the depth and quality of available DNA sequencing data. Briefly, the CHaRM study was conducted as an adjunct to the GUMLi (Growing Up Milk “Lite”) trial, a multi-center double-blind, randomized control trial (Brisbane, Australia, and Auckland, New Zealand) investigating the effect of a fortified young child formula (Active) compared to unfortified cow's milk (Control) on various outcomes in toddlers between the ages of 1 and 2 years (30).

The Active group received fortified milk supplemented with synbiotics: a formulation with 7.8 × 10⁸ cfu/100 ml of B. breve M-16V, 1.8 g of short-chain galactooligosaccharides (scGOS), and 0.2 g of long-chain fructooligosaccharides (lcFOS) per 100 ml (Supplementary Table 1). The participants consumed two 150 mL servings of formula daily. An energy-matched, non-fortified cow’s milk was used as a Control. In addition, the Active milk was fortified with 1.3 mg of Fe and 1.2 µg of vitamin D and vitamin C to increase the bioavailability of non-heme Fe (ferrous citrate). Further compositional differences between the Active and the Control milk are described in the study by Wall et al. (30) and the corresponding ethical approval can be found in the study by Matsuyama et al. (29).

Of the 48 children who completed the CHaRM study, a subset of 29 (18 subjects from the Active group and 11 from the Control group, details in Supplementary Table 2), for which a minimum of 15,000 high-quality sequencing reads were generated, was included in the current study.

Blood samples were collected at baseline and month 12 as described by Lovell et al. (28). In this study, we aimed to link the same parameters [serum iron, ferritin, and hemoglobin (Hb)] with the functional analyses of the gut microbiome.

All fecal samples were collected and preserved, and gut microbiota DNA was extracted as described previously (31). Briefly, a repeated beat-beating and a column (RBB + C) technique was used based on Yu and Morrison (32) and modified to suit the automated Maxwell 16MDx system (Promega, Madison, WI, USA). An aliquot of each fecal sample was mixed with beads, homogenized with lysis buffer, and homogenized in a Precellys 24 homogenizer (Bertin Corp, Rockville, MD, USA). After the protocol was completed, the supernatant was transferred into Maxwell 16 MDx cartridges for elution. At this step, a non-template control was used to check the quality of each batch of buffers (lysis and elution). DNA samples were normalized to sample concentration before the creation of 16S rRNA gene amplicon libraries following Illumina library preparation instructions, using the primers of Klindworth et al. (33) to amplify the V3-V4 hypervariable regions of the bacterial 16s rRNA gene region. The libraries were then sequenced for 300 paired-end cycles using the Illumina MiSeq platform.

In parallel, we used four targeted quantitative polymerase chain reaction (q-PCR) assays to detect total bacteria, total Bifidobacterium, B. breve group, and B. breve M-16V strain. Briefly, the normalized DNA samples were used for all the amplified quantitative real-time PCR amplification performed using an ABI Prism 7900HT (Applied Biosystems, California, USA) using the TaqMan Universal Master Mix (Applied Biosystems, Austin, TX, USA). Standards for each target were generated from their genomic DNA and were amplified with specific primers, and purified using a MinElute PCR Purification Kit (Qiagen, Valencia, CA, USA) as per the manufacturer’s protocol. Final standard concentrations were quantified using a NanoDrop 2000 (Thermo Scientific, Waltham, MA, USA). For a detailed description of the methods, including the oligonucleotide sequences of the primers and probes used for the q-PCR analyses, please see Chua et al. (34).

Samples were prepared using a Microlab NIMBUS (Hamilton Robotics, Reno, Nevada, USA) in a customized setup. The q-PCR conditions were 1 cycle of 95°C for 20 s, followed by amplification at 95°C for 1 s, and 62°C for 20s for 40 cycles, and then 1 cycle of 95°C for 15 s, followed by 60°C for 15 s, and 95°C for 15 s. SDS 2.4 (Applied Biosystems) was used to visualize and check the abnormality of the curves deviating from standard amplification. Raw data were then exported into Microsoft Excel, where the Ct values were transformed to log copy numbers per gram of feces for the statistical analyses.

The data analysis performed in this sub-study differs from that of Matsuyama et al. (29) by using a higher sequence coverage cut-off and optimized sequence data quality control. Briefly, pre-processing and filtering of the raw 16S rRNA gene sequence data resulted in an average of 18,916 high-quality reads per sample, and each rarefaction analysis showed slopes reaching a plateau, indicating that a sufficient sequencing depth had been achieved. An adaptation of the “Quantitative Insights into Microbial Ecology” (QIIME) v1.9.1 package (35) was used for the processing and analysis as described in (29). The paired-end reads were demultiplexed and trimmed (q > 20) before being merged using QIIME. These were dereplicated and counted using mothur (36) and reads with a low abundance (less than 2 reads across all samples) were discarded. Chimeras were removed using VSEARCH (37), using the RDP gold database (38) as a reference. Reads that contained PhiX or adapters as defined in Deblur (part of QIIME2) (39) were eliminated. Taxonomic assignment was performed using the RDP classifier (40) against the SILVA_138.1 database (41), from which results where the sequences were aggregated at the genus and phylum level, were further explored. Reads with eukaryotic assignments and reads with a low relative abundance of up to 0.0005% across all samples were excluded from further downstream analysis. Samples were rarefied, and α-diversity was calculated using the phyloseq (42) and vegan (43) packages in R software (version 3.5.1) (44).

In addition to the amplicon-based analysis, we used a metagenome functional content bioinformatics package that provides potential functional category abundances for a microbial community from marker genes extracted from full genomes using 16S operational taxonomic unit (OTU) profiles as implemented in the PICRUSt software package (45) using the default parameter. Briefly, the 16S OTUs were normalized by copy number and mapped onto a phylogenetic tree created from an ancestral state reconstruction (ASR) reference genome database with more than 10× coverage. The predicted metagenome functional abundance counts per sample were then collapsed into functional categories based on KEGG Gene Orthology pathways for analyses. We investigated the change in predicted functionality abundances using the delta difference in the counts at baseline and at month 12.

PRISM GraphPad 8 (version 8.4.3) was used for all the statistical analyses. The change from baseline within each group was analyzed using a non-parametric pair-wise Mann–Whitney test. The baseline and end-of-intervention parameters of the groups were compared using a non-parametric Wilcoxon test. The effect size was calculated using Cohen's d, where the mean difference between the two groups was divided by the pooled standard deviation. The interaction of the intestinal microbiota composition with the intervention was visualized spatially using distance matrices in the principal component analysis (PCA), an unconstrained ordination method, using CANOCO 5 (2012) software. The Benjamini–Hochberg procedure (46) was used to reduce the false discovery rate (FDR) for all comparisons (within each group, at baseline, and at the end of the intervention) made using the predicted microbial functionalities.

The same findings in the main GUMLi study were observed in this subset of samples from the CHaRM study. There was no statistically significant difference in serum iron at baseline between the groups (p = 0.111), and significant differences were detected at the end of the study within the Active group (p < 0.0008, effect size d = 1.390574) but not within the Control group (p = 0.898) (Figure 1). Similarly, there was no difference at baseline between the groups for transferrin, 25(OH)D and ferritin, (p > 0.05) amongst this subset of samples (data not shown). We did not observe a significant difference between the groups in hemoglobin levels at baseline (p = 0.346) or at month 12 (p = 0.580).

The hemoglobin levels were significantly increased from baseline to month 12 in the Active group (p < 0.0001, effect size = 1.125876), while it was significant in the Control group for the same comparison (p = 0.024). The Vitamin D indicator, 25(OH)D, was significantly reduced from baseline to month 12 in both the Active (p = 0.003, effect size = 0.826) and the Control groups (p = 0.001, effect size d = 1.399) although they were still within the healthy range according to the Endocrine Society’s Clinical Practice Guidelines (>50 nmol/L) (47).

At the end of the study, there were significantly higher transferrin levels in the Active group from baseline to month 12 (p = 0.0029) compared to the Control group from baseline to month 12 (p = 0.527). In addition, serum ferritin concentration was significantly lower at month 12 compared to the baseline in the Control group (p = 0.001). A significant difference was only observed at month 12 between the Active and Control groups (p = 0.0002) (Figure 1).

We quantified the total bacteria [copy number in feces (log₁₀/g)] in the fecal samples and observed an increase that was significantly different from baseline to month 12 in the Active group [mean difference 0.84 (log₁₀/g), p = 0.001] and in the Control group [mean difference 0.30 (log₁₀/g), p = 0.032]. No significant difference was observed at baseline or at month 12 between the intervention groups. We observed no significant difference in the total count of Bifidobacterium copy number in feces (log₁₀/g) at baseline or at month 12 between the Active and Control groups. There was a significant decrease in the B. breve group in the Control group from baseline to month 12 [0.79 (log₁₀/g), p = 0.049]. This decrease was not observed in the Active group. The supplemented probiotic strain, B. breve M-16V, was only detected in the Active group and we observed a significant increase [1.66 (log₁₀/g), p ≤ 0.001] between baseline and month 12 (Supplementary Figure 1).

The major bacterial phyla present in this subset of samples were the same as previously reported (29), namely, Firmicutes, Bacteroidetes, and Actinobacteria, characteristic of healthy human gut microbiota. Increased α-diversity among this subset of samples was also demonstrated (Supplementary Figure 2). Here, we report the most abundant bacteria genera showing significant differences between baseline and month 12. Bifidobacterium decreased in the Control group from 20.6% at baseline to 7.9% at month 12 (p = 0.020). The Escherichia-Shigella group significantly decreased during the intervention in the Active group (from 4.17% to 0.18%) (p = 0.0017) and the Control group (from 1.8% to 0.30%) (p = 0.004) (Supplementary Figures 3A,C). For the other genera investigated, such as Bacteroides, Collinsella, and Faecalibacterium, no significant differences were observed (Supplementary Figure 3). Furthermore, the genus Veillonella showed a significant decrease as observed previously (29), most noticeably in the Active group.

In Figure 2, the 10 most abundant predicted functionalities were subjected to PCA. PCA Axis 1 explained 41.13% of the variation, with the most significant predicted functions pointing (red arrows) to the right quadrants. However, the majority of the samples, regardless of group or timepoint, congregated and overlapped in the center of the graph. Functional abundance predicted by PICRUSt was compared between the intervention groups at either baseline, month 12 or within each intervention group from baseline to month 12. No significant differences were observed at baseline or month 12 after applying multiple correction (FDR). In the Active (baseline to month 12) and Control group comparison, several functions remained significant after FDR correction (refer to Supplementary Table 3). Among the listed functions, enterobactin biosynthesis was reduced in both the Active (p = 0.00182) and Control (p = 0.047) groups at the end of the study.

For participants in the Active group at month 12, we observed a trend toward the right quadrant of the PCA plot that was in line with the explained genera and blood parameters (iron, transferrin, and Hb). This separation explained 25.41% of the observed variation in the microbiomes during the intervention (in Supplementary Figure 4). In the Control group, every blood parameter except hemoglobin trended in the opposite direction, although the bacteria genera were also shifted toward the participants at month 12, explaining 46.95% of the variation observed. In addition, Bifidobacterium pointed in the opposite direction (toward baseline), while the remaining top 10 genera pointed toward month 12, indicating a possible normal microbiota maturation with age. We evaluated the influence of the intervention on the microbiota over time by using the delta change in the relative abundance from baseline to month 12. The PCA analysis of the 10 most abundant bacteria genera explained 31.29% of the observed variation, as shown in Figure 3. The figure shows that the clustering in the Active group was mainly driven by Bifidobacterium and Bacteroides. The clustering in the Control group was mainly driven by Prevotella 9 and Sutterella. A combined unconstrained PCA plot is presented in Figure 4. PCA axis 1 explained 25.59% of the variation seen in the microbiota, with the baseline samples clustered in the bottom left quadrant of the plot. The Escherichia-Shigella group pointed in the opposite direction to the other influencing genera that clustered at month 12 (for both Active and Control group samples). On the same plot, serum ferritin and 250HD pointed toward the left quadrants (higher levels at baseline) while serum iron, transferrin, and hemoglobin pointed toward the month 12 samples (increased levels at the end of the intervention). This is consistent with the results shown in Figure 1.