The association between breastmilk oligosaccharides and faecal microbiota in healthy breastfed infants at two, six, and twelve weeks of age

Several factors affect gut microbiota development in early life, among which breastfeeding plays a key role. We followed 24 mother-infant pairs to investigate the associations between concentrations of selected human milk oligosaccharides (HMOs) in breastmilk, infant faeces, and the faecal microbiota composition in healthy, breastfed infants at two, six and 12 weeks of age. Lactation duration had a significant effect on breastmilk HMO content, which decreased with time, except for 3-fucosyllactose (3FL) and Lacto-N-fucopentaose III (LNFP III). We confirmed that microbiota composition was strongly influenced by infant age and was associated with mode of delivery and breastmilk LNFP III concentration at two weeks, with infant sex, delivery mode, and concentrations of 3′sialyllactose (3′SL) in milk at six weeks, and infant sex and Lacto-N-hexaose (LNH) in milk at 12 weeks of age. Correlations between levels of individual breastmilk HMOs and relative abundance of OTUs found in infant faeces, including the most predominant Bifidobacterium OTUs, were weak and varied with age. The faecal concentration of HMOs decreased with age and were strongly and negatively correlated with relative abundance of OTUs within genera Bifidobacterium, Parabacteroides, Escherichia-Shigella, Bacteroides, Actinomyces, Veillonella, Lachnospiraceae Incertae Sedis, and Erysipelotrichaceae Incertae Sedis, indicating the likely importance of these taxa for HMO metabolism in vivo.


Results
The infants included in this study were delivered vaginally (n = 22) or via C-Section (n = 2). For one of the infants, data on sex, weight and place of delivery were not reported. The remaining infants included 13 boys and ten girls, and their average body weight at birth was 3,536 g and ranged from 3,024 g to 4,140 g. Five infants were born at a clinic, seven at home, nine at a hospital, and for two infants birth started at home but was completed at a hospital. Six infants had a common cold and one had diarrhoea during the study. The actual ages in days ranged within The Illumina HiSeq sequencing resulted in a total number of 8,550,719 (range: 3,383-421,482 per sample, M = 118,760, SD = 86,261, SE = 10,166) sequencing reads that passed the quality filtering (Fig. S1). A total number of 411 Operational Taxonomic Units (OTUs) were identified of which 83 OTUs were found in more than 5%, or at least four samples (Table S2).
The development of infant faecal microbiota. Alpha diversity analyses showed no statistically significant differences in richness and diversity in faecal microbiota of infants at different time points when tested with Kruskal-Wallis and Wilcoxon tests (Fig. S2). The Principal Component Analysis (PCA) showed that composition of faecal microbiota of individual infants was diverse, but despite the high inter-individual variation, the observed changes in microbial profiles were directional and progressed with time towards bifidobacteria-dominated communities (Fig. 1). Redundancy analysis (RDA) was performed on the OTU level infants based on their faecal microbiota composition at two, six and 12 weeks of age and the twenty best fitting OTUs; (b) Average relative abundance of genus level taxa in the 24 infants at different time points. In case the taxonomic assignment could not be made at genus level, the lowest classifiable taxonomy assignment is used, and the unidentified genus is indicated with "_g".
Milk HMo content and changes during lactation. HMO profiles of each breastmilk sample were analysed (Table S3a), and samples from six mothers revealed consistent lack of 2′FL and low levels of Lacto-N-fucopentaose I (LNFP I), Lacto-N-difucohexaose I (LNDFH I) and Difucosyllactose (DFL) (non-secretors). PCA of milk HMOs showed a clear separation of samples in relation to collection time point postpartum and secretor status (Fig. 2a). This finding was confirmed with RDA, which showed that collection time point and secretor status had a significant effect on the types and levels of the HMOs present, explaining respectively 16.7% and 2.5% of variation in the data. When analysed separately for each time point, neither the mode of delivery nor maternal stress measured via saliva cortisol, were significantly associated with the HMO levels in mothers' milk (data not shown). The average concentration (µg/ml) of the HMOs measured in the milk decreased with time of lactation, except for the concentrations of 3FL and LNFP III, which were secreted at higher amounts at later time points (Fig. 2b, Table S3a). After adjusting the HMO concentrations for the literature-based averages 19 of the amounts of milk consumed (480 g, 580 g and 630 g at weeks two, six and 12, respectively), the estimated average HMO intake remained stable in the first three months of life (Fig. 2c). correlation between milk HMos and infant faecal microbiota. The estimated daily intake (g/day) of each HMO in each infant, and the relative abundance of OTUs for the corresponding faecal samples were correlated to examine the potential link between the levels of ingested HMOs and the abundance of different microbial taxa in infant faeces (Fig. 3). The strongest positive correlations were detected between few of the Staphylococcus OTUs and milk levels of 6′SL, Sialyl-lacto-N-tetraose c (LST c), Sialyl-lacto-N-tetraose a (LST a), LNH, Lacto-N-neohexaose (LNnH), and LNFP I. Streptococcus OTUs were positively correlated with 3FL, LNFP III, and Para-lacto-N-hexaose (pLNH). Positive correlations were also found between OTUs within the genus Actinomyces and 3′SL, and Enterococcus OTUs and 3FL, Lacto-N-fucopentaose II (LNFP II), LNFP III, Lacto-N-fucopentaose V (LNFP V) and LNDFH II. No significant positive correlations were found between the most predominant Bifidobacterium OTU 1263 and any of the milk HMOs. Spearman correlation analysis at each separate time point (Fig. S3) showed that at two weeks of age three bifidobacterial OTUs were positively associated with LNH, pLNH and LNFP I, and the main Bifidobacterium OTU 1263 was negatively associated with LNFP II and LNFP III (Fig. S3a). At six weeks three different low abundance bifidobacterial OTUs were positively associated with LNFP V, LNnH, 6′SL, Sialyl-lacto-N-tetraose b (LST b) and LST c (Fig. S3b), and at 12 weeks LNnH, LNFP II, and LNDFH II showed a low positive correlation with two low abundance bifidobacterial OTUs. Overall the associations did not exceed ±0.6, there was no consistency in the type (positive or negative), or strength (passing threshold of ±0.3, p < 0.05) of associations between specific OTU-HMO pairs when data from all time points was combined, or when individual time points were analysed separately. correlation between infant faecal microbiota and HMos excreted in infant faeces. Undigested milk HMOs were secreted in infant faeces and their concentrations varied between infants and across study time points (Table S3b). Our hypothesis was that aside of being directly dependent on the amounts of the HMOs ingested with the milk, the faecal HMO concentrations could serve as an indicator of the utilisation of different HMOs by infant GI tract microbiota. PCA (Fig. 4a) showed that samples clustered by infant age, and to lesser extent also by the maternal secretor status. RDA confirmed these findings showing a significant (FDR < 0.05) effect of infant age (7.5% explained variation) and maternal secretor status (7.3%) on the HMO concentrations in infant faeces (Figure not shown). The concentrations of all HMOs in faeces decreased with infant age, with the exception of 3FL (Fig. 4b). Spearman correlation between the faecal HMO concentrations and faecal microbiota showed that the majority of statistically significant negative associations were between bifidobacterial OTUs and thirteen different HMOs (Fig. 5). The main Bifidobacterium OTU 1263 was negatively correlated with nine different HMOs, of which LNH, LNFP V, and the LNT/LNnT group showed strongest correlations, suggesting an important role of Bifidobacterium OTU 1263 in consumption of these HMOs in the infant GI tract. Positive associations were also observed between various OTUs of Streptococcus, Staphylococcus, Clostridium and Escherichia-Shigella and thirteen different HMOs. Spearman correlation analysis at each separate time point showed that at two weeks of age the strongest negative correlations were found between various HMOs and highly abundant Bifidobacterium OTU 1263 and 681 (Fig. S4a). At week six, an additional six bifidobacterial OTUs showed negative correlations with pLNH, LNH, LNFP II and III. At both time points Bifidobacterium OTU 1263 was the only one correlating negatively with 3FL, possibly highlighting the unique link between this HMO and the major Bifidobacterium OTU during the initial stages of the development of GI tract microbiota (Fig. S4b). At 12 weeks of age, the HMO concentrations in faeces were very low or no longer within the detectable range, leading to overall less clear correlation patterns between bifidobacteria and faecal HMOs. Still, at 12 weeks of age Bifidobacterium OTU 1263 correlated strongly with 2′FL, LNDFH II and LNFP V (Fig. S4c).
In order to account for the initial availability of the different HMOs in the milk, we calculated ratios between each faecal and milk HMO, for each mother-infant pair at each time point. The resulting ratios were then used to estimate the consumption level for each HMO as either "high", "medium" or "low" based on tertiles. RDA showed that high consumption was associated with older infant age and higher levels of bifidobacteria. Low consumption of pLNH, 2′FL, LNH, LNnH, LNTandLNnT, LNFP I and V, 3′SL and LNDFH I, medium consumption of LNDFH II, LNFP II and 6′SL, and high consumption of pLNH, LNH, LNnH, LNFP I and V, LST a, LNTandLNnT, and 3′SL www.nature.com/scientificreports www.nature.com/scientificreports/ were significantly associated with the microbial composition (FDR < 0.05; Fig. 6). Overall, high consumption was detected in association with various bifidobacterial OTUs including the most predominant Bifidobacterium OTU 1263, and several OTUs within genera Parabacteroides, Escherichia-Shigella, Bacteroides, Actinomyces, Veillonella, and Erysipelotrichaceae Incertae Sedis (Figs. 5 and 6). The result was confirmed when relative abundance of OTUs between infants assigned into a low and high consumption groups for each HMO were compared (Table S4) showing significantly higher relative abundance of various bifidobacterial OTUs, and in most cases the most  www.nature.com/scientificreports www.nature.com/scientificreports/ predominant OTU Bifidobacterium 1263 in the high consumption group (Fig. 7). Interestingly, high consumption of 6′SL, DFL, LST a, and LST c was not significantly linked with the presence of bifidobacteria, but instead was linked with other bacterial groups (p < 0.05; Table S4).

Discussion
In the cohort described in this study we observed directional changes in both HMO concentrations in breast milk and in infant faces, and in the composition of infant faecal microbiota. HMO concentrations differed between individuals and depended on lactation time and infant age. Earlier studies showed that faecal microbiota structure is highly dynamic during infancy and is associated with multiple factors, including diet [20][21][22][23] . The variation in the faecal microbiota of infants in our study could be explained by infant age, sex, place and mode of delivery and certain milk HMOs, namely 6′SL, 2′FL, 3FL, 3′SL, LNDFH II (p < 0.05). The effect of factors varied with age. Mode of delivery and breastmilk LNFP III concentrations showed a significant association with faecal microbiota structure at two weeks of age. The important role of the mode of delivery in the initial seeding of the GI tract has been previously reported 24,25 and has been linked to various health outcomes, both in infancy and beyond 26,27 . At six weeks the significant effect of the mode of delivery could still be detected in our data, and also infant sex and milk 3′SL seemed to play a role. A recent study using animal models showed that sex specific, microbiota-independent differences in immunity may lead to the selection of a sex-specific GI microbiota in adult germ free mice 28 . Sex-related differences in faecal microbiota have been also reported in adults 29 , in pre-term 30 , and term infants 31 , yet the timing and the possible mechanisms that might underlie the sex differentiation of GI microbiota in healthy full term infants remain largely unknown. As male infants tend to have a higher daily milk intake, it is possible that they receive higher doses of microbial and HMO components from the mother's milk 32 . At 12 weeks of age, infant sex and LNH were significantly associated with microbiota, however, the significant association between microbiota composition and mode of delivery was no longer present. With the accumulating evidence linking mode of delivery with various health effects later in life 26,27,33 , it is likely that the microbiota related programming of the host happens soon after birth, or during specific "windows of opportunity", thus, even when microbiota recovers to its normal state, the long term health effects of such disturbance might persist throughout life 33 . www.nature.com/scientificreports www.nature.com/scientificreports/ One of the crucial factors shaping the development of GI microbial community during infancy is the type of feeding that an infant receives, with formula-and breastfed infants showing very different microbial profiles [34][35][36][37][38] . Human milk is a source of prebiotic HMOs which support microbial colonisation in the infant GI tract 7,16,39 . The HMO content in breastmilk is variable and influenced by maternal, environmental and infant feeding practices 12 . The HMO content is highest in colostrum, and the concentrations of HMOs decrease in mature milk [8][9][10] . Our data agrees with these findings showing that concentrations of the HMOs measured per mL of milk varied between mothers, and decreased between two and 12 weeks of lactation, except for 3FL and LNFP III, which increased in concentration as lactation progressed. Earlier studies showed that the daily intake of milk increases in the first months of life 19 likely compensating for the decreased HMO concentrations. Considering the daily milk intake, we could conclude that in our study the amounts of ingested HMOs remained relatively stable during the first 12 weeks of life. The exception was in the intake of 3FL and LNFP III which gradually increased in time. Interestingly, this increase in the 3FL and LNFP III intake corresponded with the increase in the faecal concentrations of these two HMOs, suggesting that, on average, their supply likely exceeded the ability of the infant GI tract microbiota to utilise these two HMOs. Furthermore, the role of these structures might go beyond their prebiotic effect, as both 3FL and LNFP III can act as decoys against pathogen binding and infection [40][41][42] .
One of the signature bacterial groups found in faeces of breastfed infants is the genus Bifidobacterium 2 . In vitro studies showed that Bifidobacterium bifidum JCM1254 39 , Bifidobacterium longum subsp. infantis, Bacteroides fragilis and Bacteroides vulgatus 16 grow well on HMOs as sole carbon source. Thus, we expected to find positive correlations between certain breastmilk HMOs and the bifidobacterial OTUs. We did not measure the absolute abundance of bacteria in infant faeces and are aware of the limitation of using relative abundance data. However, culture-based studies showed that the total faecal bacterial load, as well as Bifidobacterium counts tend to increase in the first month of life 35 . Thus, the observed increase in the number of sequencing reads and the relative abundance of bifidobacteria in our data likely reflect the actual increase in the abundance of this group (Fig. S1). We hypothesized that quantity of HMOs might selectively promote growth of either the primary or secondary HMO degraders, leading to increase in their abundance within the microbial community. However, when using data from the three time points combined, we saw an opposite trend -as the predicted daily intake of most HMOs was stable in time, the relative abundance of bifidobacteria was increasing. The analyses repeated for individual time points showed few positive correlations and the results varied between time points (Fig. S3a-c). This could be due to the fact that, in addition to individual HMOs promoting growth of selected microbes capable to utilise this carbon source, there may be other mechanisms controlling the microbial community structure. For example, presence of other breastmilk components, such as lysozyme, secretory IgA and other endogenous factors can suppress growth of certain members of the community and thus, indirectly allow other bacterial species to dominate the infant GI tract ecosystem 35 . Furthermore, the possible effect of HMOs on microbiota may already start in the breastmilk itself 43 . A recent study on human milk investigated the associations between the HMO content and microbiota composition in colostrum and reported strong positive correlations between different HMOs and various microbial groups, including streptococci, staphylococci, enterococci and bifidobacteria, in particular Bifidobacterium breve and LNFP III 44 . In the present study, RDA showed a significant association of milk LNFP III with infant faecal microbiota in all time points combined. There was a strong positive association www.nature.com/scientificreports www.nature.com/scientificreports/ of milk LNFP III with OTUs belonging to the genera Veillonella, Enterococcus and Streptococcus (Fig. S3). Positive associations at all time points for breastmilk 3′SL and unidentified OTUs within family Enterobacteriaceae and Actinomyces 695 were also observed. 6′SL was positively associated with clostridia -Clostridium 1789 at week two and Clostridium 1639 at week 12. Finally, LNFP III was positively associated with Enterococcus 1698 at six and 12 weeks. At two and at six weeks Lachnospiraceae Incertae Sedis 876 was negatively associated with LNDFH II, Lactobacillus 1718 was negatively associated with 3′SL, and Bifidobacterium 1295 with LNFP I. Studies on mature breastmilk microbiota and the microbial transfer of microbiota from mother to infant show that breastmilk contains a distinct microbial community and that breastfed infants receive on average nearly 30% of the bacteria from breastmilk and 10% from areolar skin in the first 30 days of life 45 . The study also concluded that the association was lower in older infants, and it was proportional to the frequency of breastfeeding that an infant received 45 . Thus, it is then likely that some of the correlations observed here were due to a combined effect of the HMOs modulating the microbiota of the mother's milk and the infant GI tract, as well as due to a direct transfer of bacteria during breastfeeding 43 .
Infant GI microbiota plays an important role in energy metabolism via utilising otherwise indigestible HMOs. Our data showed that the average concentrations of faecal HMOs decreased with age, suggesting that microbiota of older infants was more adapted and efficient in degrading these compounds (Fig. 5b). Furthermore, we noted that the increase in efficiency was correlated with the increase in the relative abundance of several bifidobacterial OTUs, but also Parabacteroides, Escherichia-Shigella, Bacteroides, Actinomyces, Veillonella, and Erysipelotrichaceae Incertae Sedis (Fig. 6). High ability to degrade a wide range of HMOs was associated with higher relative abundance of one or more Bifidobacterium OTUs confirming the important role of this bacterial group in the HMO metabolism (Fig. 6). Thirteen faecal HMOs negatively correlated with nine different Bifidobacterium OTUs, and the highly abundant Bifidobacterium OTU 1263 was negatively correlated with nine different HMOs in faeces, especially LNH, LNT and LNnT and LNFP V. Aside of bifidobacteria, members of the genus Bacteroides were significantly more abundant in infants who were good degraders of 2′FL, LNFP I, II, V, and pLNH, and Parabacteroides in the high degraders of 3FL, LNFP V, LNH, LNT and LNnT, indicating that these microbial groups might have a mutualistic or symbiotic relationship degrading those compounds. In addition, Halomonas, Enterococcus, Lactobacillus, Staphylococcus, Suterella, Varibaculum, Veillonella, Streptococcus, Actinomyces, and Lachnospiraceae Incertae Sedis were also associated with degradation of the same HMOs as bifidobacteria. Interestingly, high levels of degradation of 6′SL, DFL, LST a and LST c were not associated with high levels of any of the bifidobacterial OTUs, but with various OTUs belonging to Bacteroides, Streptococcus, Varibacullum (6′SL), Actinomyces, Clostridium, Collinsella and Streptococcus (LST c), and Haemophilus, Veillonella (DFL), and Lachnospiraceae Incertae Sedis, and Halomonas (LST a).
Negative associations were also observed for LNFP II and Bacteroides, and for LNFP V and Parabacteroides suggesting the role of these bacteria in the HMO degradation. The fact that Bacteroides and Parabacteroides (formerly also Bacteroides) were identified in our analysis is in line with earlier studies showing growth of Bacteroides spp. on selected milk glycans by activating the mucus degradation pathway 46 . Finally, LNDFH I in both milk and faeces was negatively associated with Enterococcus OTU 1702, but the association was stronger in faeces. Even though in vitro studies showed that in a monoculture Enterococcus was not able to grow on milk HMOs 16 , another study showed that this group was found in breastmilk 44 , that its abundance in infant faeces could be predicted from the maternal HMO profile and that it was positively correlated with the abundance of Bifidobacterium, Streptococcus and Veillonella 7 . One of the suggested explanations was that Enterococcus can cross feed on HMO fermentation products or HMO breakdown by-products that are released in the ecosystem by HMO degrading bifidobacteria or Bacteroides spp. 7 .
The correlation analysis of infant faecal HMOs and infant faecal OTUs for all time points combined also detected numerous significant positive associations between various HMOs and Streptococcus, Staphylococcus, Escherichia-Shigella, and Clostridium OTUs (Fig. 5). In both, milk and faeces LST c, 6′SL and LNnH showed significant correlation with staphylococci, while LNFP III and 3FL were positively correlated with streptococci. Earlier studies showed that neither Streptococcus, Staphylococcus, Escherichia-Shigella 7,16 , nor Clostridium 16 could effectively utilise and grow on milk HMOs. However, all these bacterial groups are members of the microbiota of breastmilk and areolar skin 45,47,48 , and the positive link might be due to breastfeeding practises, for example with more frequent feedings resulting in higher ingested doses of the bacteria and the HMOs. If the HMOs are not well digested, the positive associations may still persist in the faeces.
Our study has few important limitations. The number of participating mother-infant pairs was relatively small and future studies should be done on larger sample sizes to verify our findings. Larger sample size would also allow for stratification of data based on maternal secretor status. Another limitation of this study is the lack of measurements of the actual volumes of the breast milk ingested and the mass of the faeces produced by each infant in a 24-hour sampling period. Therefore, our calculations had to be based on estimates from other studies that measured the daily milk intake in infants of similar ages. Actual intakes, although likely increasing over time, can vary greatly between infants, and thus, our calculated values are only an estimate and should therefore be considered as such. Measurements of volumes or weights of breastmilk samples and faecal samples should be included in future studies. Finally, with our methods we could not identify species or strains of the bacteria, nor could we quantify their actual numbers which could further support our findings.
In this study we found only weak associations between selected breastmilk HMOs and faecal microbiota community structure in breastfed infants. However, we found a strong link between changes in levels of various HMOs in infant faeces and specific microbial groups, in particular different species of Bifidobacterium. Earlier studies showed that different bifidobacterial species vary in their ability to break down HMOs, and some species can degrade HMOs without experiencing a detectable population growth 49,50 . Thus, including metatranscriptome or metaproteome analyses in this set would have been very helpful in understanding the community dynamics in regard of HMO metabolism in the infant GI tract. Our findings could provide the basis for assembling Scientific RepoRtS | (2020) 10:4270 | https://doi.org/10.1038/s41598-020-61024-z www.nature.com/scientificreports www.nature.com/scientificreports/ simple synthetic communities to study microbial interactions and community structure changes, which are centred around degradation of different HMO structures. In vitro fermentation studies incorporating purified compounds would also help to eliminate confounders, such as presence of milk's own microbiota and presence of milk components, which have a regulatory effect on microbiota in both, milk and in the infant GI tract.

conclusion
Faecal microbiota of breastfed infants during the first 12 weeks of life is highly diverse, dynamic and influenced by age and other factors. The effect of mode of delivery disappeared after six weeks of age, whereas the effect of infant sex became detectable over time. Overall, microbiota development in this cohort followed a normal colonization pattern resulting in faecal microbial communities dominated by Bifidobacterium, in particular the most predominant Bifidobacterium OTU 1263. Breastmilk HMO analyses showed that the composition of the 18 HMOs that were measured varied between mothers and throughout the duration of lactation. In our analysis we did not observe strong and consistent positive correlations between the HMOs in maternal breastmilk and specific microbial OTUs including bifidobacteria in infants' faeces. Thus, we believe that HMO composition is only one of many factors regulating colonization and structure of the infant GI microbial community. However, the findings of this study supported the key role of bifidobacteria in the infants' ability to utilize most of the measured HMOs, in addition to indicating the role of other microbial taxa in the degradation or metabolism of specific HMOs.

Methods
Sample collection. The infants included in this study were born at term (after 37 weeks of gestation) from single pregnancies, were healthy and had no congenital abnormalities related to growth, did not receive oral or systemic antibiotic treatments, were exclusively breastfed during the study period, and had the maternal breastmilk and infant faecal samples available from each study time point. The samples originated from the BINGO (Dutch acronym for Biological Influences on Baby's Health and Development) cohort, which is an ongoing longitudinal study investigating prenatal predictors of infant health and development. This study and all its experimental protocols were approved by and carried in accordance with the ethical committee of the Faculty of Social Sciences of the Radboud University [ECSW2014-1003-189]. Parents and/or legal guardians of all participating infants were asked to sign an informed consent form and were free to withdraw from the study at any point. The study design and infant recruitment criteria can be found at http://www.bingo-onderzoek.nl/deelname/. Both, the infant faecal samples and the maternal breastmilk samples were collected between years 2015-2016, within a 48-hour period, by the mothers at home, at two, six, and 12 weeks post-partum. Breastmilk samples (approximately 20 ml) were collected in the morning before feeding the infant, into clean, sterile collection cups as previously described 51 . Mothers were asked to wash hands, breasts and nipples before collecting the sample by hand, or in case of mechanical collection, to first sterilize the breast pump compartments by boiling. Stool samples were collected from infants' diapers using sterile stool collection vials (80 × 16.  (Table S1). The HMOs were extracted, purified by solid phase extraction (SPE), and quantified by using porous graphitized carbon-ultra high-performance liquid chromatography -mass spectrometry (PGC-UPLC-MS) 52-54 , or by high performance anion exchange chromatography with pulsed amperometric detection (HPAEC-PAD) 55 .
DNA extraction, amplification of 16S rRNA genes and sequencing. Total bacterial DNA was extracted from 0.1-0.15 g of faeces using the double bead-beating procedure and the Maxwell ® 16 Total RNA system (Promega; Wisconsin, USA) customized with Stool Transport and Recovery Buffer (STAR; Roche Diagnostics Corporation, Indianapolis, IN) for the use with the DNA samples, as previously described 56 . The resulting DNA templates (20 ng) were used for subsequent PCR amplification of the V4 region of 16S ribosomal RNA (rRNA) genes using barcoded primers 515F-n (5′-GTGCCAGCMGCCGCGGTAA-) and 806R-n (5′-GGACTACHVGGGTWTCTAAT) following conditions described previously 56 . Seventy unique barcode tags were used in each library 57 . Negative control blanks were included during DNA extraction and PCR to check for any possible contamination. The negative controls samples gave no PCR products when visualised on agarose gels, and subsequently were not sequenced. All remaining PCR products were purified, 100 ng of each barcoded sample was added to an amplicon pool and the pools were adjusted to 100 ng/µL final concentration. The libraries were sent for adapter ligation and Illumina HiSeq 2000 sequencing at GATC-Biotech, Konstanz, Germany 56 . Data analysis. The 16S rRNA sequencing data analysis was carried out using the NG-Tax analysis pipeline with standard parameters 57 . Filtered libraries contained only read pairs with perfectly matching barcodes that were subsequently used to separate reads by sample. OTUs were assigned using an open reference approach and the SILVA_111_SSU 16S rRNA gene reference database (https://www.arb-silva.de/) 58 . Microbial composition data was expressed as a relative abundance of each OTU obtained with NG-Tax. Milk  www.nature.com/scientificreports www.nature.com/scientificreports/ variability due to naturally occurring differences in milk dilutions. Changes in HMO concentrations between study time points were assessed using one-way repeated measures ANOVA for HMOs for normally distributed data sets as assessed with Shapiro-Wilk's normality test. Differences between not-normally distributed HMO sets were tested using nonparametric Kruskal-Wallis test. We estimated average daily volumes of ingested breastmilk based on literature data 19 to be 480 g at week two, 580 g at week six, and 630 g at week 12 of life. These values were then used to calculate the daily amount of each HMO consumed by infants in our study at each time point.

Statistical analyses.
Microbial composition data was expressed as relative abundance (RA) of each OTU obtained in the NG-Tax pipeline. OTUs which had a prevalence of less than 5% across all samples were removed and their values were summarized as "Other OTUs". Alpha diversity index determinations (Shannon, Chao1, and PD Whole Tree) were carried out in QIIME on rarefied data with OTU cut-off of 3380 reads per sample 59 . Spearman correlations were calculated using R to evaluate associations between OTU members of the faecal microbial community, and between the faecal OTUs and milk or faecal HMO concentrations. We only reported associations which passed the correlation threshold of ±0.3 and significance of p < 0.05. Unconstrained (PCA) and constrained (RDA) multivariate analyses were carried out in Canoco5 on log transformed HMO concentrations data and on microbial OTU level relative abundance data with the significance assessed using the Monte Carlo permutation test at 499 random permutations 60 . For better image clarity, only the top twenty best fitting taxa were displayed on the PCA and RDA plots. These taxa were those for which the highest percentage of variation in the relative abundance data was explained by the ordination axes. The vectors corresponding to these taxa point towards the ordination plane where samples containing higher relative abundance of these taxa are located and the lengths of these vectors equal to R-squared calculated by dividing the taxa scores by their SD 60 . The explanatory variables used in the RDA included infant age at the time of collection, estimated amounts of the 18 milk HMO ingested during a 24 h period (mg/24 h) or faecal HMOs (µg/g of faeces), infant sex, place and mode of delivery, maternal secretor status, and if an infant was sick at the time of sample collection, as recalled by the mother. Degradation of each breastmilk HMO was calculated as a ratio of HMO concentration in infants' faeces and the concentration of the same HMO measured in mothers' milk. If a given HMO was not detected in milk, the consumption score was not included in the analysis, and if the concentration in faeces exceeded the amount detected in milk, the infant was assigned to the "low" category for that HMO. Resulting values (ratios) and the natural breaks in the data were then used for separating and assign infants to either a "low", "medium", or a "high" consumption category for each HMO at each time point. The association between faecal microbiota composition and the assignment of each infant to a "low", "medium", or "high" consumption category for each HMO were investigated with RDA in Canoco5, with significance assessed using a permutation test 60 . Differentially abundant OTUs were then identified by comparing the microbial abundance data of infants assigned into "high" and "low" consumption groups and for each individual HMO using Kruskal-Wallis analysis in QIIME with a significance cut-off set at FDR < 0.05 59 .
Nucleotide sequences. BINGO data sets cannot be made publicly available due to the data being part of an ongoing longitudinal study. Parts of the data are available to the research community for scientific collaborations upon request to Prof. Dr. C. de Weerth at: Department of Cognitive Neuroscience, Donders Institute for Brain, Cognition and Behaviour, Radboud University Medical Center, 6525 HR Nijmegen, the Netherlands, e-mail: Carolina.deWeerth@radboudumc.nl.