Metabolism-related microRNAs in maternal breast milk are influenced by premature delivery.

BackgroundMaternal breast milk (MBM) is enriched in microRNAs, factors that regulate protein translation throughout the human body. MBM from mothers of term and preterm infants differs in nutrient, hormone, and bioactive-factor composition, but the microRNA differences between these groups have not been compared. We hypothesized that gestational age at delivery influences microRNA in MBM, particularly microRNAs involved in immunologic and metabolic regulation.MethodsMBM from mothers of premature infants (pMBM) obtained 3-4 weeks post delivery was compared with MBM from mothers of term infants obtained at birth (tColostrum) and 3-4 weeks post delivery (tMBM). The microRNA profile in lipid and skim fractions of each sample was evaluated with high-throughput sequencing.ResultsThe expression profiles of nine microRNAs in lipid and skim pMBM differed from those in tMBM. Gene targets of these microRNAs were functionally related to elemental metabolism and lipid biosynthesis. The microRNA profile of tColostrum was also distinct from that of pMBM, but it clustered closely with tMBM. Twenty-one microRNAs correlated with gestational age demonstrated limited relationships with method of delivery, but not other maternal-infant factors.ConclusionPremature delivery results in a unique MBM microRNA profile with metabolic targets. This suggests that preterm milk may have adaptive functions for growth in premature infants.

BACKGROUND

Maternal breast milk (MBM) provides optimal nutrition for growth and development of infants. However, premature infants face a variety of physiological challenges that differ from infants born at term and therefore have different nutritional needs from maternal milk. When compared to the milk of mothers with term infants (tMBM), the milk of mothers with preterm infants (pMBM) has higher concentrations of protein, nitrogen, fat, phospholipids and elements such as sodium, chloride, and iron (1–2). Conversely, there are lower levels of lactose and zinc in pMBM compared with tMBM (1–3). Macronutrient profiles in MBM also vary with gestation of delivery. An example of this includes lipids: medium chain fatty acids are higher in pMBM, but levels of oleic acid are lower (1). Analysis of the nutritional sufficiency of pMBM shows that nutrients received per volume of MBM exceed the intrauterine requirements for a given gestational age. In addition, preterm infants fed pMBM receive more nutrients than if they were fed tMBM (1–2). The importance of these differences in pMBM and tMBM is highlighted by findings that preterm infants fed pMBM have improved health outcomes compared to preterm infants fed tMBM (4).

Micro ribonucleic acids (miRNAs) are non-coding nucleic acid sequences derived from inter- and intragenic DNA regions. They reduce the efficiency of gene expression, inhibiting protein translation through targeted degradation of mRNA transcripts. A single miRNA can bind to multiple mRNAs exponentially increasing its regulatory activity (5). Although miRNAs regulate numerous processes and are present in biofluids throughout the body, they are found in the greatest concentrations in MBM (6). The miRNAs found in MBM originate primarily from mammary gland epithelium, where they play a role in mammary epithelial cell function (7). A portion may also be derived from extra mammary tissue via the bloodstream (5,8–9). In MBM, miRNAs are packaged into vesicles such as exosomes (6). This accounts for the higher proportion of miRNAs in the lipid fraction of MBM compared to the skim fraction (5,9). Recent studies have identified 281–602 unique miRNAs in human MBM, most of which are found in lipids(5,8,10). Over 65% of MBM miRNAs are related to immune (10) or metabolic function (9,11–13). For example, miR-181a and miR-155 are abundant in breast milk and work together to increase B- and T- cell proliferation, resulting in enhanced adaptive immunity (5,9–10). In addition, miR-33 and miR-122 (which are expressed in MBM) regulate macromolecule synthesis via cholesterol and lipid homeostasis (9,12). Therefore, investigators have suggested that MBM miRNAs affect infant health and development (5,9–10). This idea is supported by a study in which arthritic mice fed milk rich in immune related miRNAs (miR-30a, miR-223, miR-92a) exhibited a reduction in cartilage depletion and inflammation (14).

Although many components of pMBM have been described, no study has examined how pMBM differs from tMBM in its composition of miRNA. When considering nutrient composition differences of pMBM and the proposed importance of miRNAs on infant health, a functional difference in miRNA content seems likely. This study tests the hypothesis that miRNA profiles differ between pMBM and tMBM samples collected one month post-delivery and that miRNAs “altered” in pMBM will have functional relevance to infant metabolism and immune regulation.

METHODS

This study was approved by the Institutional Review Board at the Penn State College of Medicine. All participants provided informed consent before participation.

Term Milk Recruitment and Collection

tMBM was collected as part of a follow up study to the Intervention Nurses Start Infants Growing on Healthy Trajectories (INSIGHT) investigation (15), which examined the influence of genetics and environment on the growth and development of second born siblings. This cohort included multiparous mothers ≥20 years of age that had given birth to singleton, full term infants (≥38 weeks gestation). Exclusion criteria for this group included medical complications with delivery, prenatal ultrasound suggesting intrauterine growth retardation, or newborn birth weight <2,500g. Basic health and demographic information was obtained for each mother infant dyad, including maternal age (years), pre pregnancy weight (kg), pregnancy weight gain (kg), basic metabolic index (kg/m2) at delivery, maternal race/ethnicity, presence or absence of maternal hypertension, delivery method (cesarean section or vaginal delivery), infant gestational age (weeks), and infant sex (Table 1). The 36 tMBM samples utilized from INSIGHT study participants were chosen at random and researchers were blinded to the maternal infant characteristics. The medical and demographic characteristics of the present cohort are reflective of those present in the larger INSIGHT investigation (15). Term MBM samples were collected during the first 48 hours after delivery (tColostrum; n=13) and 3–4 weeks after delivery (tMBM) both before feeding (foremilk; n=10) and after feeding (hindmilk; n=13). No study has yet identified differences in the miRNA content of foremilk and hindmilk MBM and none were identified in our subset of samples (Supplementary Table S8; online). Therefore foremilk and hindmilk samples were combined as part of the tMBM analysis. MBM samples were manually expressed by participants into 50 ml RNAse free conical collection tubes (ThermoFisher Scientific) after self-cleaning the nipple with soap and water. Samples were stored at −4°C for up to 24 hours prior to collection by research staff and storage at −80°C while awaiting analysis.

Table 1

Medical and demographic characteristics of maternal-infant dyads

	Maternal Characteristics				Infant Characteristics
	Age(years)	Pre-pregnancy BMI(kg/m2)	BMI at time of Delivery(kg/m2)	White Race/Ethnicity(%)	Hyper-tension(%)	C-Section Delivery(%)	Gestational Age(weeks)	Male (%)
Term Group (n=13)	32 ± 4	26 ± 8	34± 12	94	6	24	39 ± 1	47
Preterm Group (n=31)	29 ± 5	N/A	N/A	86	43*	60*	32 ± 3*	55

Abbreviations: N/A = Not available;

denotes p<0.05

Preterm Milk Recruitment and Collection

pMBM was collected as part of an ongoing larger study to identify relationships between genetic variation, diet, and gastrointestinal microbiota in preterm infants. This cohort includes both primiparous and multiparous women, 18–40 years of age giving birth to premature infants (28–37 weeks gestation) that require care in the Neonatal Intensive Care Unit. All breastfeeding mothers involved in the initial study were asked to provide an MBM sample as part of the ongoing microbiota investigation. A subset of those samples was chosen at random for inclusion in the current study by researchers blinded to maternal infant characteristics. Mother infant dyads were not excluded from the pMBM cohort based upon birth weight, or medical complications such as emergent Cesarean delivery. The pMBM samples (n = 31) were all collected 3–4 weeks after delivery by expressing milk from both breasts using an electric Medela Lactina or Symphony breast pump. No distinction was made between fore- and hind-milk collection in this cohort. A 10 ml aliquot of MBM was stored in a sterile container at 4°C for up to 24 hours and then frozen at −80°C until samples could be retrieved for analysis. Basic health and demographic information was obtained for each mother-infant dyad, including maternal age, maternal race/ethnicity, presence of maternal hypertension, delivery method (cesarean section or vaginal), infant gestational age (weeks), and infant sex (Table 1).

RNA Extraction

RNA was extracted with a Norgen Circulating and Exosomal RNA Purification Kit (Norgen Biotek, Ontario, Canada) per manufacturer instructions. Samples were flash-thawed and spun for 20 minutes at 4°C at 800 RPM to separate the lipid and skim fractions. For each sample, 50 μL of lipid and 200μL of skim milk were pipeted into separate 50 mL tubes. When there was not a sufficient amount of lipid or skim, the maximum amount that could be extracted was used. All samples were diluted to 1mL using PBS. A two-step elution procedure was employed using 50 μL of Elution solution centrifuged at 2,000 RPM for 2 minutes and then 3 minutes at 14,000 RPM. This was followed by a second elution using the 50μL of initial flowthrough centrifuged for only 3 minutes at 14,000 RPM. The final RNA concentration was determined on a Nanodrop Spectrophotmeter and the extracted RNA was stored at −80°C prior to sequencing.

RNA sequencing and alignment

Yield and quality of the RNA samples was assessed using an Agilent 2100 Bioanalyzer prior to library construction with the NEXTflex Small RNA-Seq Kit v3 (Bioo Scientific; Austin, Texas). Multiplexed samples were run on an Illumina HiSeq 2500 Instrument at a targeted depth of one million reads per sample. FastQ outputs were clipped, trimmed, and filtered to a maximum read length of 30 base pairs using the FASTX Toolkit Module in Mobaxterm. Reads were aligned to the hg38 build of the human genome in Partek Flow (Partek; St. Louis, Missouri) using the Bowtie 2 algorithm. Total miRNA counts within each sample were quantified with miRBase mature microRNAs v21 and read counts were normalized across samples using a trimmed mean of M-values (TMM) method. Only miRNAs with raw read counts greater than five in at least 50 % of samples were evaluated in the differential expression analysis. A multi-model approach employing gene set analysis (GSA) in Partek Flow was used to identify individual miRNAs with differential expression between tMBM and pMBM samples collected 3–4 weeks post-delivery (in lipid and skim fractions respectively). Differential expression was defined as a Benjamini Hochberg False Discovery Rate (FDR) < 0.05. The data set supporting the results of this article will be available in the NCBI Sequence Read Archive.

Functional annotation

The miRNAs with differential expression across tMBM and pMBM fractions underwent functional annotation analysis in DIANA mirPath v3 online software (http://snf-15788.vm.okeanos.grnet.gr/). The microT CDS algorithm was used to identify mRNA targets for each miRNA of interest. MicroT-CDS uses miRBase v18 and Ensembl v69 for computational prediction of interactions between miRNAs and protein-coding genes in a species-specific manner (16). DIANA mirPath identified KEGG pathways and gene ontology (GO) categories with significant enrichment (FDR < 0.05) for the miRNA targets in each MBM fraction. A list of high confidence mRNA targets from each MBM fraction (microT CDS score > 0.98) was interrogated for protein-protein interaction networks using moderate stringency settings (interaction score > 0.40) in String v10 software (http://string-db.org). String employs BIP, BioGRID, HPRD, IntAct, MINT, and PDP databases to identify biochemical interactions between lists of mRNA-protein products (17). Biological pathways over-represented in high confidence target gene sets were reported. Finally, for the lipid fraction of pMBM (which showed functional enrichment for metabolic processes) a list of 22 genes highly-implicated in human obesity (18) was interrogated for target overlap with pMBM miRNAs of interest using miRDB software (http://www.mirdb.org).

Statistical analysis

Analysis of maternal/infant medical and demographic data with a two-tailed student’s t-test was performed to identify differences between term and preterm groups. A power analysis determined the sample size necessary for differentiating miRNAs among pMBM and tMBM cohorts. Assuming a logarithmic distribution of miRNA expression, it was determined that 22 MBM samples per group provided 80% power to detect 1.8-fold differences among 500 miRNAs (per-gene alpha =0.002). Based upon these a priori calculations we collected 23 tMBM and 31 pMBM samples and examined differences across 465 miRNAs in lipid MBM fractions. Individual miRNAs with significant differences between tMBM and pMBM fractions were identified using GSA analysis with Benjamini Hochberg FDR correction. Those miRNAs with FDR < 0.05 were considered to have significant changes between groups. Fold change in miRNA levels between groups were reported as log2 values. The mRNA targets and their biological pathways were analyzed for functional enrichment with Fisher’s Exact Tests, using a 2×2 table to compare the ratio of identified mRNA targets to the number of targets expected by chance. A non-parametric Kruskal-Wallis test was used to examine differences in individual miRNAs across lipid fractions of tMBM, pMBM, and tColostrum. Low volumes of colostrum produced following delivery do not allow for separation of skim and lipid fractions. Thus, tColostrum samples were compared only with lipid fractions of tMBM and pMBM (given their similar high-fat content). This comparison provided insight into how MBM miRNA might evolve following delivery (tColostrum versus tMBM at one month) while also investigating the influence of maternal and infant characteristics on MBM miRNA profiles. Spatial relationships between total miRNA profiles for each group were examined with a two-dimensional partial least squares discriminant analysis (PLS-DA) in MetaboAnalyst software. Hierarchical clustering was performed for the three groups using a Pearson distance measure of the top 20 miRNAs on Kruskal-Wallis testing with a complete clustering algorithm. Finally, the relationship between maternal/infant medical and demographic characteristics and MBM miRNAs of interest were examined with Pearson correlation coefficients.

RESULTS

Medical and demographic characteristics

There was no significant difference in maternal age (p=0.2), maternal race/ethnicity (% white; p=0.4) or infant sex (p=0.6) between the pMBM and tMBM groups. Within the pMBM group, there was a higher percentage of mothers with gestational hypertension (p=0.0008) and a higher percentage of Cesarean deliveries (p=0.008). Gestational age at the time of delivery was significantly lower in the pMBM group (p= 4E-19) (Table 1).

Comparison of tMBM and pMBM Lipid Fractions

Of the 2588 mature human miRNAs interrogated for expression levels in pMBM and tMBM, 465 were measured in the lipid fractions (Supplementary Table S1; online). There were 113 miRNAs with significant expression differences (FDR < 0.05) between tMBM and pMBM lipid samples. Sixty eight of the miRNAs were down regulated in pMBM lipid fractions and 45 were-up regulated. Eighty-two of these miRNAs have been described in previous studies of MBM (6, 14). The 15 miRNAs with the most significant differences across tMBM and pMBM lipid fractions (Table 2) underwent functional analysis. Of these, two pairs of miRNAs have overlapping seed sequences (miR-1260a/miR-1260b and miR378a-3p/miR378c) and six miRNAs have been previously described in breast tissue, breast milk, or breast discharge (19–24).

Table 2

The 15 miRNAs most “altered” in the lipid fraction of pMBM relative to lipid tMBM one month post-delivery

Top 15 miRNAs(lipid fraction)	FDR	Fold change(pMBM)	Seed sequence	Overlapping miRNA Seeds	Target mRNAs(Total/High- confidence)	Breast-Related Expression
hsa-miR-4470	4.31E-06	6.75E+01	GGCAAAC	None	633/63
hsa-miR-4687-3p	4.31E-06	−3.54E+00	GGCUGUU	None	428/46	Nipple discharge from intraductal papillomaa
hsa-miR-1260a	9.25E-06	−2.10E+00	UCCCACC	hsa-miR-1260b	39013
hsa-miR-1260b	1.70E-05	−2.01E+00	UCCCACC	hsa-miR-1260a	73/2	MCF-7 breast cancer cellsb
hsa-miR-378a-3p	1.70E-05	−1.97E+00	CUGGACU	hsa-miR-378b/c/d/e/f/h/I hsa- miR-422a	231/9	Breast milk, MCF-7 breast cancer cellsc
hsa-miR-378g	3.35E-05	−2.10E+00	CUGGGCU	None	405/16
hsa-miR-378c	5.36E-05	−1.88E+00	CUGGACU	hsa-miR-378a/b/d/e/f/h/I hsa- miR-422a	255/12
hsa-miR-4474-5p	5.94E-05	3.33E+01	UAGUCUC	None	227/14	Plasma in unexplained recurrent spontaneous abortiond
hsa-miR-6763-3p	5.94E-05	−2.19E+00	UCCCCGG	None	17/1
hsa-miR-5585-5p	7.16E-05	−1.70E+00	GAAGUAC	None	342/18
hsa-miR-1200	3.39E-04	−1.64E+00	UCCUGAG	None	597/26	MCF-7 breast cancer cellse
hsa-miR-4783-5p	3.44E-04	−1.99E+00	GCGCGCC	None	0/0
hsa-miR-6772-3p	4.39E-04	4.97E+00	UGCUCCU	None	88/7
hsa-miR-187-5p	5.25E-04	−1.78E+00	GCUACAA	None	173/9	BRCA1 triple-negative breast cancer cellsf
hsa-miR-6510-5p	5.25E-04	−1.99E+00	AGCAGGG	None	542/17

Target mRNAs identified with DIANA-miRPATH software. High-confidence targets included those with target score > 0.98. a. Zhang et al., 2015, b. Gonul et al., 2015, c. Ikeda et al., 2015, d. Qin et al., 2016, e. Shah et al., 2011, f. Matamala et al.,2016; Abbreviations: pMBM = premature maternal breast milk (n=31), tMBM = term maternal breast milk (n=23).

The 15 miRNAs most altered in the lipid fraction of pMBM had a total of 4401 predicted mRNA targets in DIANA mirPath v3. These mRNA targets were associated with 90 GO categories and 24 KEGG pathways (Supplementary Table S2; online). Notable GO categories with targeted enrichment included cellular nitrogen compound metabolic process, biosynthetic process, catabolic process, symbiosis, and viral process. Filtering down to the top 247 high-confidence mRNA targets with target scores ≥ 0.980 revealed 15 mRNAs targeted by more than one miRNA on the list. Analysis of the 247 mRNAs in String v10 revealed a significant protein protein interaction network (p= 8.1E-4) containing 230 nodes and 148 edges with a clustering coefficient of 0.752 (Supplementary Figure S1; online). Analysis of this network identified 13 biologic processes with enrichment (Table 3). Of note, three of the top five GO processes were involved in negative regulation of cellular macromolecule biosynthesis and each included more than 30 genes. In addition, an investigation of 22 unique genes implicated in human obesity (18) revealed that 19 genes had predicted miRNA interactions, and 13 of these were targeted by a miRNA of interest in lipid pMBM. For two obesity related genes (Adrenergic, beta-receptor 3, ARDB3; and Nuclear receptor sub-family 3, group C, member 1, NR3C1) there was consistent up-regulation of all target miRNAs present in lipid pMBM. For three other obesity-related genes (Leptin, LEP; Tumor necrosis factor, TNF; and Uncoupling protein 3, UCP3) there was consistent down-regulation of all target miRNAs present in lipid pMBM.

Table 3

The miRNAs differentially expressed in lipid pMBM target gene ontology (GO) pathways related to macromolecule and cellular biosynthetic metabolism

pathway ID	pathway description	Gene count	FDR	Target mRNAs with protein products in network
GO.0045892	negative regulation of transcription, DNA-templated	31	0.0046	ARID5A,BMI1,CEBPG,CTBP1,CYP1B1,EGLN1,EHMT1,FBXW11,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,TFCP2L1,TSG101,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.2000113	negative regulation of cellular macromolecule biosynthetic process	34	0.0046	ARID5A,BMI1,CEBPG,CELF1,CTBP1,CYP1B1,EGLN1,EHMT1,FBXW11,FBXW7,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,METTL3,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,TFCP2L1,TSG101,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0010629	negative regulation of gene expression	35	0.0069	ARID5A,BMI1,CEBPG,CELF1,CTBP1,CYP1B1,DYRK1A,EGLN1,EHMT1,FBXW11,FLOT2,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,METTL3,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,TFCP2L1,TSG101,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0010558	negative regulation of macromolecule biosynthetic process	32	0.0103	ARID5A,BMI1,CEBPG,CELF1,CTBP1,CYP1B1,EGLN1,EHMT1,FBXW11,FBXW7,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,METTL3,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0031327	negative regulation of cellular biosynthetic process	33	0.0103	ADRA2C,ARID5A,BMI1,CEBPG,CELF1,CTBP1,CYP1B1,EGLN1,EHMT1,FBXW11,FBXW7,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,METTL3,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0045934	negative regulation of nucleobase-containing compound metabolic process	32	0.0103	ADRA2C,ARID5A,BMI1,CEBPG,CTBP1,CYP1B1,DYRK1A,EGLN1,EHMT1,FBXW11,FBXW7,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0050794	regulation of cellular process	122	0.0103	ABI2,ADCYAP1,ADRA2C,ALMS1,ATAD2,ATF6B,AZIN1,BCL2L2,BMI1,BMPR1A,C15orf23,CACNA2D4,CAMSAP2,CCNY,CD47,CDC42BPG,CEBPG,CELF1,CENPB,COBL,COLQ,CPLX1,CTBP1,CYP1B1,DACT2,DSCAM,DYRK1A,EGLN1,EHMT1,EIF3D,ELL2,ESPL1,ETV1,FBXW7,FLOT2,FNDC5,FOXO1,GCC2,GDI1,GLS2,GPR110,H3F3B,HIPK3,HNRNPH1,IL36RN,ITSN2,JAZF1,KCNA4,KDM6A,KHDRBS2,KPNA1,LECT1,LTF,MAP1A,MAP3K12,METTL3,MKX,MYO18A,MYO1C,MYRIP,NAB1,NAMPT,NR2F2,NREP,OSTF1,PABPN1,PAX5,PDE10A,PDE7A,PHACTR4,PHF21A,PITPNC1,PKD1,PLCZ1,POLB,POU2F1,PPM1L,PPP1R3F,PPP3CB,PRKAR2A,PTPRG,PTX3,PXDN,RAB3B,RAB6C,RABL3,RALBP1,RC3H2,RFWD2,RIPK2,RPAP2,SEC62,SEMA3A,SESN3,SLC2A8,SLC30A1,SLITRK5,SNCAIP,SNX27,SRF,STK4,SYNJ1,TACC1,TDG,TEX15,TFAP2A,TIRAP3,TMEM64,TMEM67,TMX2,UBAC2,UBN2,UNC13C,USP31,WWC2,YY1AP1,ZBTB44,ZEB2,ZNF217,ZNF286B,ZNF451,ZNF675
GO.0051172	negative regulation of nitrogen compound metabolic process	34	0.0103	ADRA2C,ARID5A,BMI1,CEBPG,CELF1,CTBP1,CYP1B1,DYRK1A,EGLN1,EHMT1,FBXW11,FBXW7,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,METTL3,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0051253	negative regulation of RNA metabolic process	30	0.0103	ARID5A,BMI1,CEBPG,CTBP1,CYP1B1,DYRK1A,EGLN1,EHMT1,FBXW11,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.1903507	negative regulation of nucleic acid-templated transcription	29	0.0103	ARID5A,BMI1,CEBPG,CTBP1,CYP1B1,EGLN1,EHMT1,FBXW11,FOXF2,FOXG1,H3F3B,IKZF4,JAZF1,MECOM,MKX,NAB1,NR2F2,PAX5,PHF21A,POU2F1,TDG,TFAP2A,UBE2D2,UBE2D3,WWC2,ZEB2,ZFHX3,ZNF217,ZNF675
GO.0050789	regulation of biological process	124	0.0181	ABI2,ADCYAP1,ADRA2C,ALMS1,APC,ATAD2,ATF6B,ATG12,AZIN1,BCL2L2,BMI1,BMPR1A,C15orf23,CACNA2D4,CAMSAP2,CCNY,CD47,CDC42BPG,CEBPG,CELF2,CENPB,COBL,COLQ,CPLX1,CTBP1,CYP1B1,DACT2,DBF4,DSCAM,DYRK1A,EGLN1,EHMT1,EIF3D,ELL2,ENPEP,ESPL1,ETV1,FBXW7,FLOT2,FNDC5,FOXO1,GCC2,GCH1,GDI1,GLS2,GPR110,H3F3B,HIPK3,HNRNPH1,ITSN2,JAZF1,KCNA4,KDM6A,KHDRBS2,KPNA1,LIPC,LTF,MAP1A,MAP3K12,METTL3,MKX,MYO18A,MYO1C,MYRIP,NAMPT,NR2F2,NREP,OSTF1,PABPN1,PAX5,PDE10A,PDE7A,PHACTR4,PHF21A,PITPNC1,PKD1,PLCZ1,POLB,POU2F1,PPM1L,PPP1R3F,PPP3CB,PRKAR2A,PTPRG,PTX3,PXDN,RAB3B,RAB6C,RABL3,RALBP1,RC3H2,RFWD2,RIPK2,RPAP2,SEC62,SEMA3A,SESN3,SLC2A8,SLC30A1,SLITRK5,SNCAIP,SNX27,SRF,STK4,SYNJ1,TACC1,TEX15,TFAP2A,TMEM64,TMEM67,TMX2,UBAC2,UBN2,UBXN2B,UNC13C,USP31,WWC2,YY1AP1,ZBTB44,ZEB2,ZNF217,ZNF286B,ZNF451,ZNF675
GO.0048513	organ development	49	0.0261	AARD,APC,APLN,ARID5A,BCL2L2,BMI1,BMPR1A,CBFA2T3,CEBPG,COBL,CYP1B1,DSCAM,EGLN1,ENPEP,ETV1,FBXW7,FOXG1,FOXO1,FRS2,H3F3B,IMPAD1,KDM6A,LECT1,LTF,MAGI2,MBNL1,MEGF11,MKX,NAB1,NR2F2,PGF,PHACTR4,PKD1,POLB,POU2F1,PPP3CB,PRRC2C,RC3H2,SEMA3A,SLITRK5,SRF,STK4,TEX15,TFAP2A,TFCP2L1,TMEM143,TSG101,WDR62,ZFHX3
GO.0060341	regulation of cellular localization	28	0.0377	ADRA2C,APC,APLN,BMPR1A,C15orf23,CPLX1,FBXW11,FBXW7,GCC2,GDI1,GLS2,IL36RN,KPNA1,MAGI2,MYO18A,MYO1C,MYRIP,NPSR1,NUMB,PKD1,PPP3CB,PRKAR2A,RAB3B,SLC30A1,SNCAIP,SYNJ1,TSG101,UBAC2

Abbreviations: FDR = Benjamini Hochberg False Discovery Rate, pMBM = premature maternal breast milk

Comparison of tMBM and pMBM Skim Milk Fractions

There were fewer miRNAs identified in skim milk fractions. Thirty-two miRNAs showed robust expression in the skim milk fractions of both tMBM and pMBM (Supplementary Table S3; online). All 32 of the miRNAs measured in skim milk fractions were also found in lipid fractions. In the skim fractions, there were 12 miRNAs with significant expression differences between pMBM and tMBM samples (Table 4). Seven of the 12 miRNAs were down-regulated in skim fractions of pMBM and five were upregulated. Three of these miRNAs have been described in previous studies of MBM and nine are novel (6, 9–10, 14). Of these 12 miRNAs, two pairs have overlapping seed sequences (miR-1260a/miR-1260b and miR378a-3p/miR378c), and five have been previously described in breast tissue or breast milk (19–20, 25–27).

Table 4

The 12 miRNAs “altered” in the skim fraction of pMBM relative to skim tMBM one month post-delivery

Top 12 miRNAs(skim fraction)	FDR	Fold change(pMBM)	Seed sequence	Overlapping miRNA Seeds	Target mRNAs(Total/High-confidence)	Breast-Related Expression
hsa-miR-378g	1.58E-04	−2.12E+00	CUGGGCU		405/16
hsa-miR-378a-3p	6.40E-04	−1.90E+00	CUGGACU	hsa-miR-378b/c/d/e/f/h/Ihsa-miR-422a	231/9	Breast milk, MCF-7 breast cancer cellsa
hsa-miR-5787	6.40E-04	2.50E+00	GGCUGGG	hsa-miR-4505	728/29
hsa-miR-1260a	7.47E-04	−1.77E+00	UCCCACC	hsa-miR-1260b	390/13
hsa-miR-1260b	7.84E-04	−1.66E+00	UCCCACC	hsa-miR-1260a	73/2	MCF-7 breast cancer cellsb
hsa-miR-378c	7.84E-04	−1.69E+00	CUGGACU	hsa-miR-378b/c/d/e/f/h/I hsa-miR-422a	255/12
hsa-miR-7975	7.84E-04	−1.67E+00	UCCUAGU		50/0
hsa-miR-7704	9.03E-04	2.12E+00	GGGGUCG		13/1
hsa-miR-4784	1.58E-02	1.76E+00	GAGGAGA	hsa-miR-3150b-3p	482/39	BRCA1 triple negative breast cancer cellsc
hsa-miR-4294	1.69E-02	1.69E+00	GGAGUCU		114/1	Breast cancer cellsd
hsa-miR-4783-3p	3.11E-02	2.24E+00	CCCGGUG		21/2
hsa-miR-4454	4.97E-02	−1.39E+00	GAUCCGA		5/0	Inflammatory breast cancere

Target mRNAs identified with DIANA-miRPATH software. High-confidence targets included those with target score > 0.98. a. Ikeda et al., 2015, b. Gonul et al., 2015, c. Yang et al., 2015, d. Kumari, 2014, e. Maltseva et al., 2014; Abbreviations: pMBM = premature maternal breast milk (n=31), tMBM = term maternal breast milk (n=23).

The 12 miRNAs with significant differences in expression in pMBM had a total of 2767 predicted mRNA targets in DIANA mirPath v3. These mRNAs showed enrichment for 10 GO categories and 9 KEGG pathways (Supplementary Table S4; online). The most significantly enriched GO class was molecules involved in nitrogen metabolism. The top three KEGG pathways included genes involved in glycosphingolipid biosynthesis, lysine degradation, and calcium signaling. There were 124 mRNA targets that survived the high confidence filter and 10 of these were targeted by more than one miRNA on the list. Functional analysis of the 124 mRNAs in String v10 demonstrated a significant protein protein interaction network (p=0.034) containing 111 nodes and 23 edges with a clustering coefficient of 0.892 (Supplementary Figure S2; online), however GO analysis of this network failed to identify a protein-protein interaction pathway with significant enrichment.

Combined analysis of pMBM Skim and Lipid Fractions

When comparing both lipid and skim fractions of pMBM with tMBM fractions, we identified nine differentially expressed miRNAs (miR-378a 3p, miR378c, miR-378g, miR-1260a, miR-1260b, miR-4783-5p, miR-4784, miR-5787, miR-7975). All nine miRNAs were altered in the same direction in skim and lipid fractions of pMBM with the exception of miR-4783-5p (which was up-regulated in the skim pMBM fraction but down regulated in the lipid pMBM fraction). Two miRNAs were up regulated in both fractions of pMBM (miR-5787 and miR-4784) and six were down-regulated. Functional analysis of the nine miRNAs with differential expression patterns in both the lipid and skim fractions of pMBM yielded 2614 predicted mRNA targets in DIANA mirPATH v3. These mRNAs showed enrichment for 42 GO categories and 5 KEGG pathways (Supplementary Table S5; online). Top GO and KEGG targets included networks of mRNAs related to micronutrient metabolism, such as biosynthetic processes, cellular nitrogen metabolism, glycosphingolipid biosynthesis, and lysine degradation. Eight of the nine miRNAs also targeted 62 genes involved in viral process.

Comparison of pMBM with tMBM and tColostrum

A Kruskal Wallis test was used to examine differences in individual miRNAs between lipid tMBM (collected 3–4 weeks post-delivery), lipid pMBM (collected 3–4 weeks post-delivery), and tColostrum samples (collected within 48 hours of delivery). There were 127 miRNAs with significant differences among these groups (Supplementary Table S6; online). For most of the miRNAs interrogated, expression patterns in tMBM and tColostrum trended together in the opposite direction of pMBM (Supplementary Figure S3; online). PLS-DA analysis achieved partial separation of the three sample types based upon total miRNA profile (Figure 1). A hierarchical clustering analysis using the 20 miRNAs with the most significant changes across the three groups demonstrated complete separation of pMBM samples from both tMBM and tColostrum samples (Figure 2). Note that duplicate (foremilk and hindmilk) tMBM samples (e.g. tMBM 1a and tMBM 1b) tended to cluster together. In addition, tColostrum and tMBM samples from the same mother (e.g. tColostrum 7 and tMBM 7a) often cluster near one another.

Figure 1

Partial Least Squares Discriminant Analysis.

A PLS DA of the total microRNA profile for pMBM (class 2, +, n=31), tMBM (class 3, ×, n=23) and tColostrum (class 0, ∆, n=10) achieved partial separation using two dimensions and accounted for 23.6% of the variance between samples.

Figure 2

Hierarchical clustering (HC) analysis.

HC of the 20 miRNAs with the most significant changes across the three groups showed clustering of pMBM (class 2; n=31) from tMBM (class 3; n=23) and tColostrum (class 0; n=13). Note that tMBM fore milk (n=10) and hind milk (n=13) samples (denoted a and b respectively), as well as tMBM and tColostrum samples (taken from the same mother one month apart) were also spatially clustered. Gray scale values indicate average Z score of normalized abundance for each miRNA.

Relationship of medical/demographic characteristics and MBM miRNA

A Pearson correlation analysis was used to explore the relationships between medical/demographic variables and miRNA expression. This analysis identified significant correlations between gestational age (measured as a continuous variable) and 21 of the 26 miRNAs “altered” in lipid or skim fractions of pMBM samples (Supplementary Table S7; online). Six of the 26 miRNAs were correlated with method of delivery, but none correlated with maternal race/ethnicity, presence of maternal hypertension, twin gestation, or infant sex. Of the six miRNAs correlated with delivery method, three were positively correlated, three were negatively correlated, and all six showed inverse patterns compared to gestational age. Heatmap visualization of Pearson’s correlations demonstrated hierarchical clustering of gestational age with lipid and skim fractions of miR-1260a/b as well as the lipid fraction of miR-6763-3p and the skim fraction of miR-7975 (Figure 3). All these factors were positively correlated.

Figure 3

Pearson correlations.

Heatmap visualization of Pearson’s correlation for the 26 miRNAs altered in lipid or skim pMBM demonstrates correlations between gestational age (as a continuous variable) and 21 miRNAs. Six of these miRNAs were correlated with delivery method, but no correlation was found between miRNAs and other variables. Color represents Pearson R values.

DISCUSSION

The expression of miRNAs in the lipid and skim milk fractions of pMBM differs significantly from tMBM fractions one month after delivery. There are nine miRNAs “altered” across both fractions and these miRNAs target a number of transcripts involved in metabolic processes. The pathway with the most significant enrichment in miRNA targets from pMBM is glycosphingolipid biosynthesis. Glycosphingolipids are important for neurodevelopment (28), are essential to a functioning cell membrane, and play a major role in signal transduction of lipid rafts (29). These same lipid domains may aid in the protected delivery of miRNAs to the infant gut, regulating gastrointestinal function, metabolism and energy utilization.

Potential mechanisms for a unique miRNA profile in pMBM

Differences in miRNA composition of pMBM and tMBM may result from regulation within the mammary epithelial cell nucleus (through pri miRNA expression or pre miRNA processing), or they may represent a response to environmental changes (such as abrupt premature delivery) which alters the way miRNAs are packaged and extruded into MBM. MiRNAs are packaged in a variety of ways, including vesicles, shedding, exosomal transfer, and RICS-complex protein-binding (30). Because individual miRNAs show preference for specific packaging mechanisms (31) a change in the ratio of carriers (e.g. exosomes or nucleoproteins) could affect the secretion of certain miRNAs in pMBM. Given that tMBM and pMBM differ in macronutrient and micronutrient composition (1–2), its plausible that they might also differ in types of miRNA carriers.

Alternatively, premature delivery might influence production of miRNAs within the cell nucleus. Given the milieu of hormonal changes occurring pre and post-partum, a shift in maternal hormones could influence miRNA transcription. For example, lactogenic hormones alter expression and secretion of miRNAs in cultured mammary cells (7). Mothers of preterm infants have lower prolactin levels (32) and this may alter expression and secretion of miRNAs. The miRNA profile in breast cells is altered by exposure to estrogen (33), and exposure to progesterone increases the levels of miRNA processing machinery, which could result in an altered miRNA profile (34). Since estrogen and progesterone levels increase throughout pregnancy, mothers of premature babies have lower circulating levels of both hormones at delivery. This could have a lasting impact on the miRNAs in pMBM, providing potential evolutionary benefits for the premature neonate.

Implications for infant health

Premature infants face increased risk for a number of medical co morbidities, including failure to thrive, sepsis, necrotizing enterocolitis (NEC), and neurodevelopmental delays. These risks may be potentiated or mitigated by maternal-infant transfer of MBM miRNAs. Notably, premature infants have a vastly different growth trajectory than term infants and metabolic pathways, such as the glycolytic pathway, are affected by miRNAs “altered” in pMBM (Table 3, Supplementary Table S2; online). An investigation of miRNAs that target obesity-related genes (18) reveals changes relevant to growth and metabolism in the premature newborn. Interestingly, the miRNAs that target UCP3, LEP, and TNF are down regulated in lipid pMBM. UCP3 alters metabolism in response to fuel depletion, while TNF increases LEP production and adipocyte insulin resistance (35). Together, polymorphisms in UCP3 and TNF act synergistically to increase risk for type 2 diabetes mellitus (36). Thus, widespread down regulation of miRNAs targeting these genes might influence glucose homeostasis in the premature infant. Conversely, up-regulation of miRNAs in lipid pMBM that target ADRB3 and NR3C1 could lead to down-regulation of their protein products in the premature infant. A decrease in ADRB3 and NR3C1, which are both negative regulators of adipogenesis (37–38), could increase fat production have functional benefits in the cachectic premature newborn.

Premature infants face increased risk of viral and bacterial infections. Previous studies have described enrichment of immune-related miRNAs in MBM (5, 10) and the current investigation identifies a pair of immune-related miRNAs (miR-1260a and miR-1260b) that are down regulated in pMBM compared with tMBM (Tables 2, 4). Human activated B cell lines demonstrate robust miR-1260 expression in comparison to indolent germinal center cells, suggesting miR-1260 may influence B-cell proliferation (39). Thus, down regulation of miR 1260 in pMBM may not serve a protective function in premature infants, but could contribute to their immune susceptible state. Coupled with decreased miRNA-targeting of TNF in pMBM these changes could have pathophysiologic implications in NEC. Conversely, robust expression of miR-1260a/miR-1260b in tMBM may impart immune protection on the infant gut and play a role in the decreased incidence of gastroenteritis in breast fed infants (40).

The relative influence of prematurity on MBM miRNA

To account for the vast differences between a 35-week and a 27-week pre-term infant, the current study investigated relationships between MBM miRNA and gestational age as a continuous variable. This approach reaffirmed that 21 of the 26 miRNAs altered in pMBM versus tMBM samples were significantly related to premature delivery. Furthermore, none of the 26 miRNAs influenced by gestational age were correlated with maternal race/ethnicity, maternal hypertension, twin gestation, or infant sex (Figure 3). Six of these miRNAs were correlated with cesarean section delivery and this may represent an important health factor contributing to MBM miRNA shifts following premature delivery. These results suggest that maternal factors may influence MBM miRNA content. This is consistent with a previous study (13), showing that expression of miRNAs in MBM was perturbed by a maternal high fat diet. In fact, the miRNA composition of MBM may be influenced more by maternal factors than length of lactation (i.e. time since delivery). This idea is supported by Figure 2, showing that miRNA profiles for tMBM preferentially cluster alongside tColostrum from maternal pairs, rather than pMBM samples taken the same duration post-delivery. Prematurity delivery, acting as a health factor, may influence the miRNA profile of MBM for a substantial period of lactation. This has implications for infant growth trajectory throughout the breast-feeding period. It may explain why premature infants grow more appropriately on their mother’s own milk compared to donor milk from mothers delivering at term.

Limitations

The current study is limited by utilization of pMBM and tMBM samples collected from two separate investigations. This approach introduces differences in MBM collection methods (e.g. manual versus machine-pump) which could influence miRNA content. MBM miRNA arises primarily from the mammary gland epithelium rather than maternal circulation (9) and the compressive nature of manual pumping could disrupt this epithelial layer, resulting in the collection of intracellular miRNA. To control for this, we have performed differential centrifugation to eliminate epithelial components in both pMBM and tMBM samples, focusing the analysis on extracellular miRNAs in lipid and skim milk layers. Furthermore, the miRNAs in the present analysis included only those miRNAs present in >50% of samples. Thus, intracellular miRNAs introduced by a manual pumping technique should effectively be excluded from consideration.

An additional limitation of this study is its reliance on both foremilk and hindmilk samples within the tMBM cohort. Hindmilk, which is enriched in lipid content, but not lactose or protein, could demonstrate increased miRNA diversity through the stability of lipid-bound miRNAs. However, raw RNA read counts in hindmilk tMBM samples did not exceed those observed in foremilk samples. In addition, such differences should be eliminated by differential centrifugation and analysis of skim fractions, which removes lipid bound miRNAs. One previous study examining differences in pre and post-feed miRNAs within human milk cells identified 33 miRNAs with “alterations” post-feed (11). Only three of these miRNAs (miR-1260b, miR-4454, and miR-671-3p) show differences in our comparison of pMBM and tMBM. This small overlap is likely related to removal of human milk cells by centrifugation. This technique focuses on MBM fractions with the highest likelihood for maternal-infant transfer.

A final limitation of the present investigation is its sample size. Previous studies have utilized 4–10 MBM samples for exploring miRNA profiles, but our a priori power analysis required a minimum of 22 samples per group. To power a comparison of tMBM and pMBM groups we included 23 and 31 samples from these cohorts, respectively. Thus, our investigation is well powered for an analysis of 465 miRNAs across term and preterm cohorts, but is somewhat underpowered for comparison of tColostrum samples (n=13) and comparisons between foremilk (n=10) and hindmilk (n=13) samples. Thus, these findings should be interpreted with caution.

Conclusion

This study has identified nine miRNAs with differences in lipid and skim fractions of pMBM compared with tMBM collected 3–4 weeks post-delivery. Many of the genes targeted by these miRNAs influence metabolic pathways and could affect infant growth trajectory via maternal-fetal transfer. A subset of miRNAs was correlated with delivery method and this may represent one medical factor which triggers changes in the MBM miRNA profile. These findings underscore how MBM, which is generally considered the optimum nutrition for infants, may be individually adapted to specific infant needs through the regulation of factors such as miRNA.