Comprehensive Targeted and Quantitative Profiling of the Human Milk Metabolome: Impact of Delivery Mode, Breastfeeding Practices, and Maternal Diet.

J.C.‐L., R.C.‐R., and C.L. contributed equally to this work and share first authorship. The authors’ responsibilities were as follows—M.C.C. and T.B.: designed the study; C.M.‐C.: responsible of clinical study, recruitment, and collection of the milk samples; C.L. and T.B. analyzed milk samples for metabolite composition; M.S.R.; sample management and metadata processing; S.G.: responsible for the dietary information and analysis; R.C.; responsible for the biocomputational analysis and data report; J.C.L.: integrated data and drafted the first draft; and all authors contributed to the data interpretation, and also all have read and approved the final manuscript.

From 125 HM samples, a final total number of 123 HM samples with available maternal‐infant clinical and maternal dietary data availability were analyzed ( Figure S1 , Supporting Information). The Biocrates MxP Quant 500 kit platform quantified 432 metabolites in the HM samples out of the 630 that is capable of identifying. These metabolites were categorized into 20 analyte classes including acyl‐carnitines, alkaloids, amine oxides, amino acids, amino acid‐related compounds, bile acids, biogenic amines, carboxylic acids, ceramides, cholesterol esters, cresols, diacylglycerols, fatty acids, glycerol‐phospholipids, glycosyl‐ceramides, indoles derivatives, sphingolipids, sugars, triacylglycerols, and vitamins and cofactors ( Table   1 ). Based on the analytes quantified by the platform, the most abundant group was triacylglycerols (median 12 660 µM), followed by amino acids (2087 µM), fatty acids (1983 µM), and sugars (1443 µM). Conversely, those present in the lowest concentrations (less than 1 µM) included bile acids, indoles derivatives, and alkaloids (Table  1 ). Out of the 630 compounds that the method is able to quantify, 198 were not detected. The number of unidentified compounds were distributed among the analyte classes as follows: acyl‐carnitines (29/40), alkaloids (0/1), amine oxides (0/1), amino acids (0/20), amino acid‐related compounds (14/30), bile acids (10/14), biogenic amines (4/9), carboxylic acids (4/7), ceramides (10/28), cholesterol esters (15/22), cresols (0/1), diacylglycerols (24/44), fatty acids (3/12), glycerol‐phospholipids (28/90), glycosyl‐ceramides (16/34), indoles derivatives (3/4), sphingolipids (0/15), sugars (0/1), triacylglycerols (25/242) and vitamins and cofactors (0/1). In addition, the analyte classes from which no metabolite was detected were dihydroceramides (eight available compounds), or nucleobases related (1), hormones (4). Quantification of metabolite groups in the whole cohort and in the three clusters of HM metabolome profiles (median and IQR) and FDR p ‐value. Clustering analysis based on metabolic profiles defined three distinct clusters ( Figure   1 A,B ), revealing statistically significant differences in metabolite groups among them. The main difference between clusters (Table  1 ) was observed in triacylglycerols (FDR p ‐value < 0.00001), with concentrations of 21 814, 11 419, and 4258 µM in clusters 1, 2, and 3, respectively. Additionally, other lipid‐related metabolite groups exhibited a similar decreasing trend from Cluster 1 to 3, including cholesterol esters (FDR p ‐value <0.0001), ceramides (FDR p ‐value <0.0001), diacyl‐glycerols (FDR p ‐value = 0.002), and glycerol‐phospholipids (FDR p ‐value < 0.0001). Conversely, amino acids were significantly lower in cluster 1 compared to the other clusters (FDR p ‐value of cluster 1 vs 2 = 0.035; cluster 1 vs 3 = 0.042). HM metabolomic profile and differential metabolite groups. A) k‐means clustering with the algorithm PAM. The optimal centers of k‐means according to the method were established as 3 groups. B) Optimal number of clusters based on the Elbow method. Further examination of specific metabolites ( n = 432) revealed significant differences among the three clusters, with 268 metabolites showing statistical significance (FDR p ‐value < 0.05) (Table S1 , Supporting Information). No differences were found in maternal‐infant clinical characteristics among the three clusters ( Table   2 ) except for a lower prevalence of C‐section births in cluster 3 compared to clusters 1 and 2 (FDR p ‐value = 0.035). We then examined potential differences at the metabolite level between the two birth modes ( Figure   2 A ). We observed higher levels of the essential fatty acid docosahexaenoic acid (DHA, p ‐= 0.001, FDR p ‐value > 0.05) in human milk from women who had vaginal births compared to those who underwent C‐sections (Figure  2A ), although the observation was not significant after adjustment due the high variability in the concentrations, we highlight this result due the relevance of DHA for infant's development. Regarding breastfeeding practices (Figure  2B ), hexoses (sugars) showed significantly higher levels in milk from women following exclusive breastfeeding compared to mixed feeding (FDR p ‐value = 0.0098). Clinical and anthropometrical characteristics. Pre‐gestational weight (kg) Mean (SD) Gestational weight gain (kg) Mean (SD) Pre‐gestational BMI (kg m −2 ) Mean (SD) Gestation time (weeks) Mean (SD) Mode of delivery (vaginal, %) Infant gender (female/male) Mode of feeding (mixed/breastfeeding) Infant WFA at birth Median (IQR) Infant WFA at 1 month Median (IQR) Infant WFA at 1 year Median (IQR) Infant LFA at birth Median (IQR) Infant LFA at 1 month Median (IQR) Infant LFA at 1 year Median (IQR) Accumulated infant infections at 1 year ( n ) 0 infections 1 infection >1 infection 89/123 19/123 15/123 31/44 8/44 5/44 39/51 5/51 7/51 20/28 6/28 2/28 Associations between HM metabolites and perinatal factors. Top 25 metabolites associated to (A) vaginal (blue) and C‐section (pink) deliveries and also. B) Breastfeeding practices: exclusive (blue) and mixed breastfeeding (pink). Spearman's rank correlations using pattern search function in MetaboAnalyst 5.0 between factors were tested and p ‐values indicate statistically significant difference between the most contributing metabolite. p ‐values were also adjusted for multiple comparisons using False Discovery Rate (FDR) based on Benjamini–Hochberg. Concerning metabolic Ratios or/ and Sums of HM metabolites obtained using MetaboINDICATOR, no significant differences were identified according to mode of delivery, feeding type, and maternal dietary cluster, however, differences between metabolite profile clusters were found ( Figure   3 A–F ). Associations between sums and ratios of HM metabolites according to metabolite cluster. Top 25 ratios HM associated to metabolome cluster 1 versus cluster 2 (A), to cluster 1 versus cluster 3 (B); and cluster 2 versus cluster 3 (C). Top 25 Sums HM metabolites associated to metabolome cluster 1 versus cluster 2 (D); cluster 1 versus cluster 3 (E); and cluster 2 versus cluster 3 (F). Spearman's rank correlations using pattern search function in MetaboAnalyst 5.0 between factors were tested and p ‐values indicate statistically significant difference between the most contributing metabolite. p ‐values were also adjusted for multiple comparisons using False Discovery Rate (FDR) based on Benjamini–Hochberg. In terms of metabolic activity indicators ( Figure   4 ), same observations were found. No significant differences between exclusive breastfeeding and mixed feeding (Figure  4A ) or between vaginal and C‐section deliveries (Figure  4B ) were identified. However, several significant differences were detected when comparing metabolome clusters. Notably, comparisons between cluster 1 and cluster 2 (Figure  4C ) revealed significant enrichment ratios (FDR p ‐values) in nicotinate and nicotinamide metabolism (FDR p  = 2.3 × 10 −11 ), nitrogen metabolism (FDR p ‐value = 4.8 × 10 −11 ), and d ‐glutamine and d ‐glutamate metabolism (FDR p ‐value = 4.8 × 10 −11 ). Comparisons between cluster 1 and cluster 3 (Figure  4D ) demonstrated significant differences in fatty acid biosynthesis (FDR p ‐value = 9.4 × 10 −10 ), primary bile acid biosynthesis (FDR p ‐value = 2.3 × 10 −9 ), and taurine and hypotaurine metabolism (FDR p  = 9.4 × 10 −10 ). Likewise, comparisons between cluster 2 and cluster 3 (Figure  4E ) indicated the highest disparities in fatty acid biosynthesis (FDR p ‐value = 6.0 × 10 −22 ), followed by selenocompound metabolism (FDR p ‐value = 1.5 × 10 −19 ) and histidine metabolism (FDR p ‐value = 5.2 × 10 −21 ). Differences in metabolism indicators according to perinatal factors and metabolite cluster. A) Breastfeeding practices (exclusive breast feeding vs mixed feeding); B) mode of delivery (vaginal vs C‐section); C) metabolome cluster 1 versus cluster 2; D) metabolome cluster 1 versus cluster 3; and e) metabolome cluster 2 versus cluster 3. The Enrichment Analysis using the “ CalculateGlobalTestScore ” function in MetaboAnalyst 5.0 between factors were tested with the Q ‐stat. p ‐Values were also adjusted for multiple comparisons using False Discovery Rate (FDR) based on Benjamini–Hochberg. Using MetaboAnalystR Package approach to analyze the metabolic pathway analysis, we confirmed the results identified with the metabolic activity indicators obtained with biocrates’ MetaboINDICATOR tool (Figure S2 , Supporting Information) where no significant differences according breastfeeding type and mode of delivery, but several significant differences were observed between metabolite clusters (Table S3 , Supporting Information). We next examined whether the HM metabolome was associated to maternal diet. Dietary components significantly associated with HM metabolome clusters are depicted in Figure   5 . Cluster 3 was characterized by significantly higher amounts of most macronutrients, minerals, and vitamins. Cluster 2 exhibited higher amounts of bioactive compounds such as catechin, lycopene, caffeic acid, or quercetin. In contrast, cluster 1 was distinguished by lower amounts of total and digestible carbohydrates and vitamin B12, along with the highest intake of trans fatty acids. Specific Dietary components are associated to the specific HM metabolome cluster. LDA effect size (LEfSe) analysis was used to identify the specific dietary compounds: a) macronutrients and dietary fiber; b) minerals and vitamins; and c) bioactive compounds associated to HM metabolome cluster. Blue color indicates the metabolome cluster with significantly higher value compared to the other two; red color indicates the cluster with significantly lower value compared to the other two clusters. Non‐parametric factorial Kruskal–Wallis sum‐rank p < 0.05 was considered significant. The potential impact of the diet on the HM metabolome was also analyzed independently of the metabolome cluster, against the dietary patterns of the subjects previously reported [ 16 ] dietary pattern 1 (high vegetal protein and high fiber intake) and dietary pattern 2 (high animal protein and high lipids intake). We observed significant associations between dietary patterns and three metabolite groups. Specifically, HM from women adhering to dietary pattern 2 (high animal protein and high lipid intake) had higher levels of amino acids ( p = 0.024, FDR p = 0.185), sphingolipids ( p = 0.009, FDR p = 0.105), and glycosyl‐ceramides ( p = 0.005, FDR p = 0.105) compared to women in dietary pattern 1 (high vegetal protein and fiber intake) ( Figure   6 ). Associations between HM metabolites and maternal dietary patterns. Spearman correlation using pattern search function in MetaboAnalyst 5.0 was performed and p ‐values indicate statistically significant differences between dietary patterns. p ‐Values were also adjusted for multiple comparisons using False Discovery Rate (FDR) based on Benjamini–Hochberg. Pattern 1: high vegetal protein and high fiber intake; Pattern 2: high animal protein and high lipids intake.

This study represents the first comprehensive exploration of the targeted and quantitative HM metabolome using a multi‐nature compound identification technique capable of quantifying over 400 metabolites simultaneously. Our findings revealed different metabolomic patterns, with certain factors such as birth mode potentially explaining variations in HM composition. The relevance of this study lies in its comprehensive approach, which encompasses the analysis of various metabolite families that had not been previously examined as a whole, but as individual smaller metabolite groups. Previous studies on the metabolic profile of human milk have typically focused on particular groups of metabolites due to specific research priorities and the need for the multiple methodological tools required to assess all of them. By assessing these components collectively, our study elucidates specific metabolome patterns in HM across a wide range of components. Within our sample population, statistical analysis revealed three distinct clusters or groups in the HM metabolomic profile. Triacylglycerols were the metabolite group showing the greatest differences among the three clusters, suggesting that this type of lipid, which represents 98% of the HM total lipids [ 25 ] are highly variable. This finding is consistent with existing literature indicating that lipid content and structure are the most variable of the HM components, depending on factors such as geographic location, diet, and other environmental conditions. [ 26 ] On the other hand, fatty acids, carboxylic acids, and indole derivatives seemed to be evenly distributed among the subjects in terms of concentration. Assessing dietary elements in relation to HM metabolome clusters revealed significant correlations. Notable differences were observed in intake of some macronutrients, minerals, and bioactive compounds between the three metabolome clusters. It is well documented that HM composition is resilient to dietary factors. However, in the present study we did not only assess nutrient‐related metabolites in HM but also a wide range of metabolites, defining three metabolomic clusters. In this sense, our study contributes to the growing body of evidence on the effects of diet on HM composition. Previous studies focused on elucidating the impact of diet on the HM nutrient (not metabolome) composition only reported indirect associations, and only specific correlations (e.g., DHA or Vitamin C). [ 27 ] Moving to an extreme example, a systematic review concluded that all non‐vegetarian, vegetarian, and vegan mothers produce breast milk of comparable nutritional value. [ 28 ] When addressing the effect of diet from a different perspective, we also found some coherent results. Defining the diet into two patterns, i.e., high‐animal protein and lipids on one hand, and high‐vegetal protein and fiber on the other, we observed that three metabolite groups tended to correlate with the first of the dietary patterns: glyceroceramides (FDR p ‐value = 0.105), sphingolipids (FDR p ‐value = 0.105), and aminoacids (FDR p ‐value = 0.185), which is coherent with the high proportion of dietary lipids and protein compounds. No significant impact of mode of birth on the global metabolite profile was observed, although we found that C‐section was associated to lower DHA concentrations than in those levels observed in human milk samples from women who had vaginal delivery in agreement with other reports. [ 29 ] In addition, the other maternal factor showing a statistically significant difference in HM composition was type of infant feeding: hexoses (or sugars) were higher in women giving exclusive breastfeeding than in those on mixed feeding with formula (FDR p ‐value = 0.0098). This result could be supported by previous reports showing that HM lactose was positively correlated to the volume of intake [ 30 ] and those infants on the exclusive breastfeeding regime are supposed to have higher intake than in the case of mixed feeding. [ 31 ] While methodologies allowing for multi‐compound determination have been available for years, they were not applied to HM analysis until recently. Specifically, the Biocrates Quant 500 platform used in this study enables the quantification of a large series of different compounds covering a broad range of pathways. In a recent study in which the same Biocrates Quant 500 method was applied to serum samples up to 600 compounds were identified and classified in groups in a similar way than in this study, with the focus on identifying pollutants and their possible association with endometriosis. [ 32 ] Furthermore, the technique was first applied in HM to assess its effectiveness on this biological fluid, by using a pooled sample, concluding it was successful in identifying 400 metabolites. [ 33 ] In this field, another study focused on assessing HM metabolome related to maternal nutritional status (obesity), and its impact on infants’ body composition and weight gain after birth. [ 34 ] The study relied on untargeted liquid chromatography–gas chromatography–mass spectrometry and identified 275 metabolites from 35 breast milk samples providing a semiquantitative analysis. Our study with the Biocrates Quant 500 platform contributes to research by testing its application in 123 subjects, and by increasing the number of detected and quantified components to 432, thanks to the currently available kit that includes more specific standards. The present cross‐sectional study enabled the identification of patterns in the HM metabolome and their associations with maternal and perinatal characteristics. These results encourage the application of the newly used analytical tools in other research contexts, especially in interventional studies in which the assessment of the effect of specific diets or other variables on the HM metabolome could be more comprehensively explored. Also, this methodology would be relevant to compare HM composition among different groups (preterm vs term infants; geographical locations; etc.), which up to now have been focused in one component family or in a small range of components. [ 2 , 35 , 36 , 37 ] The application of this analytical approach to knowledge generation will contribute to the enrichment of the actual evidence and knowledge on HM metabolites, while simultaneously promoting the HM research area. Despite the novel data, our study has some potential limitations. First, the cross‐sectional nature of our study may not fully reflect the HM metabolome, which may change during lactation. Therefore, longitudinal studies are needed to gain a more comprehensive understanding. In this sense, HM composition is known to vary depending on perinatal factors including preterm birth, maternal diet (fat intake, caloric intake), and other maternal health conditions (caquexia, obesity, and others). Our study did not consider these factors and they could be of relevance. In addition, HM composition is known to vary in macronutrient composition throughout the day [ 38 ] so other metabolites could be also modified, although this possibility was not considered in our study. Nevertheless, despite the aforementioned limitations, our study represents the first to quantify a collection of more than 400 metabolites in over 120 HM samples. Associations were established between the composition and characteristics of the subjects.

In conclusion, the characterization of the HM metabolome including components from different families, demonstrated that the main differences in the metabolomic pattern are influenced by the concentration of triacylglycerols, diacylglycerols, cholesterol esters, ceramides, and amino acids. Maternal factors, such as mode of delivery or maternal diet partially explain differences in the metabolome profile and, breastfeeding practices (exclusive vs mixed) resulted in a significant difference in the content of hexoses (sugars). The progress in the characterization of HM targeted and quantitative metabolites and metabolite sums and ratios achieved in this study provides a solid foundation for further applications in interventional prospective studies.

In this cross‐sectional study, HM samples from 125 healthy women at 1 month postpartum from the MAMI birth cohort [ 13 ] were included in the study and subjected to targeted metabolite quantification. None of the participants had any diagnosed pathologies or were undergoing drug or pro‐ and prebiotic treatments, except for some who had received antibiotics during pregnancy or at delivery. The study was approved by the Ethical Committee of the Hospital Clínico Universitario Valencia, Spain and by CSIC Ethics Committee. Written informed consent was obtained from all the participants prior to enrollment in the study. The study was registered on the ClinicalTrials.gov platform (registration number NCT03552939 ). Anthropometrical and clinical parameters of mother–infant pairs were available, including maternal age, gestational age, pre‐gestational body mass index (BMI), weight gain during pregnancy, mode of delivery (vaginal or C‐section), antibiotic exposure, and type of feeding (exclusive breastfeeding or mixed feeding). Maternal and neonatal anthropometric data were also collected in the subjects and their children, and weight‐for‐age (WFA) and length for age (LFA) z‐scores were calculated according to the WHO references ( www.who.int/childgrowth/software/en/ ). Maternal dietary information covering the pregnancy period was collected and characterized in terms of energy, nutrients, and other specific compounds as previously described. [ 14 ] In brief, dietary records were collected by a nutritionist using a 140‐item validated food frequency questionnaire (FFQ). FFQ information was analyzed for the energy and daily intakes of macro‐ and micro‐nutrients by using the nutrient Food Composition Tables developed by the Centro de Enseñanza Superior de Nutrición Humana y Dietética (CESNID). The intake of specific dietary fiber, as soluble and insoluble fiber types, was completed from the Marlett food composition tables. [ 15 ] In addition, maternal dietary patterns, previously established by means of Jensen–Shannon distance and partitioning around medoid clustering [ 16 ] were also used to explore correlations with the HM metabolites. HM samples from healthy women were collected following a previously established protocol. [ 13 ] Women were trained and received detailed instructions to self‐collect HM samples. Morning collection and the use of sterile pump were recommended to standardize the protocol. Briefly, participants cleaned their nipples and the surrounding breast skin before collection. The initial drops of milk were discarded to minimize skin contamination, and the remaining milk was collected into sterile tubes using a sterile pump (approximately 10 mL were collected) and transferred to a biobank, aliquoted, and stored at −80 °C for subsequent analysis. An aliquot of 10 µL of whole milk was used for the targeted quantitative metabolomic analysis. Targeted metabolomic analysis was based on ultra‐high‐performance liquid chromatography tandem mass spectrometry (UHPLC‐MS/MS) using Shimadzu Nexera chromatography platform (Shimadzu Corporation, Kyoto, Japan) coupled to Sciex QTrap 5500 mass spectrometer (AB Sciex LLC, Framingham, MA, USA). [ 17 ] The study applied the Biocrates MxP Quant 500 targeted kit (Biocrates Life Sciences AG, Innsbruck, Austria), potentially quantitating 106 small molecules and free fatty acids in chromatography mode and 524 complex lipids in positive flow injection mode (FIA‐MS/MS). [ 17 ] The study followed the manufacturer's instructions to perform the HM analysis. In brief, 10 µL of whole HM sample was transferred onto a kit plate with pre‐injected internal standards and dried down. The plate also included blanks, quality control samples (three levels), and 7‐point calibrator standards. The rest of the assay included derivatization with 5% phenylisothiocyanate in pyridine, ethanol, and water (1:1:1), and subsequent extraction with 5 mM ammonium acetate in methanol. Chromatography was done with 0.2% formic acid in acetonitrile (organic mobile phase) and 0.2% formic acid in water (inorganic mobile phase). Flow‐injection analysis was performed with methanol and Biocrates MxP Quant 500 additive of undisclosed composition. In addition to the metabolites, the study used biocrates’ MetaboINDICATOR tool to obtain metabolism indicators ( n = 232 potential metabolites sums and ratios) to extract biological significance from targeted metabolomic and providing links to biochemical pathways. Peak areas and concentrations were obtained as described elsewhere. [ 17 ] Mass spectrometry signal was acquired in Sciex Analyst v1.6.24 (AB Sciex LLC, Framingham, MA, USA) and chromatographic peaks were identified and integrated in Biocrates MetIDQ Oxygen‐DB110‐3005 (Biocrates Life Sciences AG, Innsbruck, Austria). Plates were normalized (per metabolite) by scaling through median normalization: for a given metabolite, values of reference samples in each plate were scaled by such a factor so that their median was equivalent to the median of values of all reference samples before normalization. LOD was calculated as 2× median signal in blank samples. Metabolites >50% of their values <LOD across all subject groups were excluded from the analysis. Values <LOD were not adjusted, as they represented the best estimate of the true values. However, values that were entirely zero were interpolated as half of the minimum non‐zero value for that specific metabolite to avoid absolute zeros. Metabolites Raw data were converted to mzXML format using ProteoWizard software. Then, a series of data preprocessing methods was implemented by an in‐house R package (based on XCMS), including peak detection, extraction, alignment, and integration. Data were summarized by median and 1 st and 3 rd quartile ( Q ) due to the non‐normal distribution, and shown in the tables of descriptive results (clinical data and metabolite quantification into metabolite groups). In order to identify metabolites that were similar, k‐means clustering was performed to divide this set of metabolites into a set of K groups, also in order that this clustering was not so sensitive to outliers (affecting the cluster assignment), the “Partitioning Around Medoids” (PAM) algorithm was used, which provided a more robust algorithm. All this was performed with the “factoextra” R package, [ 18 ] which determined the optimal number of clusters and k‐means clustering of the data. The optimal centers of k‐means according to the method [ 19 ] were established as three groups. This clusterization abundance was confirmed by executing “discretize” function of “ arules ” R package. [ 20 ] To compare differences between the metabolome clusters, the study obtained p ‐values with the False discovery ratio method to adjust for multiple comparisons. Differences in clinical parameters between metabolite clusters were analyzed with multi‐variable ANOVA, and Spearman correlation to assess differences in individual metabolites and metabolite groups between the three metabolite clusters, Kruskal–Wallis non‐parametric one‐way ANOVA (including multiple comparisons between the clusters) were applied on the log‐transformed raw data; then, the p ‐values were adjusted for multiple comparisons using False Discovery Rate (FDR) based on the method of Benjamini–Hochberg. [ 21 ] The dietary components that differed significantly between metabolite groups were identified using linear discriminant effect size analysis (LEfSe). This analysis integrated statistical significance (non‐parametric factorial Kruskal–Wallis, p < 0.05) with effect size estimation (linear discriminant analysis). Using this method, each dietary constituent was identified as being significantly more or less abundant in each of the three metabolite groups. The analysis of the indicators of metabolic activity was performed with MetaboAnalystR Package [ 22 ] with statistical analyses that took the covariates. The metabolite table was first filtered using Interquartile range (IQR) to identify and eliminate variables that were probably not useful when modeling the data. [ 23 ] Following a subsequent normalization step, was carried out in three steps: 1) the normalization of the sample was for adjusting the systematic differences between the samples; 2) the data transformation applied a mathematical transformation to the individual values themselves. Performing a transformation to Log. 3) The data scale adjusts each variable/characteristic by a scale factor calculated based on the dispersion of the variable. This process was performed with MetaboAnalystR Package [ 22 ] to know the indicators of metabolic activity. Spearman correlation analysis was performed for a given feature. Differences in the metabolite indicators between women’ characteristics (mode of delivery, breastfeeding practices, and metabolome cluster) were analyzed by the metabolite set enrichment analysis (MSEA) containing 80 metabolite sets based on KEGG human metabolic pathways [ 24 ] from the MetaboAnalystR Package [ 22 ] with the same normalization steps as described above in the methods section, in order to relate the metabolites functionally and whether they were significantly enriched. Enrichment tests were based on the well‐established global test to test associations between metabolite sets and the outcome. The algorithm used a generalized linear model to compute a “Q‐stat” for each metabolite set. The Q‐stat was calculated as the average of the Q values calculated for each single metabolite; while the Q value was the squared covariance between the metabolite and the outcome. In addition, the study carried out the metabolic pathway analysis by use of MetaboAnalystR Package (integrating pathway enrichment analysis and pathway topology analysis) and visual exploration for >120 species. For this analysis, the study used the same methodology protocol than enrichment analysis.

Human milk (HM) is a complex fluid comprising not only essential nutrients but also a wide array of bioactive compounds. [ 1 , 2 ] Given the established impact of these compounds on infant growth and development, various methodologies for their identification and quantification have been developed and are now widely utilized. While straightforward analytical procedures are of commonly used for determining immunoglobulins or microbiota, a broader spectrum of methodologies is required for quantifying nutrients and other bioactive compounds. [ 3 , 4 ] Oligosaccharides are among the HM bioactive compounds that attract the most attention due to their role in microbiota development as prebiotics. [ 5 ] Typically, their determination involves reverse‐phase HPLC after derivatization. Lipids, another important constituent, have also been extensively studied given that the content and composition of lipids are particularly sensitive to dietary and environmental influences. [ 6 ] The analysis of HM lipid species has relied on methods including gas chromatography, involving prior sample treatment (extraction, derivatization, transesterification, etc.), and nuclear magnetic resonance. [ 7 , 8 , 9 , 10 ] Amino acids and related components are also of interest due to their pivotal role in infant growth, with techniques like ultra‐performance liquid chromatography tandem mass spectrometry and ion‐paired liquid chromatography being commonly employed for their quantification. [ 2 , 11 , 12 ] However, several other HM components have roles that are yet to be fully established. Advances in understanding the role of specific constituents of human milk are critical to understanding their impact on future infant development. This understanding would greatly benefit from holistic approaches that allow the comprehensive characterization of HM composition, thereby contributing to robust evidence‐based knowledge. As HM composition is known to vary depending on many factors, a comparative analysis including a wide range of compounds across a large sample size would be of great interest. Therefore, the objective of the present study was to perform a comprehensive targeted metabolite analysis of HM using the Biocrates MxP Quant 500 kit. This analysis was designed to provide an overview of the overall composition and quantification of specific HM metabolites and to determine associations with maternal and perinatal characteristics.

The authors declare no conflict of interest.

Supporting Information‐Figures S1_S2 Supporting Information‐Table S1 Supporting Information‐Table S2 Supporting Information‐Table S3