Associations between maternal lifetime stressors and negative events in pregnancy and breast milk-derived extracellular vesicle microRNAs in the programming of intergenerational stress mechanisms (PRISM) pregnancy cohort

ABSTRACT

Maternal stress is associated with adverse child health. Breast milk microRNAs encapsulated in extracellular vesicles (EVs) are involved in mother-infant biochemical communication during early-life programming. We leverage the PRogramming of Intergenerational Stress Mechanisms (PRISM) pregnancy cohort to investigate associations between maternal stress and breast milk EV-microRNAs. Lifetime stress and negative life events (NLEs) during pregnancy were assessed using the Life Stressor Checklist-Revised (LSCR) and the Crisis in Family Systems-Revised surveys, respectively. RNA was extracted from breast milk EVs (N = 80; collected 6.1 ± 5.9 weeks postnatally), and microRNAs were profiled using the TaqMan OpenArray Human miRNA panel. Associations between stress scores and detection (yes/no) of 173 microRNAs identified in 20–80% of samples were assessed using logistic regression; associations with expression levels of 205 EV-microRNAs identified in >50% of samples were assessed using linear regression. In adjusted models, detection of 60 and 44 EV-microRNAs was associated with higher LSCR and NLE scores, respectively (p < 0.05). Expression level of 8 and 17 EV-microRNAs was associated with LSCR and NLE scores, respectively, at our a priori criteria of p < 0.05 and |B_regression|>0.2. Enriched KEGG pathways for microRNAs associated with stress scores included fatty acid metabolism and the Hippo signaling pathway. Maternal lifetime stress and NLEs during pregnancy were both associated with detection and expression level of breast milk EV-microRNAs, although associations with microRNA profiles differed between stress measures. Further research is needed to identify biological pathways impacted by associated microRNAs and investigate relationships with child health outcomes.

Abbreviations: EV: extracellular vesicle; PRISM: PRogramming of Intergenerational Stress Mechanisms pregnancy cohort; LSCR: Life Stressor Checklist-Revised survey; NLE: negative life event; CRISYS-R: Crisis in Family Systems-Revised survey; KEGG: Kyoto Encyclopaedia of Genes and Genomes; NYC: New York City; SD: standard deviation; IQR: interquartile range; C_q: relative cycle threshold values; PCA: principal component analysis

KEYWORDS:

• Breast milk
• microRNA
• extracellular vesicles
• negative life events
• stress

Introduction

Maternal psychosocial stressors prior to and during pregnancy may have intergenerational health consequences. Maternal stress has been associated with an array of adverse early life and childhood health outcomes including reduced neurobehavioral and cognitive development [1–3], neurobehavioral disorders[4], lower BMI Z-scores[5], reduced lung function[6], wheeze [7,8], asthma[9], and other atopic diseases[10]. Effects have been observed with stressful and traumatic events that occur both prior to conception [4] and during pregnancy [3–5]. The biological mechanisms through which maternal stress is associated with health outcomes of offspring are not fully understood and are likely multifactorial.

Recent research has focused on the effects of the gestational uterine environment on epigenetic programming. The association between maternal stress and child health may also be mediated postnatally by biochemical communication through breastfeeding [11,12]. There is a large body of research establishing the association between breastfeeding and child health outcomes such as reductions in lower respiratory tract infections, asthma, obesity, and type 1 and 2 diabetes[13]. In addition to nutrients, human breast milk contains a range of regulatory biomolecules including stress hormones, growth factors, cytokines[11], and non-coding RNAs [14] which may be transferred to the infant.

MicroRNAs are small non-coding RNA molecules (~22 nucleotides) that are involved in epigenetic gene regulation[15]. MicroRNAs are transcribed into primary microRNAs and are processed by the ribonuclease Dicer to form mature microRNAs which associate with the RNA-Induced Silencing Complex (RISC). MicroRNAs can pair with mRNA and control post-transcriptional gene repression by inhibiting translation or facilitating degradation of a target mRNA molecule[16]. An individual microRNA may base-pair with the microRNA-recognition elements of hundreds of mRNA molecules[17]. The profile of breast milk microRNAs is different from the profile of microRNAs isolated from maternal plasma, suggesting that they are primarily derived from the mammary epithelium[18]. Breast milk microRNAs can be encapsulated in extracellular vesicles (EVs) derived from breast cells[19]. The term EVs can refer to vesicles on the scale of nanometres of varying size and cellular origin and including exosomes and microvesicles [20].

Breast milk microRNAs encapsulated in EVs are resistant to degradation in the gastrointestinal tract and can be taken up by intestinal cells through endocytotic processes [21,22], and may subsequently affect gene regulation. The bioavailability of milk-derived exosomes and microRNA cargo has been demonstrated in animal models. Ingested milk-derived exosomes accumulated in the liver, spleen, and brain, and labelled microRNAs transfected into milk exosomes accumulated in the intestinal mucosa, spleen, liver, heart, and brain of mice. Biological responses following the ingestion of milk-derived microRNAs have also been identified. Bovine-derived EVs containing RNAs have been associated with changes in the gut microbiome[23], skeletal muscle growth [24], and spatial learning and memory in rodents [25]. It is thought that a large proportion of breast milk-derived microRNAs are related to immune function. Conserved immune-related microRNAs have been identified in milk samples across species including humans, cows, and goats [26]. In mouse models, EVs derived from cow’s milk were found to modify the association between dust exposure and inflammatory responses [27] and mitigate intestinal inflammation due to ulcerative colitis [28]. EVs containing microRNAs may specifically target immune cells through binding of antigen-specific surface antibody light chains [29], facilitating influence on gene regulation with low microRNA concentrations [30].

In human studies, microRNAs related to T- and B-cell maturation and regulation and innate immune response have been found to be highly abundant in breast milk [31–33], including 68% of characterized pre-microRNAs related to immune function [32]. Breast milk-derived microRNAs, including miR-148a and miR-155, may specifically affect the programming of immune response through regulatory T-cell differentiation by modifying the expression of the transcription factor FOXP3 [34,35]. MicroRNAs associated with metabolic pathways have also been identified in human breast milk [33,36].

Previous research characterizing breast milk-derived microRNAs has primarily focused on microRNAs not unique to EVs, and there is limited research addressing factors that influence EV-microRNA composition or investigating associations between EV-microRNA profiles and infant and child health outcomes. The current study sought to contribute to the understanding of breast milk-derived EV-microRNAs as a potential biochemical communication pathway that may mediate associations between maternal stressors and child health. The objectives of this study were to investigate associations between stressful events experienced over a mother’s lifetime or during pregnancy and profiles of EV-microRNAs isolated from breast milk in an ethnically mixed urban pregnancy cohort. Specifically, we assessed associations between maternal lifetime or prenatal stressful events and (1) detection of individual microRNAs (yes/no) and (2) expression levels of individual microRNAs (continuous). We also sought to explore biological pathways regulated by differentially expressed EV-microRNAs by conducting Kyoto Encyclopaedia of Genes and Genomes (KEGG) [37] pathway analyses implemented in DIANA-miRPath v.3.

Methods

Study design and sample selection

The PRogramming of Intergenerational Stress Mechanisms (PRISM) study has been described previously [38]. Briefly, PRISM is a prospective birth cohort of mother-child pairs that was designed to investigate the effects of perinatal stress and other environmental exposures on child development. These analyses include pregnant women recruited from the Icahn School of Medicine at Mount Sinai in New York City (NYC) from November 2012 to August 2014. Inclusion criteria were ≥18 years old, being in the first or second trimester of pregnancy (mean ± standard deviation (SD): 27 ± 8 weeks of gestation), having a singleton birth, and English or Spanish speaking. Women who were HIV+ or those who reported drinking ≥7 alcoholic drinks/week prior to pregnancy recognition or any alcohol during pregnancy were excluded.

Ethics

The protocol was approved by the Institutional Review Board of the Icahn School of Medicine at Mount Sinai. All participants provided written informed consent in their primary language.

Maternal stress measurement

Within two weeks of enrolment, women completed stress questionnaires assessing maternal lifetime exposure to stressful experiences and negative life events (NLEs) during pregnancy. The Life Stressor Checklist-Revised (LSCR) survey ascertains lifetime traumatic and non-traumatic stressful life experiences (e.g., natural disasters, accidents, physical or sexual assault, incarceration, and loss of a child or relative) and includes events specific to women (e.g., abortion) [39]. The questionnaire includes 30 items in a yes/no format, and an open-ended question for any event not listed. The sum of the number of endorsed events was used to calculate the LSCR score with a possible range of 0 to 31. The LSCR has been evaluated for test–retest reliability and validity in diverse populations [39,40].

The Crisis in Family Systems-Revised (CRISYS-R) survey assesses NLEs occurring in the previous six months across 11 domains (financial, legal, career, relationships, safety in the home, safety in the community, medical issues pertaining to self, medical issues pertaining to others, home issues, authority, and prejudice), with multiple items in each domain. The CRYSIS-R has been demonstrated to have good test/retest reliability, and is validated in English and Spanish, and is suitable for lower-income populations [41,42]. Women endorsed items as being positive, negative, or neutral. The number of domains with one or more negative events was summed to create a NLE domain score with a possible range of 0–11 [41].

Sample collection, EV isolation, and RNA extraction

Breast milk was collected from mothers (N = 80) within 6.1 ± 5.9 weeks postnatally. Samples were collected in women’s homes before the first morning feeding. Prior to breast milk collection, the women washed their breast, and 3–10 mL of breast milk was collected using a manual pump. Milk was stored at 4°C and transported to the laboratory the same day. Samples were centrifuged and a layer of fat was removed. Samples were centrifuged a second time and the supernatant was stored at −80°C until assayed when aliquots were thawed on ice prior to EV isolation.

EVs were isolated from breast milk samples collected from 80 women (Supplemental Figure 1) using the exoEasy Maxi Kit (Qiagen, Valencia, CA), which isolates exosomes and other EVs using membrane affinity spin columns[43], according to the manufacturer’s instructions. Briefly, samples were centrifuged and the supernatant was aspirated. The supernatant was mixed with Buffer XBP (1:1 ratio), centrifuged on an exoEasy spin column, and residual liquid was removed. Bound EVs were washed with Buffer XWP, centrifuged, and the flow through was discarded. EVs were eluded using Buffer XE.

Figure 1. A. Volcano plot of associations between EV-microRNA expression and LSCR scores. Estimates from robust linear regression adjusted for infant sex, maternal race, and maternal education. The 8 microRNAs meeting the a priori criteria for significance of p < 0.05 and |B_regression| > 0.2 are plotted in orange. B. Scatter plots of the 8 significant EV-microRNAs and LSCR scores with trend lines from linear models

Total RNA was extracted from EVs using the Qiagen miRNeasy serum/plasma kit. Briefly, QIAzol was used to lyse EVs. miRNeasy Serum/Plasma Spike in Control and chloroform were added and samples were centrifuged to recover the aqueous phase. The aqueous phase was mixed with 100% ethanol and centrifuged on a RNeasy MinElute column. The sample was washed by adding buffer RWT, centrifuging, and discarding the flow through. This process was repeated twice with Buffer RPE. Eighty per cent ethanol was added and following centrifugation the flow through was discarded. RNA was eluted using RNase-free water.

EV-microRNA profiling

RNA was analysed at the Genomics Core Facility at the University of Utah. Prior to analysis, samples were processed using Zymo RNA Clean and Concentrator kits (Zymo Research, Irvine, CA). RNA concentration and quality were measured using the Bioanalyzer (Agilent, Santa Clara, CA). There was sufficient RNA for analysis from 75 women (Supplemental Figure 1). RNA samples were normalized to 34 ng/μl; eleven samples did not have enough RNA for normalization (4.76–33 ng RNA/μl) and were not diluted. RNA was analysed on the TaqMan OpenArray Human miRNA panel (Applied Biosystems, Carlsbad, CA) using 5 μl for each of two primer pools (170 ng). For each of two primer pools, MegaPlex RT reactions and PreAmplification reactions were performed in a 96-well plate on an Applied Biosystems 9700 PCR thermal cycler following the manufacturer’s instructions. The PreAmplification product was diluted to 1:40, mixed with an equal volume of TaqMan OpenArray Real-time PCR Mastermix, and transferred to an OpenArray 384 well plate, with 5 μl of each sample transferred to 8 wells. TaqMan OpenArrays were run on the Applied Biosystems QuantStudio 12 K Flex Real-time PCR System in batches of four plates, with each plate containing three samples. All samples were run on the same day with the exception of nine samples which were rerun due to a defective array.

After removing control probes (ath-miR159a_000338, RNU44_001094, RNU48_001006, and U6 rRNA_001973) and misannotated probes measuring tRNAs (hsa-miR-1274A_002883 and hsa-miR-1274B_002884), the TaqMan OpenArray Human miRNA panel measures expression levels of 752 known human microRNAs (Supplemental Spreadsheet S1).

Covariates

Sociodemographics were collected at the enrolment visit. Self-identified maternal race/ethnicity was classified as Black, Hispanic, or White/other. Maternal education was classified as having a high school diploma or less or completing at least some college. Exposure to tobacco smoke during pregnancy was defined as self-reported active smoking or being exposed to environmental tobacco smoke for ≥ one hour/week during pregnancy. Pre-pregnancy BMI was calculated using self-reported average weight prior to the current pregnancy, and women were classified as not obese (pre-pregnancy BMI < 30 kg/m²) or obese (≥30 kg/m²).

Data analysis

A total of 74 mother-infant pairs were included in analyses due to missing maternal LSCR and NLE variables for one participant (Supplemental Figure 1). Descriptive statistics for study participant characteristics were calculated (median and interquartile range (IQR) for continuous variables and frequency for categorical variables). Differences between mother-infant pairs with breast milk EV-microRNA data and participants in the PRISM cohort overall were assessed using the Mann–Whitney test for continuous variables and Chi-squared test for categorical variables.

Standard quality control procedures were applied to raw qPCR data. Relative cycle threshold values (C_q) associated with amplification scores <1.1 and C_q confidence <0.8 were treated as missing. In addition, Cq values less than the control probe U6 rRNA_001973 (i.e., <10.1) or greater than 35 were treated as missing. A total of 544 EV-microRNAs were detected in ≥ one sample and passed quality control procedures. The number of samples that passed quality control and median (IQR) of the raw C_q values for each EV-microRNA is shown in Supplemental Spreadsheet 1.

To increase statistical power, our primary analyses were restricted to microRNAs detected in >50% of samples (N = 205). Missing C_q values were assumed to be due to non-detectable EV-microRNA levels and replaced with 35 (N = 2,995; 19.4%; the proportion of C_q values imputed in each 205 EV-microRNAs used in analyses is provided in Supplemental Spreadsheet 2). Reference genes were not available and therefore the global mean method [41] was used to normalize data and calculate Delta C_q values (ΔC_q) (i.e., for microRNA_i in sample j with n microRNAs measured, ΔC_{q i,j} = C_{q i,j} – geometric mean of C_{q 1:n,j}). -ΔC_q values were used as the primary unit of analysis, with negative values indicating lower expression relative to the geometric mean and positive values indicating higher expression relative to the geometric mean within a sample. We performed principal component analysis (PCA) to assess the effects of technical variables. The top five PCs were plotted by date run, batch, and low RNA concentration (Supplemental Figures S2-S4). We did not observe notable clustering by batch, date, or low RNA concentration.

Descriptive statistics (median, IQR, mean, and SD) were calculated for each of the microRNAs (Supplemental Spreadsheet 2). Distributions of -ΔC_q values were examined by plotting histograms for each EV-microRNA, and Spearman correlations between EV-microRNAs were calculated. Mann–Whitney and Kruskal–Wallis tests were used to test for associations between characteristics of mother-infant pairs and number of EV-microRNAs detected per sample, and Spearman correlations were used to test for associations between maternal stress variables and number of EV-microRNAs detected. Among microRNAs detected at >20% and <80% of samples, associations between microRNA detection and maternal stress were examined using logistic regression adjusted for infant sex, maternal race, maternal education, and postpartum week of breast milk collection.

Univariate relationships between -ΔC_q values and covariates, i.e., infant sex, maternal age at delivery (<30 years old, ≥30 years old), race/ethnicity (White/other, Black, Hispanic), education (≤ high school degree or GED, > high school degree or GED), pre-pregnancy obesity (BMI ≥30 kg/m²), exposure to smoking during pregnancy (active smoking or exposure to environmental tobacco smoke ≥ one hour/week), and child age at breast milk collection were assessed using robust linear regression. Robust linear regression was also used to assess associations between EV-microRNA expression and LSCR and NLE scores. Models were adjusted for infant sex, maternal race, education, and postpartum week of breast milk collection. Robust linear regression was conducted using the rlm function in R with the MM estimator and bisquare weight function. Due to limited power, to determine the significance of associations with -ΔC_q values we used the a priori criteria incorporating a threshold for effect size [44] of p < 0.05 and |B_regression > 0.2| (i.e., an effect size corresponding to ~1% difference in the range of -ΔC_q values).

We performed enrichment analysis to predict mRNA targets and regulatory roles of differentially expressed EV-microRNAs using KEGG [37] pathway analyses implemented in DIANA-miRPath v.3[45], an online bioinformatics tool. mRNA targets were identified using the Tarbase v.7.0 database, which indexes microRNA–mRNA interactions derived from experimental data[45]. Fisher’s exact test was used to identify enriched pathways at p < 0.05 after false-discovery rate (FDR) correction. Analyses were conducted for microRNAs with expression levels identified to be associated with maternal stress variables at p < 0.05 and |B_regression > 0.2|.

Breast milk EV-microRNA profiles may be variable with postpartum week, and we, therefore, conducted sensitivity analyses excluding samples collected >12 weeks postpartum (N = 5). In addition, due to potential bias introduced by imputed -ΔC_q values, we performed global normalization and subsequent analyses including missing C_q values. All analyses were conducted using R 3.6.1[46].

Results

Characteristics of study participants

Demographic and socioeconomic characteristics of mother-infant pairs (N = 74) are presented in Table 1. The median (IQR) age of mothers at delivery was 27.5 (23.1, 32.1) years, and mothers were predominantly minorities (56.2% Black and 35.6% Hispanic). The majority of mothers (64.4%) had education beyond high school. LSCR and NLE scores ranged from 0 to 19 and from 0 to 9, with a median (IQRs) of 7 (4, 9) and 3 (1, 4), respectively. Measures of maternal stress were moderately correlated (LSCR score and NLE score: r_Spearman = 0.49; p < 0.001). Approximately half of infants were male, 14.9% were preterm, and 12.2% were low birth weight. Breast milk was collected between 0.29 and 28.4 weeks postpartum (median = 4.2 weeks; IQR = 2.6, 7.9 weeks). There were no significant differences between mother-infant pairs included in breast milk EV-microRNA analyses and NYC cohort participants overall in maternal age, race, education, infant sex, gestational age, and birth weight (p > 0.05).

Table 1. Participant characteristics (N = 74)

	N	%
Maternal
Age at delivery (years), median (IQR)	27.5	(23.1, 32.1)
Prepregnancy BMI, median (IQR)^a	25.0	(22.5, 30.9)
Prepregnancy obesity ^a.b	23	32.9%
Education^c
≤ High school diploma or GED	26	35.6%
> High school diploma or GED	47	64.4%
Prenatal tobacco smoke exposure
Environmental tobacco smoke exposure ^d	19	27.5%
Maternal smoking	11	14.9%
Race/ethnicity
White/other	6	8.2%
Black	41	56.2%
Hispanic	26	35.6%
Postpartum week of breast milk collection, median (IQR)	4.2	(2.6, 7.9)
Maternal lifetime stress
LSCR score, median (IQR)^e	7	(4, 9)
Negative life events during pregnancy
NLE score, median (IQR)^f	3	(1, 4)
Infant
Sex
Male (N, %)	41	55.4%
Female (N, %)	33	44.6%
Gestational age (weeks), median (IQR)	39.0	(38.9, 40.0)
Preterm (< 37 weeks gestation)	11	14.9%
Birth weight (g) (mean, SD)	3,050	(2,795, 3,425)
Low birth weight (< 2,500 g)	9	12.2%

a. N = 69. b. BMI ≥ 30. c. N = 73. d. N = 69; being exposed to environmental tobacco smoke for ≥ one hour/week during pregnancy. e. N = 73. f. N = 71.

Characteristics of breast milk EV-microRNAs

A total of 544 EV-microRNAs passed quality control and were detected in at least one sample (Supplemental Spreadsheet S1). The mean (SD) and median (IQR) after imputing missing -ΔC_q values for the 205 EV-microRNAs detected in >50% of samples and used in our primary analyses are presented in Supplemental Spreadsheet 2. Among all microRNA pairs, 32.9% were significantly correlated (Spearman correlation p < 0.05); the majority of significantly correlated microRNAs (62.8%) was positively correlated (Supplemental Figure S5). Between 24 and 308 microRNAs (i.e., 3.2–41.0% of measured microRNAs) were detected per sample (median = 221.5; 29.5%). Among the 205 microRNAs detected in >50% of samples, the number detected per sample ranged from 17 to 205 (median = 183).

Associations between EV-microRNA profiles and participant characteristics

The number of EV-microRNAs detected per sample was not associated with child sex, maternal age at delivery, race/ethnicity, education, obesity, tobacco smoke exposure during pregnancy, postpartum week of breast milk collection (p > 0.05) (Supplemental Table 1). We observed fewer detectable EV-microRNAs among mothers who reported smoking during pregnancy (median among mothers who reported smoking = 188; median among mothers who did not report smoking = 229; Mann–Whitney test p = 0.035).

Univariate associations between EV-microRNA expression and characteristics of mother-infant pairs were assessed using robust linear regression and are reported in Supplemental Spreadsheet S3. We identified 17 microRNAs associated with infant sex, one associated with maternal age ≥30 years, 15 associated with maternal race/ethnicity (10 associated with Black vs. White/other and 7 associated with Hispanic vs. White/other), 19 associated with maternal education, 9 associated with pre-pregnancy maternal obesity, 27 associated with prenatal exposure to environmental tobacco smoke, and 19 associated with prenatal maternal smoking (p < 0.05). Overall, there was minimal overlap in significant EV-microRNAs between variables (Supplemental Figure S6).

We examined associations between microRNA expression and postpartum week of breast milk collection. A total of 59 microRNAs were associated with week of breast milk collection (p < 0.05) (Supplemental Spreadsheet S3). The majority (N = 69; 93%) of breast milk samples were collected within the first 12 weeks postpartum; however, five samples were collected between 22 and 28 weeks postpartum. In sensitivity analyses excluding samples collected >12 weeks postpartum, 28 microRNAs were associated with week of breast milk collection (p < 0.05).

Associations between EV-microRNA profiles and maternal lifetime stress

The number of detected EV-microRNAs per sample was negatively associated with LSCR score (r_Spearman = −0.34, p = 0.003) (Supplemental Table S1). We further evaluated associations between EV-microRNA detection (detected or not detected for each microRNA) and LSCR scores using logistic regression adjusted for infant sex, maternal race, education, and postpartum week of breast milk collection. Among 173 microRNAs detected in >20% of samples and <80% of samples, detection of 60 microRNAs (35%) was associated with higher LSCR scores (p < 0.05 in logistic regression models) (Table 2, Supplemental Figure S7).

Table 2. Summary of associations between EV-microRNAs in breast milk and maternal stress variables

			Associated with LSCR score		Associated with NLE score
Outcome	Number of microRNAs analysed	Criteria for significance	N	(%)	N	(%)
microRNA detection (yes/no)^a	173 detected in > 20% and < 80% of samples	p < 0.05	60	(35%)	44	(25%)
microRNA expression levels^b	205 detected in > 50% of samples	p < 0.05 and |Bregression| > 0.2	8	(4%)	17	(8%)

a. Logistic regression with the outcome of dichotomized detected/not detected for 173 microRNAs detected in > 20% of samples and < 80% of samples and adjusted for infant sex, maternal race, education, and postpartum week of breast milk collection. b. Robust linear regression with the outcome of continuous expression levels (-ΔC_q) and adjusted for infant sex, maternal race, education, and postpartum week of breast milk collection.

In unadjusted analyses of associations between EV-microRNA expression levels and number of endorsed lifetime stressful events, 29 microRNAs were associated with higher LSCR scores (p < 0.05 in robust linear regression models) (Supplemental Spreadsheet S4). After adjusting for infant sex, maternal race/ethnicity, education, and postpartum week of breast milk collection, 21 microRNAs were associated with higher LSCR scores (p < 0.05), of which eight met our a priori criteria for significance of p < 0.05 and |B_regression| > 0.2. There was a positive association between LSCR scores and expression levels of all microRNAs significant at p < 0.05 and |B_regression| > 0.2. Among all microRNAs, we also observed a trend towards a positive association, with 143 (70%) of microRNAs positively associated with LSCR scores.

Associations between EV-microRNA profiles and negative life events during pregnancy

There was a trend towards fewer EV-microRNAs detected with higher NLE scores; although this did not achieve statistical significance (r_Spearman = −0.22, p = 0.062). Detection of 44 (25%) of the 173 EV-microRNAs detected in >20% of samples and <80% of samples was associated with NLE score (p < 0.05 in adjusted logistic regression models) (Table 2, Supplemental Figure S8).

Expression levels of 37 EV-microRNAs were associated with NLE scores in unadjusted models (p < 0.05 in robust linear regression models) (Supplemental Spreadsheet S5). In adjusted models, 23 microRNAs were nominally associated with NLE scores (p < 0.05) Figure 2, of which 17 were significant at p < 0.05 and |B_regression| > 0.2. There was a trend towards a positive association between NLE score and microRNA expression; higher NLE score was associated with increased expression of 14 (82%) of the significant microRNAs and 121 (59%) of the microRNAs overall.

Figure 2. A. Volcano plot of associations between EV-microRNA expression and NLE scores. Estimates from robust linear regression adjusted for infant sex, maternal race, and maternal education. The microRNA meeting the a priori criteria for significance of p < 0.05 and |B_regression| > 0.2 is plotted in orange. B. Scatter plot of the significant EV-microRNA and NLE scores with trend line from linear models

We observed overlap between microRNAs that were significant in both measures of maternal stress as well as consistent directions of association. Six microRNAs were positively associated with both LSCR and NLE scores (hsa-miR-99b-3p: B_LSCR = 0.63, p_LSCR = 0.007; B_NLE = 0.89, p_NLE = 0.035; hsa-miR-96-5p: B_LSCR = 0.25, p_LSCR = 0.039; B_NLE = 0.83, p_NLE < 0.001; hsa-miR-550a-5p: B_LSCR = 0.41, p_LSCR = 0.007; B_NLE = 0.80, p_NLE = 0.003: hsa-miR-616-5p: B_LSCR = 0.26, p_LSCR = 0.025; B_NLE = 0.53, p_NLE = 0.012; hsa-miR-155-5p: B_LSCR = 0.31, p_LSCR = 0.005; B_NLE = 0.49, p_NLE = 0.033; hsa-miR-604: B_LSCR = 0.31, p_LSCR = 0.040; B_NLE = 0.56, p_NLE = 0.042).

Sensitivity analyses

We conducted sensitivity analyses excluding five samples collected >12 weeks postpartum. Overall, analyses of associations between EV-microRNA expression levels and maternal stress were consistent with our primary analyses. All EV-microRNAs associated with maternal stress variables in our primary analyses achieved nominal significance (p < 0.05) with the exception of two microRNAs associated with LSCR score (hsa-miR-99b-3p: B = 0.45, p = 0.06 and hsa-miR-604: B = 0.26; p = 0.10) and one microRNA associated with NLE score (hsa-miR-155-5p: B = 0.44; p = 0.06). Effect estimates in sensitivity analyses are shown in Supplemental Table S2 and S3.

Due to the high proportion of imputed -ΔC_q values among some microRNAs, we also conducted sensitivity analyses treating imputed -ΔC_q values as missing. Data were renormalized prior to analysis and robust linear regression was adjusted for infant sex, maternal race/ethnicity, education, and postpartum week of breast milk collection. Among the eight EV-microRNAs identified to be associated with LSCR scores (p < 0.05 and |B_regression| > 0.2), two achieved nominal significance in sensitivity analyses (hsa-miR-206: B = 0.17, p = 0.014; hsa-miR-99b-3p: B = 0.24; p = 0.026) (Supplemental Table S4). However, the positive direction of association was consistent among all eight EV-microRNAs. Five of the 17 EV-microRNAs associated with NLE score were nominally significant when data were analysed with missing -ΔC_q values (hsa-miR-142-3p: B = −0.36; p = 0.049; hsa-miR-374a-5p: B = −0.20; p = 0.027; hsa-miR-497-5p: B = 0.17; p = 0.022; hsa-miR-99b-3p: B = 0.43; p = 0.01) (Supplemental Table S5).

Pathway analyses

Pathways that are predicted regulatory targets of differentially expressed EV-microRNAs were explored using DIANA-miRPath v.3. In enrichment analysis of the eight microRNAs associated with the LSCR score, 19 KEGG pathways were identified (FDR < 0.05); in analysis of the 17 microRNAs associated with the NLE score, 30 KEGG pathways were identified (FDR < 0.05). Heatmaps of enriched pathways associated with LSCR and NLE scores are provided in Supplemental Figures S9 and S10, respectively. For both LSCR and NLE scores, top enriched pathways included fatty acid biosynthesis (hsa00061) and fatty acid metabolism (hsa01212), which are predicted targets of hsa-miR-604, −16-5p, −182-5p, −374a-5p, and −155-5p; Hippo signalling pathway (hsa04390), a predicted target of hsa-miR-9-3p, −374a-5p, −550a-5p, −497-5p, −16-5p, −616-5p, −99b-3p, and −182-5p; and steroid biosynthesis (hsa00100), a predicted target of hsa-miR-16-5p, −629-5p, −23a-5p, −155-5p, −550a-5p, −9-3p, and −96-5p.

Discussion

To our knowledge, this is the largest study characterizing profiles of breast milk-derived EV-microRNAs and the first study to investigate associations between psychosocial maternal stress and EV-microRNAs. Among 74 mother-infant pairs in the PRISM pregnancy cohort, we detected 544 microRNAs. The number of detected EV-microRNAs per sample was negatively associated with maternal-reported lifetime stress (p < 0.05). Among microRNAs detected in >20% and <80% of samples, detection of 35% was associated with higher lifetime maternal stress scores and detection of 25% was associated with an increased number of negative life events reported during pregnancy.

Among 205 EV-microRNAs detected in at least 50% of samples, expression levels of eight were associated with LSCR score, and 14 were associated with NLE scores (p < 0.05 and |B_regression| > 0.2). Overall, EV-microRNA expression levels were positively associated with measures of maternal stress. Increased expression of six microRNAs (hsa-miR-99b-3p, −96-5p, −550a-5p, −616-5p, −155-5p, and −604) was associated with higher LSCR and NLE scores.

Taken together, these results suggest that maternal stress, both over the life course and during pregnancy, is associated with changes of breast milk-derived EV-microRNA levels. Stressful or traumatic events that occur more distal or proximal to the time of birth may have varying consequences on the microRNA profile: detection of a greater number of microRNAs was associated with LSCR scores whereas expression levels of a greater number of microRNAs were associated with NLE scores.

Although the biological mechanism through which maternal stress may affect breast milk-derived EV-microRNA profiles is speculative, recent research has focused on the association between maternal stress and dysregulation of the hypothalamic–pituitary–adrenal (HPA)-axis due to changes in the circulating levels of glucocorticoids, including cortisol[47]. Cortisol may directly affect child health by being transferred from mothers to infants through breast milk[47]. In longitudinal studies, cortisol levels in breast milk have been negatively associated with body mass index percentile in the first two years of life[48], and higher maternal plasma cortisol levels, a surrogate for breast milk cortisol, have been associated with fearful temperament among breastfed infants, but not among formula-fed infants [49]. MicroRNA levels have been related with stress among adults [47,50], suggesting that maternal cortisol levels may have epigenetic consequences, including dysregulation of microRNAs.

We explored potential regulatory targets of differentially expressed EV-microRNAs using KEGG enrichment analysis (Supplemental Figures S7 and S8). The Hippo signalling pathway, implicated as a pathway regulated by hsa-miR-99b-3p, is involved in cell proliferation and organ development[51]. Dysregulation of this pathway has been associated with cancer development and may be involved in a broad range of metabolic diseases[52]. Activation of glucocorticoid receptors may result in a molecular cascade to affect levels of yes-associated protein (YAP), an endpoint of the Hippo signalling pathway[53]. Additional pathways identified included fatty acid biosynthesis and fatty acid metabolism, which may be influenced by glucocorticoids. Glucocorticoids activate biological responses to stress; however, chronic elevation of glucocorticoid levels may result in adverse health consequences. Although the effects of cortisol and other glucocorticoids on fatty acid metabolism are not well characterized, in vitro and in vivo studies suggest that acute and chronic glucocorticoid exposure can influence multiple components of fatty acid metabolism in humans, including lipolysis, lipogenesis, and fatty acid oxidation[54]. Pathways related to immune response, including FOXO signalling[55], PI3K-AKT signalling[56], bacterial invasion of epithelial cells were also identified, suggesting a role of breast milk-derived EV-microRNAs in mediating associations between maternal stress and child respiratory outcomes [6–9].

Among the microRNAs associated with maternal stress, hsa-miR-96-5p may be related to stress, postnatal development, and cognitive function. mRNA targets of hsa-miR-96 include mammalian target of rapamycin (mTOR) [57] and AKT1 Substrate 1 (AKT1S) [58]. mTOR, which is regulated by AKT1S, is a kinase that forms protein complexes including mTORC1, which is involved in protein, lipid, and nucleotide synthesis, autophagy[59], and regulating metabolic activity, including insulin secretion of pancreatic β-cells[60]. Amino acids, glutamine, palmitic acid, and exosomal microRNAs found in breast milk also regulate mTORC1 and AKT1S levels in neonates [60–63]. In addition, mTOR is involved in neurological processes and is suppressed by chronic unpredictable stress in mouse models[64]. There is also evidence that hsa-miR-155-5p is related to response to psychological stress. Exposure to stress increases the production of pro-inflammatory cytokines in humans[65], and miR-155 inhibition has been shown to decrease cytokine production and oxidative stress and improve neurological function in rats following transient global ischaemia. Furthermore, miR-155 targets atopy-related genes and is involved in the pathway for regulatory T-cell differentiation[35]. Taken together, this evidence suggests that hsa-miR-96-5p and hsa-miR-155-5p may be involved in biological pathways linking maternal stress and postnatal health outcomes, including immune and metabolic functions, through breast milk-derived microRNAs[66].

We observed overlap between EV-microRNAs detected in our study and those identified in previous studies characterizing EV-microRNAs in human breast milk. van Herwijnen et al. profiled breast milk-derived EV-microRNAs from four mothers collected 3–9 months after delivery[67]. Among the 10 most abundant microRNAs, we detected seven in ≥50% of our samples (hsa-miR-148a-3p, −30a-5p, −30d-5p, −21-5p, −200a-3p, −200 c-3p, −200b-3p). Zhou at al. measured EV-microRNAs from four mothers collected two months postpartum; among the top 10 abundant microRNAs, seven were detected in at least 50% of our samples (hsa-miR-629-5p, −30b-5p, −146b-5p, −29a-3p, −141-3p, −182-5p, and −200a-3p)[68]. Simpson et al. collected breast milk 3 months postpartum from 54 mothers, and seven of the top 10 abundant EV-microRNAs were detected in at least half of our samples (hsa-miR-148a-3p, −22-3p, −30d-5p, −200a-3p, −146b-5p, −24-3p, −21-5p)[69]. Liu et al. measured EV-microRNAs in breast milk collected from 12 mothers between 1.5 and 8 months postpartum, and among the 10 abundant EV-microRNAs, nine were detected in our samples (hsa-miR-22-3p, −148a-3, −30d-5p, −181a-5p, −141-3p, −26a-5p, −30b-5p, −92a-3p, −182-5p)[70]. Differences in microRNAs identified may be due in part to participant characteristics, time of breast milk collection, and EV isolation techniques.

This study was strengthened by measuring microRNAs using the TaqMan OpenArray Human miRNA targeting 752 microRNAs, which allowed us to efficiently measure known human microRNAs in a larger number of samples than previous studies. This technology is based on real-time PCR – the gold standard for transcript measurement – and has been shown to be superior to RNA sequencing for detecting and quantifying miRNAs[71]. In addition, we measured microRNAs isolated from EVs, which are more likely to have downstream biological effects in infants. Encapsulation in EVs protects microRNAs from the environment of the gastrointestinal tract and within circulation, increasing the potential for uptake by receptor cells[72].

Some limitations of our study are also worth noting. Sample processing included removing the fat layer after centrifugation of breast milk. MicroRNAs have been identified in the lipid fraction of human breast milk[36]. Although this method may have removed a source of microRNAs, it allowed us to minimize variation in microRNA profiles due to differences in the fat:skim milk ratio. In addition, EVs were isolated using the exoEasy Maxi Kit, which isolates EVs and exosomes using membrane affinity spin columns[43]. This method does not distinguish by size or cellular origin; however, we have previously demonstrated the presence of common surface markers (CD63 and CD9) and consistent size distribution (mean diameter ~100 nm) among isolated EVs[14], and the use of the exoEasy Maxi Kit is advantageous over affinity-based methods in that it does not depend on a specific surface antigen[73]. Although we had a larger sample size than previous studies of human EV-microRNAs in breast milk, this study was limited by a small sample size. To increase power, we chose to use the a priori criteria for significance incorporating effect size [44] of p < 0.05 and |B_regression| > 0.2 for determining statistically significant associations with EV-microRNA expression levels. Although this approach incorporates effect size in determining significance, it may increase the rate of type I errors. Future studies with larger sample sizes should include false discovery rate (FDR) or Bonferroni corrections to adjust for multiple testing. In addition, breast milk was collected between 0.29 and 28.4 weeks postpartum. We observed a high rate of associations between postpartum week of breast milk collection and EV-microRNA expression levels (29% of microRNAs, p < 0.05), and overlap between microRNAs associated week of breast milk collection and maternal stress (hsa-miR-99b-3p associated with LSCR and NLE scores; hsa-miR-9-3p associated with LSCR scores; and hsa-miR-142-3p, hsa-miR-515-3p, and hsa-miR-16-5p associated with NLE scores; p < 0.05 and |B_regression| > 0.2), suggesting that week of breast milk collection could act as an effect modifier in the association between maternal stress and microRNA expression. However, sensitivity analyses excluding five samples collected >12 weeks postpartum were robust overall. The direction of association was consistent across all EV-microRNAs significant in our primary analyses. Two microRNAs associated with LSCR scores and one microRNA associated with NLE scores failed to achieve significance, possibly due to a reduced sample size.

We also conducted sensitivity analyses treating imputed -ΔC_q values as missing. It should be noted that we observed limited consistency with our primary analyses: two of the eight EV-microRNAs identified as associated with LSCR scores were nominally significant in our sensitivity analyses and five of the 17 EV-microRNAs associated with NLE scores were nominally significant (p < 0.05). This limited overlap, however, may be due to small sample sizes in individual models after excluding imputed values. Visual inspection of scatter plots of -ΔC_q values and stress scores for significant EV-microRNAs Figure 1–3 suggests that the expression of these microRNAs is linearly associated with stress scores.

Conclusions

Maternal microRNAs encapsulated in EVs may reach the infant through breast milk and are a novel biochemical communication pathway for early-life programming. However, few studies have investigated EV-microRNA profiles isolated from breast milk or factors that may influence breast milk EV-microRNAs. In this study, we find that lifetime maternal stress and stressful events experienced during pregnancy were associated with differential detection and expression of breast milk-derived EV-microRNAs. Biological pathways enriched among differentially expressed microRNAs, including fatty acid metabolism and the Hippo signalling pathway, may be affected by maternal cortisol levels. Further research is necessary to understand biological mechanisms linking maternal stress to EV-microRNA profiles and to investigate potential impacts of breast milk-derived EV-microRNAs on infant and child health.

Acknowledgments

The PRogramming of Intergenerational Stress Mechanisms (PRISM) cohort has been supported under US National Institute of Health (NIH) grants R01 HL095606, R01 HL114396, and R01 ES030302 (RJW, Principal Investigator). During the preparation of this manuscript, AGL was supported by US NIH grants R01 MD013310 and K23 HL135349. Analysis of EV-microRNAs was supported in a pilot grant under the US NIH grant P30 ES023515 (AGL, RJW).

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

The PRogramming of Intergenerational Stress Mechanisms (PRISM) cohort has been supported under US National Institute of Health (NIH) grants R01 HL095606, R01 HL114396, and R01 ES030302 (RJW, Principal Investigator). During the preparation of this manuscript, AGL was supported by US NIH grants R01 MD013310 and K23 HL135349. Analysis of EV-microRNAs was supported in a pilot grant under the US NIH grant P30 ES023515 (AGL, RJW);National Heart, Lung, and Blood Institute [R01 HL095606];National Heart, Lung, and Blood Institute [K23 HL135349];National Heart, Lung, and Blood Institute [R01 HL114396];National Institute of Environmental Health Sciences [P30 ES023515];National Institute of Environmental Health Sciences [R01 ES030302];National Institute on Minority Health and Health Disparities [R01 MD013310];