The Mother’s Milk Study is a longitudinal cohort of Latino mother-infant dyads from Southern California, which is investigating the associations between HMOs, the infant gut microbiota, and early life growth and development. Detailed methods have been previously described (40, 41). Briefly, participants were recruited from Los Angeles County maternity clinics affiliated with the University of Southern California. Individuals were eligible to participate if they were 18 years of age or older at time of delivery; had a healthy, singleton birth; were enrolled by approximately one month postpartum; and were able to read at the 5^th grade level in either Spanish or English. Individuals were excluded if they had any medical diagnoses or were taking medications known to affect physical or mental health, nutritional status, or metabolism; were current tobacco users (>1 cigarette in the past week); reported recreational drug use; had pre-term or low birth weight infants; or had infants with any clinically diagnosed fetal abnormalities. The Institutional Review Boards of the University of Southern California, Children’s Hospital of Los Angeles, and the University of Colorado Boulder approved study procedures. Written informed consent was obtained from participants prior to study enrollment.

Participants had visits at 1-, 6-, 12-, 18-, and 24-months postpartum. The Mother’s Milk Study had 219 mother-infant dyads, 209 of whom contributed a human milk sample at the 1-month visit. Of these, funding was available to assess EV-miRNA levels in 110 human milk samples. Among these 110, 11 were removed since they were classified as HMO “non-secretors” (produced under 500 nmol/mL of the HMO 2’FL). We elected to exclude non-secretors because of their low abundance in our sample. An additional participant was excluded as they produced an intermediate amount of 2’FL (559 nmol/mL, considered a “non-secretor” by some definitions), and had HMO diversity 3 standard deviations (SD) from the mean. This resulted in a final sample size of 98. Participants who were excluded from EV-miRNA analysis were similar to those who were included (40).

The clinical measures used for this study were derived from the 1-month postpartum visit. Maternal weight (kg) was measured using an electronic scale and standing height was measured using a stadiometer (m) to calculate maternal body mass index (BMI, kg/m²). Additionally, maternal age at delivery, infant sex, delivery mode, days postpartum, and number of breast feedings per day were collected at the 1-month postpartum visit. Questionnaires were used to determine breast feedings per day, where participants selected 0-1, 1, 2, 3, 4, 5, 6, 7, or ≥ 8 feedings/day and this was treated as a continuous variable. Non-consecutive 24-hour dietary recalls were performed to represent average maternal dietary intake, and from this, the healthy eating index (HEI) was calculated (42). HEI is a composite dietary measure which assesses how well dietary intake aligns with the Dietary Guidelines for Americans (43). Socioeconomic status was measured using a modified version of the Hollingshead Index as previously reported (44).

Mothers fasted for at least 1 hour before human milk samples were collected, and sample collection occurred at least 1.5 hours after the most recent feeding. Participants were instructed to provide a single full expression from the right breast using an electric breast pump, ensuring the collection of fore-, mid-, and hind-milk as previously described (45). Milk was frozen and stored at -80°C until analysis. The Bode Lab at the University of California San Diego conducted HMO analysis. For each participant, one 500μL aliquot of human milk was shipped on dry ice. HMOs are extremely stable, and do not degrade during repeated freeze/thaw cycles or pasteurization (46). The HMO analysis has been previously described in detail (23). Briefly, HPLC after fluorescent derivatization (Vanquish Quaternary HPLC with fluorescent detection, Thermo Fisher Scientific) was used for HMO analysis. Lipids, proteins, salts, etc., were removed using solid phase extraction over C18 and Carbograph. Next, the reducing end of oligosaccharides were labeled with 2-aminobenzamide and removed excess label by solid phase extraction over silica. To account for analyte loss during the extraction procedure, raffinose was added at the beginning of sample processing. Overall, 19 HMOs were identified and quantified: 2′-fucosyllactose (2’FL), 3-fucosyllactose (3FL), 3′-sialyllactose (3’SL), 6′-sialyllactose (6’SL), difucosyllactose (DFLac), difucosyllacto-N-hexaose (DFLNH), difucosyllacto-N-tetrose (DFLNT), disialyllacto-N-hexaose (DSLNH), disialyllacto-N-tetraose (DSLNT), fucodisialyllacto-N-hexaose (FDSLNH), fucosyllacto-N-hexaose (FLNH), lacto-N-fucopentaose (LNFP) I, LNFP II, LNFP III, lacto-N-hexaose (LNH), lacto-N-neotetraose (LNnT), lacto-N-tetrose (LNT), sialyl-lacto-N-tetraose b (LSTb), and sialyl-lacto-N-tetraose c (LSTc). These are 19 of the most abundant HMOs, and they represent more than 95% of the total HMO concentration and all structural features of HMOs, including chain elongation, branching, and all known types of fucosylation and sialylation. HMO concentrations are reported in nmol/mL. HMO diversity was estimated using Simpson’s diversity. The summary measure, sum of all HMOs in a sample, is the sum of all HMOs quantified in each sample. HMO-bound fucose and HMO-bound sialic acid are the sum of all sialic acid and all fucose molecules bound to HMOs in a sample, respectively (e.g., each molecule of 2’FL contains 1 molecule of fucose, and each molecule of DFLNT contains 2 molecules of fucose).

EVs were isolated from stored human milk samples as previously described (47). Milk samples were thawed on ice prior to analyses. Next, samples were centrifuged to remove the lipid layer, then again to remove cellular debris and apoptotic bodies. EVs were extracted using the ExoEasy Maxi KIT (Qiagen, Germantown, MD) and total RNA was isolated with the miRNeasy Serum/Plasma Maxi KIT (Qiagen, Germantown, MD). Samples were cleaned using the RNA Clean & Concentrator-5 Kit (Zymo Research, Irvine, CA) and sample purity and quantity were measured on an Implen NanoPhotometer spectrophotometer (München, Germany). As previously described (40), concentration, sizes, and distribution were assessed in four randomly selected samples using nanoparticle tracking analysis on the ViewSizer 3000 (Horiba Scientific, Piscataway, NJ). The Exo-Check Exosome Antibody Assay (System Biosciences, Palo Alto, CA) was used to confirm the presence and purity of isolated EVs on four sets of three pooled EV samples and three sets of three matched EV-depleted samples. All relevant EV characterization data have been submitted to the EV-TRACK knowledgebase (EV-TRACK ID: EV220416) (48).

Sequencing and library preparation was performed at the University of California San Diego. The NEBNext Small RNA Library Prep Set for Illumina (NEB, Ipswich, MA) was used to construct sequencing libraries with optimization to account for low input and cell-free RNA. Reactions were conducted at one-fifth the recommended volume, adapters were diluted 1:6, and library amplification PCR used 17 cycles. Libraries were cleaned with the DNA Clean and Concentrator Kit (Zymo Research, Irvine, CA) and the concentrations were quantified using the Quant-iT PicoGreen dsDNA Assay (Invitrogen, Waltham, MA). Samples were pooled with equal volumes. The pool’s size distribution was determined with a DNA HS Chip on a BioAnalyzer (Agilent Technologies, Santa Clara, CA) before size selection (115-150 base pairs [bp]) on a Pippin Prep instrument (Sage, Beverly, MA) to remove adapter dimers and fragments larger than the target miRNA population. Libraries were sequenced to ~1 million total reads per pool using a MiSeq instrument with a Nano flow cell (Illumina Inc, San Diego, CA). This sequencing data was used to balance the samples into new pools for deeper sequencing on a HiSeq4000 instrument using single-end 75 bp runs.

Sequencing data were mapped using the ExceRpt small RNA sequencing data analysis pipeline on the Genboree Workbench (49). Samples were mapped using default parameters, except for filtering to a minimum read length of 15 nucleotides with zero mismatches. Quality control was performed according to External RNA Controls Consortium guidelines (49). One sample had <100,000 transcriptome reads and was removed from subsequent analysis. Raw EV-miRNA read counts were normalized using the trimmed mean of M (TMM) method from the EdgeR package (50).

Descriptive statistics were calculated using the mean and standard deviation (SD) for continuous variables and frequencies for categorical variables. We examined the associations between human milk EV-miRNAs and HMOs using two approaches. First, due to the high-dimensional EV-miRNA data, we summarized EV-miRNA levels using principal components (PC) analysis. PCs were then used to assess whether EV-miRNA profiles were associated with HMOs (i.e., diversity, HMO-bound sialic acid, HMO-bound fucose, and individual HMO concentrations). To accomplish this, multivariable linear regression was used to estimate the independent associations between our primary predictors (PC1 and PC2) with HMO measures. Prior to principal component transformation, values were centered and scaled to have unit variance. As a second approach, multivariable linear regression models were used to estimate the associations between individual EV-miRNAs with human milk HMOs. Based on our previous work, all models adjusted for technical covariates (i.e., proportion of rRNA, volume of skim milk) as well as days postpartum, time of day of milk collection, and number of breast feedings per day (40). As a sensitivity analysis, we additionally adjusted for maternal variables age, BMI, and HEI. Unadjusted p-values are reported; however, all analyses were also adjusted for multiple testing using the Benjamini-Hochberg procedure with a threshold of P_BH < 0.10 (51).

Pathway analysis was used to characterize putative mRNA targets of EV-miRNAs that were associated with HMO expression using DIANA MirPath version 4 (52). For pathway analysis, we converted all precursor EV-miRNAs into their mature counterparts, which applied only to hsa-mir-378c. Tarbase v8.0 (53), a catalogue of experimentally validated miRNA-gene interactions, and microT-CDS (54), which predicts in silico miRNA-gene interactions, were used to identify Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with significant enrichment.

Table 1 shows the mean physical and social characteristics of the 98 mothers included in the analysis. The Hollingshead Index was used to measure the participant’s socioeconomic status, where 54% had low or very low SES (Hollingshead score < 26.5). On average, mothers were 27.9 years of age at the time of birth (range: 18-42 years) and overweight at 1-month postpartum (BMI: 29.9 ± 4.8 kg/m²). Additionally, at 1-month postpartum, 16.3% of mothers had a healthy weight, 38.8% had overweight, and 44.9% had obesity. Mothers self-reported that they were breastfeeding an average of 6.3 times per day (range: 0-8) and were in the mature milk stage where, on average, human milk samples were collected 32.5 days after delivery. The top five most abundant EV-miRNAs were miR-148a-3p, miR-146b-5p, miR-200a-3p, let-7g-5p, and let-7b-5p, which is largely consistent with previous studies (55). These EV-miRNAs were detected in every sample and comprised 15.5% of the total reads. As expected, the concentrations of several human milk HMOs were both positively and negatively correlated with one another at 1-month postpartum ( Figure 1 ). Table 2 shows the characteristics of HMO summary measures and individual concentrations at 1-month postpartum.

PC analysis was used to assess whether EV-miRNA profiles were associated with HMO summary measures and concentrations ( Supplementary Figure 1 ). Here, we focused on PC1 and PC2, which explained a total of 34.7% of the variation in EV-miRNA measures ( Supplementary Figure 2 ). We did not include subsequent PCs because they explained little additional variance (e.g., PC3 explained 5.9%). Using multivariable linear models, we constructed models that included both PC1 and PC2 while adjusting for our a priori covariates. PC2 was inversely associated with HMO diversity (P = 0.046), positively associated with the sum of all HMOs (P = 0.0001), and positively associated with HMO-bound fucose (P = 0.02). The EV-miRNAs that were most positively associated with PC2 included miR-183-5p, miR-151a-3p, miR-511-5p, miR-99b-5p, and miR-3615. Conversely, the EV-miRNAs that were most inversely associated with PC2 included miR-502-3p, miR-629-5p, miR-146a-5p, miR-500a-3p, and miR-21-5p ( Supplementary Table 1 ). Additionally, PC2 was positively associated with the concentration of three HMOS, 2’FL, 6’SL, and FLNH ( Table 3 , all P_BH < 0.10). We did not observe any statistically significant associations between PC1 with HMO measures. PC1 was most positively associated with levels of miR-183-5p, miR-151a-3p, miR-511-5p, miR-99b-5p, and miR-3615 and was most negatively associated with levels of miR-21-5p, miR-500a-3p, miR-146a-5p, miR-629-5p, and miR-502-3p.

In sensitivity analyses, we additionally adjusted for maternal BMI, maternal healthy eating index, and maternal age ( Supplementary Table 2 ). Results from these analyses were largely similar to results from the main analyses, where the sum of all HMOs was positively associated with PC2 (P = 0.0004) and positively associated with HMO-bound fucose (P = 0.04). While not statistically significant (P = 0.09), the estimated association between PC2 and HMO diversity was negative. We also found that PC2 was positively associated with 2’FL and 6’SL.

As shown in Table 4 , multivariable linear regression analysis, which adjusted for technical covariates, days postpartum, time of day of human milk collection, and number of breast feedings per day, revealed that 26 EV-miRNAs were associated with HMO diversity (all P < 0.05), 55 with the sum of all HMOs, 9 with HMO-bound sialic acid, and 33 with HMO-bound fucose. Complete results are shown in Supplementary Table 3 . When examining individual HMO concentrations of 2’FL, 3FL, 3’SL, 6’SL, FLNH, LNFP I, and LNH, each were statistically significantly associated with levels of 7, 17, 1, 1, 4, 1, and 5 individual EV-miRNAs (respectively) after correction for multiple testing (P_BH < 0.1).

In sensitivity analyses where we additionally adjusted for maternal BMI, maternal healthy eating index (HEI), and maternal age, findings were largely consistent ( Supplementary Tables 4 , 5 ). Among the 26 EV-miRNAs that were significantly associated with HMO diversity, 19 remained significantly associated after these additional adjustments. Among the 55 EV-miRNAs significantly associated with the sum of all HMOs in the main analysis, 49 remained significant in the sensitivity analysis. Among the 9 EV-miRNAs associated with HMO-bound sialic acid in the main analysis, 4 remained significant in the sensitivity analysis. Among the 33 EV-miRNAs which were associated with HMO-bound fucose in the main analysis, 27 of those remained significantly associated in the sensitivity analysis. Among the 7 EV-miRNAs associated with 2’FL concentration, 6 remained significant in the sensitivity analysis. Among the 17 EV-miRNAs which were associated with 3’FL concentration in the main analysis, 16 remained significant in the sensitivity analysis. Results for the 3’SL, 6’SL, and LNFP I were unchanged in sensitivity analyses. However, among the 4 EV-miRNAs associated with FLNH in the main analysis, only one EV-miRNA remained significant in the sensitivity analysis. Among the 5 EV-miRNAs associated with LNH in the main analysis, 4 remained significantly associated with LNH in the sensitivity analysis.

Finally, we explored the functional pathways of predicted target genes for the EV-miRNAs that were associated with HMO expression. For example, we found that EV-miRNAs that were associated with HMO diversity, 2’FL, FLNH, LNFP I and LNH were also target genes involved in the Hippo signaling pathway ( Table 5 ). EV-miRNAs associated with sum of HMOs, HMO-bound fucose, 3FL, 6’SL, FLNH, and LNH also target genes involved in the cell cycle pathway. Using the micro-T database, similar pathways were identified as being potential targets of EV-miRNAs that were associated with HMO measures ( Supplementary Table 6 ).