The present study was performed in a subgroup of 75 healthy mother–infant pairs within the MAMI birth cohort, which is a clinical prospective mother–infant birth cohort from the Mediterranean area (48) (Clinical Trial Registry NCT03552939). The subgroup of the MAMI cohort was chosen based on breast milk sample availability at both time points (days 7 and 15) within the transition stage. Women with accessible milk samples who adhered to exclusive breastfeeding practices at the specified time points (7 and 15 days) were considered for inclusion. Our study participants were in good health, with no reported chronic health conditions. Exclusion criteria encompassed non-compliance with any of the inclusion criteria, the use of medications or drugs, the presence of chorioamnionitis, mastitis, or health complications during the gestational and breastfeeding periods, and the presence of any chronic diseases or medication usage for chronic conditions. All mothers filled out a clinical questionnaire about them and their children throughout the study (49). Maternal and neonatal clinical records, including mode of birth, neonate gender, maternal BMI, antibiotic exposure, and others, were recorded. The weight and height of the mothers and the infants were measured in medical consultation at different time points: during pregnancy, at birth, at 7 and 15 days, and at 1, 6, 12, 18, and 24 months postpartum.

Furthermore, the mothers also filled out a 140-item Food Frequency Questionnaire (FFQ) about their regular diet during pregnancy (50). The FFQ information was analyzed to obtain data on the intakes of macronutrients and micronutrients per day. For that, FFQ was previously validated in this same cohort by a 3-day questionnaire (51). The analysis of this information permitted the identification of two different dietary patterns in the MAMI cohort: Diet I and Diet II groups, by clustering using the Jensen–Shannon distance and partitioning around medoid clustering methods (51). Diet I was characterized by a higher intake of total dietary fiber, vegetable protein, and polyunsaturated fatty acids, while Diet II was based on a higher intake of animal protein and saturated fatty acids.

Moreover, 72 of the 75 mothers had their secretor status analyzed. For this, total DNA was obtained as described previously (51) and used to test the secretor status by FUT2 genotyping with the polymerase chain reaction (PCR)-random fragment length polymorphisms (RFLPs) as detailed elsewhere (52). Secretor mothers were those with the positive detection of the FUT2 gene, which encodes for the galactoside 2-alpha-L-fucosyltransferase 2.

Human milk samples for this study were obtained between 7 and 15 days after birth. All the mothers were asked to collect the samples in a specific way: in the morning before lactation, cleaning the breast skin with a solution of 0.5% chlorhexidine, and using a sterile pumper. In addition, the first 2–3 drops were discarded to normalize the collection, and approximately 10 ml were collected. Milk samples were sent to the biobank and kept at −80°C until the lactic serum was obtained by centrifugation at 2,000 × g for 15 min at 4°C. Finally, lactic serum was frozen again at −80°C in different aliquots to carry out all the determinations in one freeze/thaw cycle.

The quantification of Igs (IgE, IgG1, IgG2, IgG3, IgG4), CKs (GM-CSF, IFN-γ, IL-1β, IL-2, IL-4, IL-5, IL-6, IL-9, IL-10, IL-12p70, IL-13, IL-17A, IL-18, IL-21, IL-22, IL-23, IL-27, TNF-α), and adiponectin was performed by ProcartaPlex™ Multiplex immunoassay (Thermo Fisher Scientific, Vienna, Austria) using an Antibody Isotyping 7-Plex Human ProcartaPlex™ panel, a Th1/Th2/Th9/Th17 Cytokine 18-Plex Human ProcartaPlex™ panel, and an Adiponectin Human ProcartaPlex™ Simplex Kit, respectively, following the manufacturer’s instructions as in previous studies (49, 53–58). The plates were run on a Luminex instrument and analyzed in ProcartaPlex Analyst Software (MAGPIX® analyzer, Luminex Corporation) at the Flow Cytometry Unit of the Scientific and Technological Centers of the University of Barcelona (CCiT-UB). Assay sensitivity was as follows: 2.11 ng/ml for IgG1; 16.07 ng/ml for IgG2; 0.08 ng/ml for IgG3; 0.56 ng/ml for IgG4; 0.003 ng/ml for IgE 1.2 pg/ml for GM-CSF; 0.2 pg/ml for IFN-γ; 0.2 pg/ml for IL-1β; 0.8 pg/ml for IL-2; 1.5 pg/ml for IL-4; 0.3 pg/ml for IL-5; 0.4 pg/ml for IL-6; 0.5 pg/ml for IL-9; 0.1 pg/ml for IL-10; 0.04 pg/ml for IL-12p70; 0.1 pg/ml for IL-13; 0.1 pg/ml for IL-17A; 0.4 pg/ml for IL-18; 0.6 pg/ml for IL-21; 8.2 pg/ml for IL-22; 0.9 pg/ml for IL-23; 5.1 pg/ml for IL-27; 0.4 pg/ml for TNF-α; 4.6 pg/ml for adiponectin.

The quantification of leptin was performed by a Quantikine® Colorimetric Sandwich ELISA Kit (R&D Systems, Minneapolis, MN, USA) following the manufacturer’s instructions as in previous studies (41). Data were analyzed by Multiskan Ascent v2.6 software (Thermo Fisher Scientific, Vienna, Austria). Assay sensitivity was 7.8 pg/ml.

The quantification of IgA and IgM was performed by an ELISA kit from Bethyl Laboratories (Montgomery, TX, USA) and an ELISA kit from Cloud-Clone Corp. (Houston, TX, USA), respectively, following the manufacturer’s instructions as in previous studies (31, 59). These kits included the quantification of the secretory form of both types of immunoglobulins. Assay sensitivity was 14.1 ng/ml for IgM, and 1.03–750 ng/ml was the assay range for IgA. Data were analyzed by Multiskan Ascent v2.6 software (Thermo Fisher Scientific, Vienna, Austria).

The results are expressed as mean ± SEM (standard error of the mean) unless otherwise specified. Box plots showed the interquartile range, with the error bars representing the lowest and highest values. Shapiro–Wilk and Levene’s tests were used to determine the normality and homogeneity of the data variance, respectively. When variables were not normally distributed, a logarithmic normalization of the data was made or non-parametric tests were used. Multivariate general linear model (GLM) and repeated measures GLM analyses were conducted to assess the variance in multiple dependent variables using two-factor variables across measurements taken at different time points. For example, this approach was applied to analyze the tracking of infant weight and height up to the first month of life in the two study groups (BM immunotypes). Spearman’s correlation coefficient was used to search for correlation between variables that were not normalized. Student’s t-test and Mann–Whitney U-test were used to assess significant differences between groups, while the chi-square test (X²) compared frequencies, such as the detectability of CKs. For cytokines, a left-censored Tobit model with a threshold of 0 was fitted to test differences between times using the R software package AER (60). The p-values were adjusted for multiple comparisons using the false discovery rate (FDR) correction. A p-value of <0.05 was considered significant in all the tests.

The MFA was performed to integrate the two measurements on days 7 and 15 with the aim of providing a comprehensive understanding of the contribution of the different variables and times to the variability of the data. The MFA models were graphically displayed by plotting the projections of the samples onto the bidimensional space defined by the two first dimensions. All the data were normalized by means of log transformation prior to carrying out these statistical analyses using the R software packages FactoMineR and factoextra (61).

To discover the relationships between milk composition variables and immunotype groups, a mixed graphical model (MGM) was used. MGMs are undirected probabilistic graphical models, where each node corresponds to one variable, and the edges between two nodes represent a conditional dependency between them given all other variables in the graphical model (62). We have used ‘mgm’ R-package to estimate the network of dependencies (63). Specifications were set to allow the maximum number of interactions in the network.

All statistical analyses were performed using IBM SPSS Statistics 22 (IBM, USA) and R version 4.1.2 (R Foundation, Austria) (64).

The mean maternal weight gain during pregnancy was in the range recommended by the Institute of Medicine (65). Overall, approximately 50% of them used antibiotics during pregnancy and/or on delivery day. More than half of the participants’ deliveries were vaginal, and approximately 50% of the total participants were first-time mothers (Table 1). The table includes the breast milk (BM) immunotype classification, depending on the later results.

With regard to the infant characteristics, slightly more than half of the infants were females, 23% took antibiotics from birth to 2 years, and 11 of 32 infants developed atopy. As expected, the intake of antibiotics and the number of infections increased over time (Table 2).

The description of some immune components, particularly the characterization of the IgG isotype profile and IgE levels found in transition milk, is described for the first time in this study. The immune compounds present in milk on days 7 and 15 of lactation were analyzed and compared.

Although there were no differences in adipokine levels (Supplementary Figure 1), the immunoglobulin and cytokine concentrations showed some differences between days 7 and 15 (Figure 1; Table 3, respectively). IgA, IgM, and IgE showed their highest values on day 7. However, IgG levels did not change between milk from day 7 to day 15, and with respect to the IgG isotypes associated with the Th1 or Th2 response, only IgG2 displayed a significant decrease with time (p < 0.05) even though it did not affect the Th1/Th2 ratio. With regard to the cytokine levels, day 7 presented a higher concentration of IL-1β, IL-6, IL-9, IL-10, and IL-17 (p = 0.021, p = 0.020, p = 0.0021, p = 0.020, and p = 0.048, respectively) than day 15. On the contrary, IFN-γ and IL-18 levels were higher on day 15 than on day 7 (p = 0.021 and p = 0.028, respectively). Overall, the cytokine concentration in the samples analyzed was variable, and some cytokines were not detected in some subjects, similar to what has been previously described (66) (Table 3). The components were grouped depending on their functions to study the immune response (pro-inflammatory, anti-inflammatory, Th1, Th2, Th9, Th17, hematopoietic factors, innate immunity, and acquired immunity) (27, 58, 67–75) (Supplementary Figure 2). Some of the calculated variables were also significantly higher on day 7 than on day 15, such as anti-inflammatory response, Th2, Th17, innate immunity, and acquired immunity (Supplementary Figure 3).

A multivariate factor analysis (MFA) was performed to observe different profiles regarding the BM immune composition. The distribution of the mothers taking into consideration each day separately, on days 7 and 15 (Figures 2A,B, respectively), did not show any particular aggrupation. However, when the MFA was generated considering the two types of milk from the same mother, two clusters were identified with strong similarities among individuals (Figure 2C). It was also observed that the mothers from the same cluster had similar dynamic behavior regarding their BM immune components from d7 to d15. The clusters were based on the BM immune components; therefore, they were named BM immunotype I (BM-I) and BM immunotype II (BM-II). Furthermore, as represented in Supplementary Figure 4, considering each day separately and the existence of the two BM immunotypes, the two clusters could already be observed, but without a good separation, as noticed in Figure 3C.

The next step was to characterize these clusters, or BM immunotypes, to define the differences in composition between them and to study the association with maternal and infant characteristics to find potential sources of variability in the composition of milk immune molecules.

The variables that influenced the most separation among the participants were the Ig levels, mainly IgG subtypes, and total IgG. Observing the importance of the variables in each dimension (Dim) defining the MFA, it can be established that while Dim-1 was described mainly by IgG-15d, IgG1-15d, and IgG3-15d, Dim-2 was linked to levels of IgG-7d, IgG1-7d, and IgG2-7d (Figures 2D,E, respectively). In line with these results, these components displayed differing temporal evolutions when compared between clusters, as was observed in the interaction of variables (group and time) by conducting a multivariate general linear model (GLM; Supplementary Table 1). Additionally, differences between groups were also observed when fixing the sampling time for those components that exhibit distinct temporal evolutions. For instance, the BM immunotype II showed higher levels of IgG, IgG1, and IgG2 on day 15 with respect to the BM immunotype I, even though on day 7 they were higher in the BM immunotype I (Supplementary Table 1). The influence of these components in clustering, and of others such as IgE and IL-22 was confirmed when strong relations were directly found in the mixed graphical model (MGM) with the BM immunotype group (Supplementary Figure 5a). In addition, strong correlations were also observed among the predominant variables that described Dim-2 (IgG1-7d, IgG2-7d, IgG3-7d, IgG4-7d, IgG-7d, and IgE-7d; Supplementary Figure 5b). In addition, IgG1-15d and IgG-15d, which are important in Dim-1 were also highly correlated (Supplementary Figure 5b).

Moreover, to investigate the relations between all the components, Spearman’s correlations were made with the non-transformed data. The study of the impact of adipokine levels on the other BM immune components can be observed in Figure 3. Leptin and adiponectin showed significant positive correlations with several BM CKs and Igs (Figure 3).

To evaluate the dynamic behavior of the BM immune components within the BM immunotype group, we calculated the increase of each immune variable with the transformed data from d7 to d15 (Figure 4). Each BM immunotype group had a different dynamic pattern. Importantly, the most prominent components were the strongest in the dimensions of the MFA. As all the IgG subtype levels increased from d7 to d15 in BM-II, the total IgG concentration was also boosted in the BM-II group (Figure 4A). IgE and adiponectin concentrations also showed the same pattern (Figures 4C,D, respectively). The opposite behavior was observed in the BM-I cluster for the same components. On the contrary, while leptin levels increased after 1 week in the BM-I, they decreased in the BM-II (Figure 4D). Importantly, by assessing the IgA and IgM changes, we found that they behaved the same in both groups. Thus, these Igs decreased from d7 to d15 in all mothers (Figure 4B). In terms of cytokine concentration, IL-18, IL-22, IL-23, IFN-γ, and TNF-α also showed significantly different dynamic behaviors when comparing the two groups. For instance, IL-5, IL-18, IL-22, and IFN-γ increased over time in the BM-I group, while they decreased over time in the BM-II group.

To further investigate the potential source of variability between clusters, the maternal characteristics within each group were compared in relation to the BM immunotype. The maternal variables (e.g., pre-gestational BMI, weight gain, the secretory classification, and the presence of animals at home) and the gestational aspects (i.e., gestational age, health complications during pregnancy, and antibiotic use) were similar between the two groups, allowing us to discard these factors as contributors to the clustering of the BM immunotypes. However, mothers belonging to BM-I displayed a lower frequency of being first-time mothers (Table 1) and a tendency to increase the number of mothers practicing exclusive breastfeeding (p = 0.069; Table 1) with respect to BM-II mothers. In line with these results, more mothers started a mixed feeding before day 15 in the BM-II group (p = 0.015) than in the BM-I group. To investigate the impact of the mother’s genetic background with respect to HMO production in the BM immunotype, it was also possible to classify the mothers according to their secretor gene (secretor mother and non-secretor mother). The results showed that the genetic background did not influence the BM immune composition of the studied mothers (Supplementary Figure 6A).

Mothers were also classified into two dietary groups (Supplementary Figure 6B) according to their dietary records, as done in previous studies with the same cohort and as described in the MM section (51, 78). While the Diet I group ate foods richer in fiber and vegetable protein, the Diet II group ate foods richer in animal protein and saturated fatty acids. Although a clear association of the dietary group with the BM immunotype was not found (Supplementary Figure 6B), plant-derived metabolites, such as total polyphenols, tended to increase in the mothers belonging to the BM-II group (p = 0.066; Supplementary Figure 6C); however, vitamin D levels were lower in BM-II with respect to the BM-I group (Supplementary Figure 6D).

To investigate the impact of having different immune dynamic profiles in the transition stage, the infant characteristics between the two groups (i.e., gender, weight, height, BMI z-scores, atopy, antibiotic use, and number of infections) were also analyzed (Table 2). Although children of mothers belonging to the BM-I group seemed to exhibit higher weight and height during the first month of life, the weight trajectory over time did not differ between clusters (Table 2), and the increase in weight and height from birth to the first month of life is similar when comparing the two groups (data not shown).

To assess the influence of each BM immune component on infant growth, the Spearman correlations between these parameters were also studied (Figure 5). Surprisingly, different correlations were found in the two BM sampling days (Figures 5A,B). While negative significant correlations between leptin at 7 days and the weights and heights from birth to 1 month were found (Figure 5A), the leptin levels at 15 days only showed significant negative correlations with the weight and the height at 1 month of age (Figure 5B). Additionally, more correlations between these infant parameters and CKs were found in the BM at 7 days than in the BM at 15 days. Contrary to what happened at 7 days, at 15 days, many important negative correlations were observed between Ig types and infant growth from birth to 1 month (Figure 5B). In addition, it is important to note that most of the correlations with infant parameters were concentrated in the first month of life.

Continuing with the infant characteristics, the offspring of the BM-I group required more antibiotics from birth to 24 months (Table 1), although it is worth noting that the data decreased throughout the follow-up time. In line with these results, it was also observed that the three children who had more infections from birth to 2 years belonged to the BM-I group (Figure 6D). However, there was no association between the presence of infections and the clusters of BM immunotypes (Figures 6A–C).