Human milk oligosaccharide composition and associations with growth: results from an observational study in the US

Background: Breast milk is the recommended source of nutrients for newborns and infants. Human milk oligosaccharides (HMO) are the third most abundant solid component in human milk and their composition varies during lactation.

Objectives: Our objective was to investigate longitudinal and cross-sectional changes in HMO composition and whether these changes were associated with infant growth up to 24 months of age. Associations with maternal characteristics were also investigated.

Methods: 24 HMOs were quantified in samples taken at 2 weeks (n = 107), 6 weeks (n = 97) and 3 months (n = 76), using high performance liquid chromatography. Body length, weight, and head circumference were measured at 8 timepoints, until 24 months. Clusters of breast milk samples, reflecting different HMO profiles, were found through a data-driven approach. Longitudinal associations were investigated using functional principal component analysis (FPCA) and used to characterize patterns in the growth trajectories.

Results: Four clusters of samples with similar HMO composition were derived. Two patterns of growth were identified for length, body weight and head circumference via the FPCA approach, explaining more than 90% of the variance. The first pattern measured general growth while the second corresponded to an initial reduced velocity followed by an increased velocity (“higher velocity”). Higher velocity for weight and height was significantly associated with negative Lewis status. Concentrations of 3’GL, 3FL, 6’GL, DSNLT, LNFP-II, LNFP-III, LNT, LSTb were negatively associated with higher velocity for length.

Conclusion: We introduced novel statistical approaches to establish longitudinal associations between HMOs evolution and growth. Based on our approach we propose that HMOs may act synergistically on children growth. A possible causal relationship should be further tested in pre-clinical and clinical setting.

Introduction

Human milk oligosaccharides (HMO) are the most abundant solid component of human milk, after fat and lactose, with an estimated concentration ranging between 5 and 15 g/L (1). Although largely not digestible by the infant’s gut, HMOs are known to be associated with a range of biological functions, such as gastrointestinal development (2–8), protection against infection (9–11) and inflammation (12) and more recently in neurodevelopment (13). Breastmilk composition, and in particular the concentration of HMOs, is known to evolve during the lactation period (1, 14, 15); the total concentration of HMOs has been observed to decrease, as well as the concentration of most of the individual HMOs, with the notable exception of 3FL which increases over time (1, 15, 16). Berger et al. (17) have speculated that these dynamic changes in HMO concentrations may reflect their involvement in biological functions beyond the first few months.

The concentration of HMOs in breast milk varies among individuals and particularly depends on the expression of α1-2-fucosyltransferase (FUT2) and α1-3,4-fucosyltransferase (FUT3) genes which determine the secretor (Se) and Lewis (Le) status, respectively (18). 2′-fucosyllactose (2’FL) is absent in the milk of non-secretor women, while it is the most abundant HMO in the milk of secretors. While maternal genetic polymorphisms for the Se and Le genes explain, to a large extent, the HMO variability between mothers (18), additional maternal or infant factors like pre-pregnancy BMI, age, mode of delivery and gestational age may contribute as well to the interindividual variability (15).

The first objective of this exploratory work was to characterize the temporal and inter-individual variability of the HMO concentrations in a US cohort over the first 3 months of lactation. 24 HMOs were measured at three time points; although approximately 150 HMOs have been identified in the literature, it is known that the HMOs used in our study account for roughly 90% of the total HMO (1). We developed a novel data-driven approach to cluster breast milk samples based on the concentrations of HMOs, and tested if these clusters were associated to maternal factors, such as mode of delivery or pre-pregnancy BMI. The data-driven clustering was based on the abundances of all the measured HMOs, assigning samples to milk types sharing similar HMO profiles, and we compared the resulting clusters to the milk types defined by the presence/absence of specific HMOs (e.g., 2’FL and LNFP-I) as proxies for the Se, Le status (18). Van Leuween (19) observed that it may not always be correct to determine the milk type based on HMOs composition at one single time-point, and in our analysis we also considered the intra-individual variability, by describing trajectories of individual HMOs within each cluster.

As a second objective, we investigated the effect of HMO concentrations on growth parameters (head circumference, body weight and body length, measured at 8 time points up to 24 months of age). HMOs are hypothesized to affect growth through several mechanisms, including a possible role in functional gut maturation, therefore improving nutrient absorption (20). Moreover, in some studies, multiple HMOs (eg 2’FL, and 3FL) have been detected in the urine and plasma of breastfed infants, but not in formula-fed infants (21), suggesting the possible existence of systemic effects, yet to be identified (14, 16, 22–27).

The biological functions of many, if not most, HMOs remain unknown. Selected HMOs are increasingly added to infant formulas for their beneficial health benefits (6). It is therefore important to understand better these benefits, whether they might be supported by synergistic effects and what are the natural variations of the individual HMO concentrations over the course of lactation. This knowledge might help to target the HMOs that are the most relevant at critical windows of infant development. In turn, such investigations pose methodological challenges for the data analysis, that the present study intended to address.

Patterns of development are often non-linear, therefore linear mixed effect models that are often used in modelling growth trajectories may not describe the data accurately. Non-linear parametric approaches are an alternative to linear mixed models, however, they rely on parametric assumptions that might not be realistic. In addition, the sparsity often found in real data is a challenge for fitting the data into those models (23). Functional data analysis offers a non-parametric, non-linear approach to model sparse, non-linear longitudinal data. The FPCA approach allows extracting the patterns of growth from a set of trajectories and has been successfully applied in many domains, including neurodevelopment (28) and gene expression (29, 15).

Materials and methods

Endpoints

Table 1 summarizes the timepoints available for the milk samples and the body measures.

TABLE 1

Table 1. Overview of time points.

Description of study population

Children for this study were selected from the observational breastfeeding arm of a prospective longitudinal randomized control trial. Results for the randomized arm were analyzed in a previous publication (30). Infants for the trial were recruited and followed over 24 months at two study sites in the United States (Rhode Island Hospital in Providence, RI, and Pennington Biomedical Research Center, in Baton Rouge, LA). Maternal screenings were performed during the third trimester of pregnancy up to and including post-delivery. Following written informed consent (screening visit), sociodemographic information, medical and family histories were collected, as well as a physical and neurological examination of the infant. Withdrawal from the study was possible at any point and with no further evaluations and any additional data collection. The research ethic boards at both clinical sites approved the protocol.

In total, 107 breast milk samples were available for analysis. The demographic characteristics of the mothers are summarized in Table 2. 61 (57%) children were female, 46 (43%) were male. Mean age of the mothers at recruitment ranged from 19 to 43 years, with an average of 32 years. Most mothers were white (70%) BMI before pregnancy ranged from 19.1 to 39.9, with an average of 27.2 and a median of 26. 72% of children were born through vaginal delivery.

TABLE 2

Table 2. Maternal characteristics.

Supplementary Table 1 summarizes the anthropometric characteristics of the breastfed children involved in the study.

Exclusion criteria were: (i) Birth>41 weeks +6 days gestation as reported in medical record when available; (ii) Birth Weight < 2000 g or small for gestation age (birth weight less than the 10th percentile for the gestational age) or large for gestational age (weight, length, or head circumference above the 90th percentile); (iii) Any unsafe psychopharmacological treatment of mother using prohibited medications during pregnancy or lactation as assessed by medical interview. This included anticonvulsants, antidepressants, benzodiazepines, cytotoxic drugs, dopamine agonists, opioids.

Milk samples

Breast milk was sampled from mothers longitudinally at defined study visits, each time between 10 AM–12 PM from the right breast using a hospital grade electric breast pump. Mothers were asked to empty the right breast approximately 2 h prior to milk sampling and the time of milk sampling the complete breast was emptied using a pump (single full breast milk sampling methodology). Only a fraction of the collected milk was aliquoted for research and the rest was returned to mother to feed the baby at a later time.

Data on feeding mode and introduction of solid food were collected via questionnaires at baseline, 6 weeks, 3, 6, 9, and 12 months. Up to at least 4 months of age children received breastmilk for more than 90% of nutritional intake and no infant formula or solids for more than 10% of the nutritional intake (approximately 6 feeds per week).

HMO data (Table 3) were available for N = 107, 97, 76 individuals at 2–5 weeks, 6 weeks, 3 months (visits V0–V2), respectively. 76 mothers had complete HMO data for all three visits.

TABLE 3

Table 3. Overview of all HMOs measured, with chemical structure and abbreviation.

HMOs were analyzed by ultra-high performance liquid chromatography with fluorescence detection (UHPLC-FLD), according to the method of Austin and Benet (31). We refer to Table 3 for the full list of measured HMOs and their respective abbreviations. 2’FL, 3FL, 3’SL, 6’SL, LNT, LNnT and LNFP-I were quantified against standards of analytical quality all other HMOs were quantified against maltotriose assuming equimolar response factors. Mother’s genotypes to define Secretor and Lewis status were not available.

Statistical analysis

We first calculated the sum of all HMO concentrations at each time point in each sample and tested for significant differences between visits V0 and V1, and between V1 and V2. We compared the Shannon diversity of the HMOs between timepoints. The Shannon diversity index is low when all HMOs are present in similar concentrations. The more unequal the abundance of HMOs, the larger the corresponding Shannon entropy.

In this work, we analyzed both the proportion of each HMO, expressed as % of the sum of the 24 measured HMOs, and their concentration expressed in mg/L. Since the volume of milk per feeding was not measured, total amounts of HMOs or other milk components are not known, and only relative information was available for analysis.

We applied a data-driven method to assign mothers to clusters of HMO composition. The number of clusters was data-driven as well. Basically, we clustered together samples sharing a similar HMO profile. Our approach was based on network theory (32), where each node of the network corresponds to a milk sample and two nodes are connected based on the similarity between their HMO profiles. We first calculated the Aitchison (33) distances between all the pairs of samples, all subjects and all timepoints. The resulting matrix of distances allows us to identify the closest neighbor to a given sample. The edges in the network were weighted, with weights inversely proportional to the distance between the samples. Therefore, two samples with similar HMO composition have a short distance and therefore will be connected with a high weight in the network. We applied the leading eigenvalues method (34) to identify clusters in the network. These clusters correspond to groups of samples sharing a similar HMO composition. We applied a Kruskal-Wallis test to each HMO to test whether its concentration was significantly different between the clusters. We also made qualitative comparisons with clusters of HMOs that can be derived from presence or absence of 2’FL and LNFP-II (3, 6).

We tested for differences between HMO clusters for several covariates: mother’s BMI before pregnancy, BMI status, mode of delivery, ethnicity, and gestational age. BMI status was defined as ‘overweight’ (OW) for a BMI ≥25 kg/m (2), ‘obese’ (OB) if BMI ≥30 kg/m², ‘normal weight’ (NW) otherwise. Fisher’s exact tests were used to test the independency between two categorical variables. Permutational analysis of variance (PERMANOVA) was used to test the equality of the overall HMO composition between groups (e.g., clusters). Kruskal-Wallis test was used to check the independence between a continuous and a categorical variable, followed by a Dunn post-hoc test.

We calculated all pairwise correlations between the individual HMOs, at each timepoint, and visualized them as a heatmap, to summarize and display patterns in the data. In particular, the existence of strong correlations between pairs of HMOs is an important background information when interpreting associations arising from the application of univariate models, and suggests a possible synergistic effect of several HMOs acting in combination.

All correlations reported in this study were Spearman correlations, and a Benjamini-Hochsberg correction for multiple testing was applied whenever applicable, using a false discovery rate of 5%.

Previous work (14, 16) has explored the influence of individual HMOs on growth, using linear mixed models with growth parameters as response variable and the HMO level as predictor. We propose an alternative approach derived from functional data analysis (22–27) to describe the growth trajectories and their main modes of variation; this approach allowed us to identify distinct growth patterns in the data, that could not be properly detected using linear models. These patterns were then investigated in association with the clusters derived in the first part of the analysis.

We investigated the associations between the HMO composition and the child’s growth over time.

We plotted the growth data separately for boys and girls, using the WHO reference values.¹ Growth trajectories were analyzed within the methodological framework of functional data analysis (23). More specifically, we applied Functional Principal Component Analysis (FPCA) to describe trajectories using a limited number of numerical parameters, the FPCA scores (11). FPCA models longitudinal data as samples from smooth curves, so that the time-varying trait of the i-th subject admits a Karhunen-Loeve expansion

$X_i\left(t\right)=\mu \left(t\right)+{\displaystyle \sum }_{k\ge 1}{\xi }_{i,k}{\varphi }_k\left(t\right)$

where $\mu \left(t\right)$ represents the average trajectory and the ${\xi }_{i,k}$ are the FPCA scores. The sum in the above formula can be truncated to a fixed number of terms, and the fraction of variance explained will depend on how many terms are kept in the sum. We applied FPCA separately to the length, weight, and head circumference data.

For our applications, based on the percentage of variation explained, we selected the first 2 scores as descriptors for the trajectories, FPCA1 and FPCA2. We therefore effectively achieve a dimensionality reduction, where a trajectory, a priori defined by 8 time points, can be effectively described by just 2 parameters.

We then tested whether the FPCA scores were associated with the HMO clusters, with the single HMO concentrations or with the maternal baseline characteristics.

We then further investigated which of the HMOs concentrations were associated with higher velocity, by running linear models adjusted for clustering and with an interaction term between cluster and HMO concentration:

$\begin{array}{ll}\mathit{FPCA}2=\mathit{Intercept}+HMO+\mathit{cluster}+\mathit{cluster}\ast HMO& (1)\end{array}$

where HMO stands for the HMO concentration in mg/L, cluster is a categorical variable 1–4 and the last term accounts for the interaction.

It is expected that the body length, weight, and head circumference are correlated, therefore we report the correlations between the FPCA scores associated to the various growth measures (length, weight, head circumference).

All analyses were performed with R, version 4.0.2. Network analysis was performed with the package igraph (v 1.3.0) and functional PCA was performed using the package fdapace (v 0.5.8).

Results

Description of HMO composition over the lactation period/time

The total measured HMO concentration decreased over the first 3 months of lactation from 9182 (2013) mg/L, at V0, to 7,887 mg/L (1813) at V1 and 6,248 mg/L (1322) (mean (SD)) at V2.

Consistently with previous findings (15) we observed a decrease of almost all HMOs with time of lactation (Supplementary Table 2). A noticeable exception to this decreasing pattern was 3FL, which increased from 704 mg/L (538 mg/L) at V0 to 1,118 mg/L (698 mg/L) at V2.

Mothers can be deemed as secretors/non-secretors based on the presence/absence of α-1-2-linked fucosylated HMOs (2’FL and LNFPI). In secretor mothers, 2’FL was the most abundant HMO, representing up to 58% of the total measured HMO (Supplementary Table 2).

The Shannon diversity at V2 was lower than at V1 (Wilcoxon test, p < 0.01), and lower than at V0 (Wilcoxon test, p < 0.01). This is reflected in the fact that the cumulative proportion of the two most abundant HMOs, 2’FL and 3FL, increased from 29% at V0 to 34% at V1 and increased to 41% at V2. There was no significant difference in diversity between V0 and V1.

Pairwise correlations between the concentrations of individual HMOs are reported in the supplementary materials (Supplementary Figures 1–3): significant correlations ranged between −0.9 (2’FL and LNFP-V at V0) and 0.9 (3FL and LNFP-II at V1). In general, we observed very strong positive correlations between 2’FL, DFNLHa, LNFP-I, between 3FL and LNFP-II, and between DSLNT, LSTb, LNFP-V, LNT at all timepoints. Also, Hex4 HexNAc2 was positively correlated (Spearman ρ = 0.7 at V0) with LNnT, and 6’SL with LSTc (ρ = 0.7 at V0, V1, ρ = 0.8 at V2). 2’FL and LNFP-II were not significantly correlated with LNH, 3’GL, 6’GL, 6’SL, LSTc. LNFP-I was significantly and negatively correlated with 3FL (ρ = −0.8 at V0), LNFP-V (ρ = −0.5 at V0), LNFP-II (ρ = −0.5 at V0), LNnFP-V (ρ = −0.4 at V0), LNFP-III (ρ = −0.3 at V0). In most cases, these correlations and their statistical significance were similar across visits. Other HMOs, like A-Tetrasaccharide or 6’GL, were weakly correlated, or not significantly correlated with the other HMOs.

We calculated the pairwise Aitchison distances between all samples, resulting in a 276 × 276 matrix. Figure 1 compares the distances between samples from the same donors and samples from different donors. As expected, the intra-individual distances between samples were significantly smaller than the inter-individual distances.

FIGURE 1

Figure 1. Comparison of distances between samples. Inter-individual distances were calculated between different mothers, and intra-individual distances between samples from the same mother but at different time points.

We then defined a weighted network with the samples as nodes and edges weighted by the inverse of their distance. Each node corresponds to a subject at a given timepoint. The resulting weighted network was partitioned, producing four clusters of sizes 97(35%), 89 (32%), 62 (22%), 28 (10%), consisting of samples sharing a similar HMO composition.

The average value of each HMO at each timepoint and in each cluster was plotted, as concentrations in mg/L (Figure 2). Overall, it appeared that several HMOs (e.g., 3FL, 2’FL, LNFP-I, LNFP-II) had well separated, non-overlapping, trajectories in the different clusters, while for other HMOs (e.g., 6’SL, LSTc, LSTb) the trajectories were overlapping or even crossing. The concentrations of 6’SL decreased in all clusters, from an average of 468 mg/L at V0 to 119 mg/L at V2 in cluster 1, from 489 mg/L to 147 mg/L in cluster 2, from 535 mg/L to 155 mg/L in cluster 3 and from 588 mg/L to 203 mg/L in cluster 4. In Figure 3, we looked at the temporal consistency of clustering: whether the milk samples remained in the same cluster across visits. All subjects belonging to cluster, 3 or 4 at the first visit, remain in the same cluster afterwards. Some of the subjects starting in cluster 2 switched to cluster 1 afterwards.

FIGURE 2

Figure 2. Each dot represents the average concentration of each HMO, in each cluster and for each time point. Several clustering variables show consistent temporal patterns: for example, LNFP-V is consistently highest in cluster 4, while 3FL is consistently highest in cluster 3 (non-secretors). V0 = 2–5 weeks, V1 = 6 weeks, V2 = 3 months. For other HMOs, like 6’GL, LSTc and DSLNT, the concentration is similar between the clusters. Several HMO concentrations exhibit a non-linear trend, with an increase between the first and second visit and a decrease between the second and third visit.

FIGURE 3

Figure 3. Temporal consistency of clustering: each row corresponds to a mother, the heatmap describes whether the milk samples remained in the same cluster between visits. Grey cells are missing values. All subjects belonging to cluster, 3 or 4 at the first visit, remain in the same cluster afterwards. Some of the subjects starting in cluster 2 (green) switch to cluster 1 (orange) afterwards.

Cluster 1 was characterized as having the highest concentrations of A-TETRA, DFL, LNDFH-I and LNnFP-V and second highest proportions of 3FL, LNFP-II (Figure 2). It also had the lowest levels of MFLNH-III. Cluster 2 was characterized by the highest concentrations of LNnT and LNH, and second highest concentrations of 2’FL, A-TETRA and LNDFH-I. Also, the concentration of 3FL in this cluster was roughly stable across visits. Cluster 3 was characterized by highest concentrations of 3FL, LNFP-II, LNFP-V, LNT, MFLN-III, and lowest concentrations of 2’FL, LNDFH-I, LNFP-I, LNnT. 2’FL was actually absent in this group, which then corresponds to the non-secretor group. LNFP-I was also absent in this cluster. Cluster 4 is characterized by highest concentrations of 2’FL, DFLNHa, LNFP-I, and lowest concentrations of 3FL, LNDF-I, LNFP-II, LNFP-V. Moreover, LNDFH-I and LNFP-II were absent in this cluster. Clusters 3 and 4 had significantly smaller Shannon diversity at all time points compared to clusters 1 and 2, cluster 4 had a significantly lower diversity than cluster 2 at visits V0 and V1, and cluster 2 had lower diversity than cluster 1 at visit V1.

A majority of HMOs had significantly different concentrations between the clusters, at least one timepoint, with the exception of 3’GL, 6’GL, 6’SL, LNnDFH, LSTc (Figure 4).

FIGURE 4

Figure 4. Dark cells correspond to significant difference between the clusters (kruskal test). On the left, HMOs are expressed in relative concentrations (as % of total HMO), on the right, they are expressed in mg/L.

The total measured HMO concentration varied between the clusters and was consistently the lowest in cluster 3 (Figure 5).

FIGURE 5

Figure 5. Sum of measured HMO concentrations, at each visit, split by cluster. Cluster 3, corresponding to non-secretors, had the lowest total concentration of HMOs.

Association with maternal and baseline characteristics

HMO clustering was not associated with mother’s age at recruitment, pre-pregnancy body mass index (ppBMI) or gestational age.

We further tested the impact of BMI status on HMO concentrations, First, total HMO concentration and Shannon diversity index were not different between overweight (OW), obese (OB) and normal weight (NW) groups, at any visit. We then compared the average trajectories for each of the HMOs between the BMI groups regardless of clusters, using functional ANOVA (23): tests were significant only for 6’GL (p = 0.04), LNH (p = 0.04), LSTc (p = 0.04) (Supplementary Figure 4). We also tested if HMO concentrations were different between vaginal delivery and C-section. Specifically, the concentration of LNnT and Hex4 HexNAc2 (2) at 3 months was affected by birth mode and lower in the group of mothers who delivered through C-section (Table 4).

TABLE 4

Table 4. Concentrations of Hex4 HexNAc2 and LNnT at 3 months are lower in the C-section group.

Association with growth

Children in this study showed a normal development for length and weight (see Supplementary Figures 5–8).

Although the differences did not reach statistical significance, children from mothers in cluster 3 (non-secretors) were on average 1 cm shorter than the rest.

Length

When applied to the trajectories of body length, FPCA resulted in 97% of the variance explained by the first two scores, the first score FPCA1 explained 93% of the variance, and was positively correlated with the length z-scores at all time points (Figure 6). The second score (FPCA2) explained 4% of the variance; it was negatively correlated with length at early timepoints, and positively correlated at V7, so positive values of FPCA2 correspond to trajectories having a slow growth before V3 months and an increase at V6 and V7 (Figure 7). We will refer to FPCA1 as ‘general growth’ and to FPCA2 as ‘higher velocity’. See Figure 8 for an illustration of these patterns.

FIGURE 6

Figure 6. FPCA1 score is positively correlated with the z-scores for length (cm). Correlations: 0.54 (V0 = 2–5 wk), 0.70 (V1 = 6 wk), 0.73 (V2 = 3 mo), 0.82 (V3 = 6 mo), 0.81 (V4 = 9 mo), 0.88 (V5 = 12 mo), 0.88 (V6 = 18 mo), 0.78 (V7 = 24 mo). All p-values <0.01.

FIGURE 7

Figure 7. Correlations between FPCA2 scores and length (cm). Correlations: –0.79 (V0), –0.56 (V1), −0,37(V2), −0.05 (V3), −0.13 (V4), −0.12 (V5), 0.3 (V6), 0.35 (V7). Correlations were significant only at V0-V2 and V6-V7.

FIGURE 8

Figure 8. Average trajectories for length. Trajectories were split in 4 groups, based on median FPCA scores. A shared letter indicates a non-significant difference at the corresponding time point. For example, the ‘Low FPCA1l/high FPCA2’ is significantly lower than the ‘Hi FPCA1/low FPCA2’ at V0, V1, V2, V3, but not at V7.

General growth was higher for boys (T test, p < 0.01), indicating that in girls most trajectories of length were below the trajectories for boys. Higher velocity rates were not significantly different between males and females.

Clusters of HMOs at V2 were significantly associated with higher velocity (FPCA2 score), but not with general growth (FPCA1 score), with cluster 4 having significantly higher values than cluster 1 (Figure 9). The same trend was observed at V0 and V1, albeit not significant.

FIGURE 9

Figure 9. The second FPCA score for the child’s length is associated with the HMO clustering. Above the boxplots, are shown the p-values from a post-hoc Dunn pairwise comparison test, with multiplicity correction.

Concentrations of 3’GL, 3FL, 6’GL, DSNLT, LNFP-II, LNFP-III, LNT, LSTb were negatively associated with higher velocity (Supplementary Table 3). For 3FL, this negative association was consistent across the time of lactation (V0, V1, V2).

Weight

For the weight, the first score FPCA1 explained 91% of the variance, and was highly correlated with the weight z-scores (ρ = 0.58 at V0, 0.72 at V1, 078 at V2, 0.84 at V3, 0.9 at V4, 0.9 at V5, 0.93 at V6, 0.9 at V7, all significant). The second score FPCA2 explained 7% of the variance, and Spearman correlation coefficients with z-scores were − 0.5 at V0, −0.57 at V1, −0.59 at V2, −0.53 at V3, −0.37 at V4, all significant (p < 0.01); correlations at V5-V6 were not significant. At V7, there was a significant positive correlation between FPCA2 and the z-score (ρ = 0.21); higher values of FPCA2 correspond to weight trajectories that are below the average before 9 months and above the average afterwards. FPCA1 was positively and significantly correlated with growth velocity at V2, 3 months (ρ = 0.48); at later time points FPCA1 and velocity were not significantly correlated. FPCA2 was significantly positively correlated with velocity at V4, 9 months (ρ = 0.49) and V5, 12 months (ρ = 0.57).

As in the case of length, clusters of HMOs at V2 were significantly associated to higher velocity for weight, as measured by FPCA2 (Kruskal-Wallis, p = 0.02, average values were − 0.7 for clusters 1,2, −0.2 for cluster 3, 2.7 for cluster 4). A post-hoc Dunn comparison test resulted in significantly higher values in cluster 4 compared to the other clusters.

General growth (FPCA1) was higher for boys (T test, p = 0.02) with a mean value of 2.02 kg, versus a mean value of −1.73 for girls. FPCA2 was higher for girls (T test, p = 0.04), with a mean value of 0.33 kg, versus a mean value of −0.57 for boys.

Head circumference

For the head circumference, the first score FPCA1 explained 93% of the variance, and the second score FPCA2 explained 4% of the variance. Again, we can interpret FPCA1 as a measure of general growth and FPCA2 as a measure of higher velocity. FPCA scores for head circumference trajectories were not significantly associated with the HMO clustering.

Discussion

Our analysis aimed to explore the temporal and interindividual variability of the HMOs during the first 3 months of lactation, and to investigate the potential effect on growth. Since several groups of HMOs are highly correlated, multivariate methods may help to understand their interactions, and their cumulative effects. In this exploratory analysis, we proposed a modelling methodology to overcome some of the difficulties that occur when analyzing such data, in particular to address at the same time the non-linearity of growth patterns and their potential association with multi-dimensional, highly correlated HMO concentrations.

We applied a data-driven algorithm to assign mothers to clusters sharing similar HMO profiles. The significance of this clustering is that it is based on the concentrations of a panel of 24 HMOs, and takes into account their mutual correlations. The clustering algorithm produced 4 clusters, distinct from the 4 milk types that are usually defined based on fucosyltransferase (FUT2 and FUT3) polymorphisms. However, the two partitions overlap: for example, our clusters 1 and 2 were composed of samples from mothers belonging to the (FUT2⁺, FUT3⁺) group and they differed in their respective concentrations of several fucosylated HMOs, including 3FL (higher in cluster 1), LNnT (higher in cluster 2), DFL (higher in cluster 1). Therefore, clusters 1 and 2 split the (FUT2⁺, FUT3⁺) group in two subgroups, with distinct levels of several of the most abundant HMOs; this might reflect differences in FUT enzymatic activity with FUT3 having potentially a stronger activity in cluster 1 compared to cluster 2. Cluster 3 corresponded to the (FUT2⁻, FUT3⁺) type. Cluster 4 was analog to the milk group 3 (FUT2⁺, FUT3⁻) with the highest concentrations of 2’FL, DFLNHa, LNFP-I, and lowest concentrations of 3FL, LNDF-I, LNFP-II, LNFP-V.

Previous studies (15) reported that 3FL, LNFP-II and LNnFP-V had their highest concentration in the (FUT2⁻, FUT3⁺) group, similarly to what we observed in our cluster 3. Likewise, 2’FL, LNFP-I and DFLNHa were highest in the (FUT2⁺, FUT3⁻) group, corresponding to our cluster 4. These dynamics can be explained by the substrate and enzyme availability as well as the competition between FUT2 and FUT3 enzymes. Cluster 1 had a lower FUT2 activity compared to FUT3. The opposite is true for Cluster 2. When both enzymes are necessary to synthesize an HMO like in the case of LNDFH-I, more is synthesized in the cluster with a presumably stronger FUT2.

In general, the intra-individual variability was significantly lower compared to inter-individual variability. While subjects were consistently clustered in clusters 3 and 4, independently of the time of lactation, some subjects switched from cluster 2 to 1. This might be explained by a difference in gene regulation and enzyme activity over time, in clusters 1 and 2 with the same Se and Le status, which may follow different patterns of change for each HMO. This means that the overall breastmilk composition, while remaining in the (FUT2⁺, FUT3⁺) group, can vary over the lactation period, with some HMOs decreasing more rapidly than others, changing the ratios between HMOs. We therefore suggest that it may be important to distinguish these subgroups in a statistical analysis, as these different rates of change might have a biological relevance.

When HMOs were expressed as % of total measured HMO, comparing Figure 2 with Supplementary Figure 9, patterns look very similar for several HMOs, including 3FL and LNFP-I, but differences appear for 3’SL and other HMOs. This highlights the fact that the concentration of a specific HMO can increase over time while decreasing as % of total measured HMO.

Only a few studies to date have investigated variations in the HMO composition in relation to maternal age (35), but the results are in general not conclusive. Our clusters of HMO composition were not associated with maternal characteristics such as age and maternal ethnicity. The same lack of association was observed for the levels of 2’FL and LNFP-II and, a fortiori, for the clustering based on the FUT2, FUT3 levels.

A negative association of 3’SL with ppBMI had been reported by Saben et al. (36), who suggested that HMO sialylation may be negatively associated with maternal adiposity. According to Samuel et al. overweight women from a European cohort had significantly higher concentrations of 3′SL, 6′GL and DSLNT (at day 2), 6′SL (at day 17) and LNFP-V (at 3 and 4 months), while lower concentrations of LNnT (at day 2), LNT (at 1 and 3 months) and LNFP-V (at 2 months) compared to normal weight women (p < 0.05 for all). Other studies reported non-significant associations between HMO concentrations and ppBMI (37). We observed in our data that 6’GL, LNH and LSTc were consistently higher in the overweight group. To the best of our knowledge, associations of LNH and LSTc concentrations with ppBMI have not been reported before.

Samuel et al. reported lower concentrations of 2’FL, 3’SL and 6’GL at 2 days, among women delivering through C-section (15). Although we observed the same trend for these HMOs all time points, it was not significant, possibly because this association might be stronger at a very early stage of lactation.

It was proposed that HMO composition might affect child growth by altering the composition of the gut microbiome (38). Our results suggest a possible association of the Lewis negative status with a specific longitudinal growth pattern, as discussed below.

Higher HMO diversity and evenness at 1 month have been associated to lower total and percentage fat mass at 1mo (39). In Lagström et al. (40), the concentration of LNnT (at 3 months) was inversely associated and that of 2’FL (3 months) was directly associated with child length and weight z-scores in a model adjusted for maternal pre-pregnancy BMI, mode of delivery, birthweight z-score, sex and time. In an exploratory study, Sprenger et al. (16) found that FUT2 related alterations of breast milk HMOs composition as assessed through 2’FL concentrations, did not impact growth of breastfed infants during the first 4 months of life. However, they reported a non-statistically significant trend, with males from mothers with low 2’FL milk appearing to have a slightly higher BMI at 1 month, which was not seen any more at 4 months of age when they rather had a smaller BMI and body weight gain. Our approach suggests a complex, dynamic association between growth and HMO composition, illustrated by diverging trajectories starting from 3 months. However, these associations did not address the question of how changes in growth trajectories could be associated with early exposure to higher concentrations of HMO. Therefore, we analyzed the trajectories of length, weight, and head circumference over a period ranging from 2–5 weeks to 24 months, applying a data-driven approach derived from functional data analysis (23, 28, 29). In our study, each trajectory could be accurately described by two numbers (FPCA1, FPCA2), the first describing general growth (higher scores correspond to trajectories above the average at all time points), and the second describing a pattern of increased velocity at 12–18 months. These two scores jointly describe accurately the shape of the individual trajectories and can be therefore taken as a space of coordinates describing the growth trajectories.

Children fed with breast milk of cluster 4 (high in 2’FL, DFNHA, LNFP-I, LNDFH-I, low in 3FL) at 3 months, experienced a slight growth spurt after 12 months, compared to children fed with breast milk of cluster 1 (high in DFSL, LNDF-I). Since cluster 1 belongs to the (FUT2⁺, FUT3⁺) group and cluster 4 belongs to the (FUT2⁺, FUT3⁻) group, this can be interpreted as a possible association of the Lewis negative status with this particular longitudinal growth pattern. This supports the possibility that HMO concentrations should be looked at in combination when assessing associations with development, and that the non-linear patterns in the developmental trajectories must be appropriately modelled.

To the best of our knowledge, this is the first study to investigate the longitudinal association between clusters of HMO composition and child’s physical development. Among the limitations, we acknowledge the modest sample size and the limited age range 0–24 months, limiting the interpretation of the growth patterns to early growth.

Using longitudinal data, we introduced novel methodological approaches, like network analysis of compositional data and functional PCA, to address the multivariate nature of the data. This approach may help to reveal complex longitudinal patterns, and can help future studies in the field.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving humans were approved by Ethics boards at Rhode Island Hospital in Providence, RI, and Pennington Biomedical Research Center, in Baton Rouge, LA. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

FM conducted the data analyses and wrote the manuscript draft. AB, PR, SA, SD, and NS reviewed the manuscript. SA performed the quantification of the HMOs. NS and SD designed the study. All authors contributed to the article and approved the submitted version.

Conflict of interest

FM, PR, AB, SA, and NS were employed by Nestlé Research.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnut.2023.1239349/full#supplementary-material

Footnotes

1. ^ https://www.who.int/tools/child-growth-standards/standards

References

1. Thurl, S, Munzert, M, Boehm, G, Matthews, C, and Stahl, B. Systematic review of the concentrations of oligosaccharides in human milk. Nutr Rev. (2017) 75:920–33. doi: 10.1093/nutrit/nux044

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Plaza-Díaz, J, Fontana, L, and Gil, A. Human milk oligosaccharides and immune system development. Nutrients. (2018) 10:1038. doi: 10.3390/nu10081038

CrossRef Full Text | Google Scholar

3. Sánchez, C, Fente, C, Regal, P, Lamas, A, and Lorenzo, MP. Human milk oligosaccharides (HMOs) and infant microbiota: a scoping review. Foods. (2021) 10:1429. doi: 10.3390/foods10061429

CrossRef Full Text | Google Scholar

4. Marcobal, A, and Sonnenburg, JL. Human milk oligosaccharide consumption by intestinal microbiota. Clin Microbiol Infect. (2012) 18:12–5. doi: 10.1111/j.1469-0691.2012.03863.x

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Walsh, C, Lane, JA, van Sinderen, D, and Hickey, RM. Human milk oligosaccharides: shaping the infant gut microbiota and supporting health. J Funct Foods. (2020) 72:104074. doi: 10.1016/j.jff.2020.104074

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wiciński, M, Sawicka, E, Gębalski, J, Kubiak, K, and Malinowski, B. Human milk oligosaccharides: health benefits, potential applications in infant formulas, and pharmacology. Nutrients. (2020) 12:266. doi: 10.3390/nu12010266

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Karav, S, Casaburi, G, Arslan, A, Kaplan, M, Sucu, B, and Frese, S. N-glycans from human milk glycoproteins are selectively released by an infant gut symbiont in vivo. J Funct Foods. (2019) 61:103485. doi: 10.1016/j.jff.2019.103485

CrossRef Full Text | Google Scholar

8. Karav, S, le Parc, A, Leite Nobrega de Moura Bell, JM, Frese, SA, Kirmiz, N, Block, DE, et al. Oligosaccharides released from milk glycoproteins are selective growth substrates for infant-associated bifidobacteria. Appl Environ Microbiol. (2016) 82:3622–30. doi: 10.1128/AEM.00547-16

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Doherty, AM, Lodge, CJ, Dharmage, SC, Dai, X, Bode, L, and Lowe, AJ. Human milk oligosaccharides and associations with immune-mediated disease and infection in childhood: a systematic review. Front Pediatr. (2018) 6:91. doi: 10.3389/fped.2018.00091

CrossRef Full Text | Google Scholar

10. Moore, RE, Xu, LL, and Townsend, SD. Prospecting human Milk oligosaccharides as a defense against viral infections. ACS Infect Dis. (2021) 7:254–63. doi: 10.1021/acsinfecdis.0c00807

CrossRef Full Text | Google Scholar

11. Kaplan, M, Şahutoğlu, AS, Sarıtaş, S, Duman, H, Arslan, A, Pekdemir, B, et al. Role of milk glycome in prevention, treatment, and recovery of COVID-19. Front Nutr. (2022) 9:1033779. doi: 10.3389/fnut.2022.1033779

CrossRef Full Text | Google Scholar

12. Thai, JD, and Gregory, KE. Bioactive factors in human breast milk attenuate intestinal inflammation during early life. Nutrients. (2020) 12:581. doi: 10.3390/nu12020581

CrossRef Full Text | Google Scholar

13. Cheng, YJ, and Yeung, CY. Recent advance in infant nutrition: human milk oligosaccharides. Pediatr Neonatol. (2021) 62:347–53. doi: 10.1016/j.pedneo.2020.12.013

CrossRef Full Text | Google Scholar

14. Samuel, T, Hartweg, M, Lebumfacil, JD, Buluran, KB, Lawenko, RB, Estorninos, EM, et al. Dynamics of human milk oligosaccharides in early lactation and relation with growth and appetitive traits of Filipino breastfed infants. Sci Rep. (2022) 12:17304. doi: 10.1038/s41598-022-22244-7

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Samuel, TM, Binia, A, de Castro, CA, Thakkar, SK, Billeaud, C, Agosti, M, et al. Impact of maternal characteristics on human milk oligosaccharide composition over the first 4 months of lactation in a cohort of healthy European mothers. Sci Rep. (2019) 9:11767. doi: 10.1038/s41598-019-48337-4

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Sprenger, N, Lee, LY, de Castro, CA, Steenhout, P, and Thakkar, SK. Longitudinal change of selected human milk oligosaccharides and association to infants’ growth, an observatory, single center, longitudinal cohort study. PLoS One. (2017) 12:e0171814. doi: 10.1371/journal.pone.0171814

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Berger, PK, Plows, JF, Jones, RB, Alderete, TL, Yonemitsu, C, Poulsen, M, et al. Human milk oligosaccharide 2′-fucosyllactose links feedings at 1 month to cognitive development at 24 months in infants of normal and overweight mothers. PLoS One. (2020) 15:e0228323. doi: 10.1371/journal.pone.0228323

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Lefebvre, G, Shevlyakova, M, Charpagne, A, Marquis, J, Vogel, M, Kirsten, T, et al. Time of lactation and maternal Fucosyltransferase genetic polymorphisms determine the variability in human Milk oligosaccharides. Front Nutr. (2020) 7:574459. doi: 10.3389/fnut.2020.574459

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Van Leeuwen, SS . Challenges and pitfalls in human milk oligosaccharide analysis. Nutrients. (2019) 11:2684. doi: 10.3390/nu11112684

CrossRef Full Text | Google Scholar

20. Hanson, LÅ, Korotkova, M, Lundin, S, Håversen, L, Silfverdal, SA, Mattsby-Baltzer, I, et al. The transfer of immunity from mother to child. Ann N Y Acad Sci. (2003) 987:199–206. doi: 10.1111/j.1749-6632.2003.tb06049.x

CrossRef Full Text | Google Scholar

21. Goehring, KC, Kennedy, AD, Prieto, PA, and Buck, RH. Direct evidence for the presence of human milk oligosaccharides in the circulation of breastfed infants. PLoS One. (2014) 9:e101692. doi: 10.1371/journal.pone.0101692

CrossRef Full Text | Google Scholar

22. Han, K, Hadjipantelis, PZ, Wang, JL, Kramer, MS, Yang, S, Martin, RM, et al. Functional principal component analysis for identifying multivariate patterns and archetypes of growth, and their association with long-term cognitive development. PLoS One. (2018) 13:e0207073. doi: 10.1371/journal.pone.0207073

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Chen, K, Zhang, X, Petersen, A, and Müller, HG. Quantifying infinite-dimensional data: functional data analysis in action. Stat Biosci. (2017) 9:582–604. doi: 10.1007/s12561-015-9137-5

CrossRef Full Text | Google Scholar

24. Zhang, W, and Wei, Y. Regression based principal component analysis for sparse functional data with applications to screening growth paths. Ann Appl Stat. (2015) 9:597–620. doi: 10.1214/15-AOAS811

CrossRef Full Text | Google Scholar

25. Sheppard, T, Tamblyn, R, Abrahamowicz, M, Lunt, M, Sperrin, M, and Dixon, WG. A comparison of methods for estimating the temporal change in a continuous variable: example of HbA1c in patients with diabetes. Pharmacoepidemiol Drug Saf. (2017) 26:1474–82. doi: 10.1002/pds.4273

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yao, F, Müller, HG, and Wang, JL. Functional data analysis for sparse longitudinal data. J Am Stat Assoc. (2005) 100:577–90. doi: 10.1198/016214504000001745

CrossRef Full Text | Google Scholar

27. Ullah, S, and Finch, CF. Applications of functional data analysis: a systematic review. BMC Med Res Methodol. (2013) 13:43. doi: 10.1186/1471-2288-13-43

CrossRef Full Text | Google Scholar

28. Viviani, R, Grön, G, and Spitzer, M. Functional principal component analysis of fMRI data. Hum Brain Mapp. (2005) 24:109–29. doi: 10.1002/hbm.20074

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Wu, PS, and Müller, HG. Functional embedding for the classification of gene expression profiles. Bioinformatics. (2010) 26:509–17. doi: 10.1093/bioinformatics/btp711

CrossRef Full Text | Google Scholar

30. Schneider, N, Bruchhage, MMK, O'Neill, BV, Hartweg, M, Tanguy, J, Steiner, P, et al. A nutrient formulation affects developmental myelination in term infants: a randomized clinical trial. Front Nutr. (2022) 9:823893. doi: 10.3389/fnut.2022.823893

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Austin, S, and Bénet, T. Quantitative determination of non-lactose milk oligosaccharides. Anal Chim Acta. (2018) 1010:86–96. doi: 10.1016/j.aca.2017.12.036

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Dey, AK, Tian, Y, and Gel, YR. Community detection in complex networks: from statistical foundations to data science applications. Wiley Interdiscip Rev Comput Stat. (2022) 14. doi: 10.1002/wics.1566

CrossRef Full Text | Google Scholar

33. Greenacre, M . Compositional data analysis in practice. New York: Chapman and Hall/CRC (2018).

Google Scholar

34. Newman, MEJ . Finding community structure in networks using the eigenvectors of matrices. Phys Rev E Stat Nonlinear Soft Matter Phys. (2006) 74:036104. doi: 10.1103/PhysRevE.74.036104

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Han, SM, JGB, D, Binia, A, Sprenger, N, Vickers, MH, and Cutfield, WS. Maternal and infant factors influencing human Milk oligosaccharide composition: beyond maternal genetics. J Nutr. (2021) 151:1383–93. doi: 10.1093/jn/nxab028

CrossRef Full Text | Google Scholar

36. Saben, JL, Sims, CR, Abraham, A, Bode, L, and Andres, A. Human milk oligosaccharide concentrations and infant intakes are associated with maternal overweight and obesity and predict infant growth. Nutrients. (2021) 13:446. doi: 10.3390/nu13020446

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Azad, MB, Robertson, B, Atakora, F, Becker, AB, Subbarao, P, Moraes, TJ, et al. Human milk oligosaccharide concentrations are associated with multiple fixed and modifiable maternal characteristics, environmental factors, and feeding practices. J. Nutr. (2018) 148:1733–42. doi: 10.1093/jn/nxy175

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Charbonneau, MR, O’Donnell, D, Blanton, LV, Totten, SM, Davis, JCC, Barratt, MJ, et al. Sialylated Milk oligosaccharides promote microbiota-dependent growth in models of infant undernutrition. Cells. (2016) 164:859–71. doi: 10.1016/j.cell.2016.01.024

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Alderete, TL, Autran, C, Brekke, BE, Knight, R, Bode, L, Goran, MI, et al. Associations between human milk oligosaccharides and infant body composition in the first 6 mo of life. Am J Clin Nutr. (2015) 102:1381–8. doi: 10.3945/ajcn.115.115451

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Lagström, H, Rautava, S, Ollila, H, Kaljonen, A, Turta, O, Mäkelä, J, et al. Associations between human milk oligosaccharides and growth in infancy and early childhood. Am J Clin Nutr. (2020) 111:769–78. doi: 10.1093/ajcn/nqaa010

PubMed Abstract | CrossRef Full Text | Google Scholar