Placental microRNA expression associates with birthweight through control of adipokines: results from two independent cohorts

ABSTRACT

MicroRNAs are non-coding RNAs that regulate gene expression post-transcriptionally. In the placenta, the master regulator of foetal growth and development, microRNAs shape the basic processes of trophoblast biology and specific microRNA have been associated with foetal growth. To comprehensively assess the role of microRNAs in placental function and foetal development, we have performed small RNA sequencing to profile placental microRNAs from two independent mother-infant cohorts: the Rhode Island Child Health Study (n = 225) and the New Hampshire Birth Cohort Study (n = 317). We modelled microRNA counts on infant birthweight percentile (BWP) in each cohort, while accounting for race, sex, parity, and technical factors, using negative binomial generalized linear models. We identified microRNAs that were differentially expressed (DEmiRs) with BWP at false discovery rate (FDR) less than 0.05 in both cohorts. hsa-miR-532-5p (miR-532) was positively associated with BWP in both cohorts. By integrating parallel whole transcriptome and small RNA sequencing in the RICHS cohort, we identified putative targets of miR-532. These targets are enriched for pathways involved in adipogenesis, adipocytokine signalling, energy metabolism, and hypoxia response, and included Leptin, which we further demonstrated to have a decreasing expression with increasing BWP, particularly in male infants. Overall, we have shown a robust and reproducible association of miR-532 with BWP, which could influence BWP through regulation of adipocytokines Leptin and Adiponectin.

KEYWORDS:

• placental microRNA
• placenta
• microRNA
• birthweight
• hsa-miR-532-5p
• adipokines
• leptin
• adiponectin

Introduction

MicroRNAs are 21–25 base pair non-coding RNAs that regulate gene expression post-transcriptionally. In humans, microRNAs, as part of the RNA-induced silencing complex (miRISC), target mRNA transcripts based on sequence complementarity to the 3ʹ untranslated region of the mRNA. The miRISC complex inhibits the translation of target mRNAs through one of three mechanisms: 1) through direct cleavage of the mRNA, 2) recruiting mRNA degradation machinery, or 3) blocking translational machinery. MicroRNAs are thought to regulate more than 50% of human genes [1–4].

The placenta is the master regulator of early development. During implantation, the trophectoderm, which gives rise to the placenta, invades through the extracellular matrix of the decidualized endometrium into the uterine wall [5–7]. Implantation is followed by placentation, in which the trophectoderm differentiates and proliferates into a network of branching villi – the functional units of the placenta which facilitate the exchange of gases, nutrients, and waste for the growing embryo [8,9]. Poor placental function is the most common cause of intrauterine growth restriction, which is associated with perinatal morbidity, mortality, and long-term health outcomes. Previous research, primarily conducted in cell lines and model systems, suggests placental microRNAs are pivotal to placental development and function by targeting genes that control trophoblast proliferation and differentiation, apoptosis, invasion, cellular metabolism, as well as vasculo- and angiogenesis (reviewed in [4]).

A number of studies have investigated the role of microRNA expression in placental function and birth outcomes [10–29]. Several have studied links between microRNA expression and low birthweight or intrauterine growth restriction (IUGR) [21–29]. Some of these studies have performed targeted analyses of candidate microRNAs, using single [27] or multiplex qPCR [24]. Many used microarray technology to target a larger number of microRNAs [15,23,25,26]. Most of these studies had sample sizes below 100 [21,22,24–28]. Although some of these studies divided their own cohorts into discovery and validation subsets [22,23,25,26], none of them replicated their findings in independently collected cohorts. From these studies, it is clear that microRNA plays an essential role in placental biology and function; and that adverse pregnancy outcomes are associated with differential expression of various microRNAs. However, well-powered, data-driven and reproducible epidemiologic studies investigating the association of placental microRNAs and birthweight are lacking.

To assess the role of all expressed microRNAs in placental function and foetal health at term, we have utilized small RNA sequencing to profile placental microRNAs from two independent mother-infant cohorts: the New Hampshire Birth Cohort Study (NHBCS; n = 317) and Rhode Island Child Health Study (RICHS; n = 225). We examined the relationships between microRNAs and birthweight percentiles (BWP) independently in RICHS and NHBCS, and evaluated consistency across these cohorts, then used integrated parallel -omics data to characterize a microRNA that associates with BWP at term, as well as its putative mRNA targets.

Methods

Cohorts

The New Hampshire Birth Cohort Study (NHBCS) was initiated in 2009 and is an ongoing study comprised of a cohort of mother-infant pairs. Pregnant women between 18 and 45 years of age were recruited from the study’s participating prenatal care clinics in New Hampshire, USA. Women were included in the cohort if their primary source of drinking water was from an unregulated residential well, they had resided in the same household since their last menstrual period and had no plans to move before delivery. All participants provided written informed consent in accordance with the requirements of the Committee for the Protection of Human Subjects, the Institutional Review Board (IRB) of Dartmouth College. In this study, NHBCS participants were singleton pregnancies recruited between February 2012 and September 2013. Placental microRNA transcript abundance data were provided for this cohort (n = 322). Questionnaires and medical record abstraction were utilized to collect sociodemographic, lifestyle, and anthropometric data (Table 1). For brevity, only variables that were considered in this analysis were included in Table 1.

Table 1. Cohort demographics

	NHBCS (n = 317)	RICHS (n = 225)
Maternal characteristics
Primigravida – % (n)	27% (85)	24% (54)
Pre-pregnancy BMI – mean (range)^a	26 (17–47)	26 (16–47)
Gestational weight gain category^a
Low ^b – % (n)	7% (11)	12% (27)
Normal ^c – % (n)	30% (47)	36% (82)
High ^d – % (n)	63% (99)	52% (116)
Infant characteristics
Race (White) – % (n)	100% (317)	73% (165)
Gestational Age (wks) – mean (range)	40 (31–42)	39 (37–41)
Birthweight (g) – mean (range)	3440 (1380–5216)	3540 (2030–5465)
Birthweight Group
SGA^e – % (n)	4% (14)	16% (36)
AGA^f – % (n)	86% (273)	56% (127)
LGA^g – % (n)	10% (30)	28% (62)
Sex (Female) – % (n)	50% (159)	49% (111)

^aNHBCS statistics for pre-pregnancy BMI and gestational weight gain were based on the 157 samples with available measures.

^bLow – Gestational weight gain below recommendation for pre-pregnancy BMI category.

^cNormal – Gestational weight gain within recommendation for pre-pregnancy BMI category.

^dHigh – Gestational weight gain above recommendation for pre-pregnancy BMI category.

^eSGA – Infants born with a birthweight percentile ≤ 10.

^fAGA – Infants born with a birthweight between the 10^th and 90^th percentiles.

^gLGA – Infants born with a birthweight percentile ≥ 90.

The Rhode Island Child Health Study (RICHS) enrolled mother-infant pairs from the Women & Infants Hospital in Providence, Rhode Island, between September 2010 and February 2013. Mothers were at least 18 years of age, had no life-threatening conditions, and delivered singletons free of congenital/chromosomal abnormalities at or after 37 weeks of gestation. Infants who were born small for gestational age (≤10th BWP) or large for gestational age (≥90th BWP) were oversampled. Infants who were adequate for gestational age (between the 10th and 90th BWPs), and matched for gestational age and maternal age were coincidentally enrolled. All participants provided written informed consent and all protocols were approved by the IRBs at the Women & Infants Hospital of Rhode Island and Emory University, respectively. Data provided by this study include placental microRNA transcript abundance (n = 230) and total RNA transcript abundance (RNA-seq, n = 199). Interviewer-administered questionnaires were utilized to collect sociodemographic and lifestyle data. Structured medical record review was used to collect anthropometric and medical history data. For brevity, only variables that were considered in this analysis were included in Table 1.

Data collection

Anthropomorphic measures

Both of the studies have parallel demographic and anthropomorphic measures for mothers and newborns, as well as placental microRNA transcript abundance. Birthweight (g), head circumference (cm), and birth length (cm) were abstracted from medical records. Percentiles were calculated for each, standardized by gestational age and infant sex, via Fenton growth curves [30]. Infants with standardized BWP below the 10th percentile are classified as small for gestational age (SGA); those above the 90th percentile as large for gestational age (LGA); and those between the 10th and 90th percentile as adequate for gestational age (AGA). SGA, AGA, and LGA are collectively referred to as birthweight groups (BWG). Mother’s BMI was calculated from height and weight collected from medical records. Using the mother’s pre-pregnancy BMI and gestational weight gain, we calculated her gestational weight gain category (low, normal, high), according to the Institute of Medicine/National Research Council 2009 gestational weight gain recommendations (Table 1) [31]. Briefly, recommended weight gain at term for underweight (BMI < 18.5), normal weight (BMI 18.5–24.9), overweight (25–29.9), and obese (BMI ≥30) mothers, respectively, are 28–40 lbs, 25–35 lbs, 15–25 lbs, and 11–20 lbs. Gestational weight gain below the recommended interval was coded as ‘low’, within the recommendation was coded as ‘normal’ and exceeding the recommendations was coded as ‘high’.

Additional study covariates

Foetal sex was obtained from medical records. Foetal race was self-reported and coded as ‘white’ or ‘not white’.

Tissue collection

The same placental sampling protocol was used by both cohorts. Foetal placental samples were collected within 2 hours of birth; sections were obtained two centimetres (cm) from the umbilical cord and free of maternal decidua. Collected tissue was immediately placed in RNA later solution (Life Technologies, Grand Island, NY, USA) and stored at 4°C for at least 72 hours. Subsequently, tissue segments were blotted dry and stored at −80°C.

NHBCS microRNA isolation and sequencing

Total RNA was extracted from fixed placenta using the Qiagen miRNeasy Mini Kit and a TissueLyser LT (Qiagen, Frederick, MD, USA) following the manufacturer’s protocol. Briefly, 25–35 mg of placental tissue was placed in a 2 ml round bottom tube with 700 µl of Qiazol Lysing Reagent and one 5 mm stainless steel bead. The tissue was homogenized on the TissueLyser LT for 2 minutes at 40 Hz. The resulting homogenate was processed with the Qiagen miRNeasy Mini Kit and eluted in 30 µl RNase-free water. The RNA was quantitated on a NanoDrop 2000 (Thermo Fisher, Waltham, MA, USA) and quality checked on Agilent Bioanalyzer using the Agilent RNA 6000 Nano kit (Agilent, Santa Clara, CA, USA). Single end, 1 × 75 bp next-generation sequencing of placental microRNA was performed by Qiagen Genomic Services (Frederick, Maryland).

RICHS microRNA isolation and sequencing

Total RNA was extracted from fixed, pulverized placenta using the Qiagen miRNeasy Mini Kit and a TissueLyser LT (Qiagen, Germantown, MD, USA) following the manufacturer’s protocol. Briefly, 25–35 mg of frozen, powdered placental tissue was placed in a 2 ml round bottom tube with 700 µl of Qiazol Lysing Reagent and one 5 mm stainless steel bead. The tissue was homogenized in a pre-chilled tube holder on the TissueLyser LT for two, 5-min cycles at 30 Hz. The resulting homogenate was processed with the Qiagen miRNeasy Mini Kit with on-column DNAse digestion and eluted in 50 µl RNase-free water. The RNA was quantitated on a NanoDrop (Thermo Fisher, Waltham, MA, USA) and quality checked on Agilent Bioanalyzer using the Agilent RNA 6000 Nano kit (Agilent, Santa Clara, CA, USA). Single end, 1 × 50 bp next-generation sequencing of placental microRNA was performed by Omega Bioservices (Norcross, Georgia).

RICHS total RNA isolation, sequencing, processing and quality control

Total RNA was extracted from fixed, pulverized placenta using the RNeasy Mini Kit (Qiagen) and stored at −80°C until analysis. We quantified RNA using a Nanodrop Spectrophotometer (Thermo Scientific), assessed RNA integrity via Agilent Bioanalyzer (Agilent), removed ribosomal RNA with a Ribo-Zero Kit [32], converted to cDNA using random hexamers (Thermo Scientific), and performed transcriptome-wide 50 bp single-end RNA sequencing via the HiSeq 2500 platform (Illumina) [33]. Initial quality control was performed on raw reads using the FastQC software. Two samples failed QC and were removed. Reads that passed the quality control metrics were mapped to the human reference genome (hg19) using the Spliced Transcripts Alignment to a Reference (STAR) aligner [34], with common SNPs in the reference genome masked prior to alignment. Raw reads have been deposited at the NCBI sequence read archive (SRP095910).

smallRNA-seq processing and quality control

Raw FASTQ reads obtained from a total of 552 (RICHS 230, NHBCS 322) samples were subject to adaptor trimming with cutadapt v1.16 [35]. For RICHS reads, the 3ʹ adaptor sequence were trimmed (TGGAATTCTCGGGTGCCAAGG) and then four bases were trimmed from each end of the read following vendor’s recommendation (BIOO scientific, Austin TX). For NHBCS reads, the 3ʹ adaptor sequence (AACTGTAGGCACCATCAAT) was trimmed based on vendor’s recommendation (Qiagen). After adaptor trimming, fastQC v0.11.5 was used to process the trimmed reads and QC results were aggregated using MultiQC v1.5 for visualization [36]. One sample from NHBCS failed QC and was removed. Then, we used trimmed reads and miRDeep2 to quantify microRNA [37]. In short, miRDeep2 was used to first perform alignment using bowtie1 with human genome hg38 [38]. The ‘Quantifier’ module in miRDeep2 was used to obtain raw counts of microRNAs with miRBase version 22 [39].

Sample filtering and transcript filtering

Raw counts were imported into DESeq2 for normalization and differential expression analysis. microRNAs with less than one count per million in more than 10% of samples were removed. Five samples were excluded from the RICHS microRNA and RNA data sets because they were missing race information, bringing the number of RICHS samples included in the analysis to 225 and 192 for the microRNA and RNA, respectively. Four samples in the NHBCS microRNA data set were missing birthweight and were excluded, bringing the number of NHBCS samples included in the microRNA analysis to 317. Of the 2,656 microRNA transcripts that mapped in both RICHS and NHBCS, 802 and 777 remained after filtering, respectively. Of the 50,810 total RNA transcripts that mapped in RICHS, 17,607 remained.

Normalization

Filtered microRNA and RNA raw counts were imported to DESeq2 for normalization and differential expression analysis. For all data sets, parametric estimates of dispersion were calculated, and the median ratio method was used to estimate size factors for normalization for modelling with DESeq2 [40]. Normalized counts were exported from DESeq2 for surrogate variable analysis. Count matrices normalized by the Variance Stabilizing Transformation (equivalent to normalization and log₂-transformation) were exported from DESeq2 for the microRNA and RNA datasets in RICHS, for correlation analysis.

Statistical analyses

Surrogate variable analysis

To adjust for unmeasured confounding effects, such as cell-type heterogeneity and other unmeasured sources of technical variation, we estimated surrogate variables from the normalized transcript reads (for microRNA and total RNA counts) via the svaseq function in the sva package [41,42]. In the svaseq function, the iteratively re-weighted least squares algorithm was used to estimate surrogate variables based on empirically derived control transcripts. The full model (mod argument) used for svaseq differed by application. The null model (mod0 argument) used all covariates except BWP. For RICHS and NHBCS microRNA differential expression analysis, one surrogate variable was estimated. For RICHS mRNA differential expression analysis, three surrogate variables were estimated.

Differential expression analyses

Transcript counts from the microRNA and total RNA sequencing datasets were modelled using a negative binomial generalized linear model with significance testing for differentially expressed transcripts via Wald tests in DESeq2 [43]. The microRNA transcripts were regressed on infant BWP (or BWG for follow-up analyses), with covariates for race (RICHS analyses only), sex, primigravida; and technical variables. In RICHS microRNA analysis, RIN (RNA integrity number) and sequencing lane number were included in the model as technical variables. In NHBCS microRNA analysis, RIN and sequencing pool were included in the model as technical variables. Surrogate variables were also included as covariates in our regression models (discussed in previous section). We considered microRNAs with a false discovery rate less than 5% to be primary differentially expressed (DEmiRs) with BWP.

In RICHS RNA analysis, RNA transcript abundances were regressed on BWG, with covariates for race, sex, primigravida, sequencing batch and three surrogate variables. Contrasts were used to estimate the log₂ fold change in the transcript for SGA and LGA, both relative to AGA. To assess the possible interaction of sex and birthweight, RNA transcripts were regressed on the interaction term of infant sex and BWG, with the same covariates listed previously. For these follow-up analyses, candidate transcripts with p-values less than 0.05 for either LGA or SGA were considered significant.

Differential expression comparison

Results from the RICHS and NHBCS microRNA differential expression analyses were filtered to the 719 microRNAs common to both analyses. The Pearson correlation of the RICHS and NHBCS differential expression analyses Wald statistics was calculated using the cor.test function in R.

Identifying mRNA targets of hsa-miR-532-5p

Potential mRNA targets hsa-miR-532-5p (miR-532) were identified from mirDIP, an online database of human microRNA-target predictions [44]. mirDIP integrates microRNA-target prediction across 30 different resources, providing nearly 152 million human microRNA-target predictions. Using the individual source ranking and confidence measures, mirDIP assigns a unified rank and confidence score using the quadratic function. For this analysis, the predictions with the top 5% confidence scores were returned [44]. The resulting miR-532-target mRNA pairs were carried forward in the analysis.

Putative target filtering

microRNA and total RNA counts were correlated using the Pearson method with the R library psych [45]. miR-532-target mRNA pairs that were negatively correlated (FDR < 0.05) were carried forward in the analysis and the associated mRNA were considered putative miR-532 target mRNAs. Finally, putative miR532 target mRNAs that were significantly associated with BWG (P < 0.05) in RICHS, were selected for pathway analysis.

Pathway analysis

Pathway over-representation analysis was conducted in consensuspathDB [46], which aggregates data from 12 separate pathway analysis databases. For each pathway gene set, consensus path DB calculates a p-value according to the hypergeometric test for the genes in both the miR-532 putative target genes and the pathway gene set. All mRNAs not filtered for low reads were included as the background, or null distribution, for the test.

Results

This study analysed data from 317 mother–infant pairs from the New Hampshire Birth Cohort (NHBCS) and 225 mother–infant pairs from the Rhode Island Child Health Study (RICHS). The demographics of the participants for this study shown in Table 1. A few important differences exist between the two cohorts. By design, RICHS contains a greater proportion of SGA and LGA newborns than NHBCS (Table 1). The majority of placentae in both studies were from full-term pregnancies and the median gestational ages were similar, NHBCS has a larger range of gestational ages than RICHS. Gestational age was not included in the final models because birthweight percentiles are calculated with respect to gestational age. Lastly, the RICHS cohort is more racially diverse than the NHBCS cohort (Table 1).

To analyse the associations between BWP and placental microRNA expression, we performed differential expression analysis. Within each cohort and for each microRNA transcript individually, transcript abundances were regressed on corresponding newborn BWP in DESeq2 [43], while also adjusting for primigravida, infant sex, race (white/not), technical covariates (RNA integrity number, sequencing lane, library pool) and one surrogate variable [42]. In NHBCS, 777 microRNA transcripts passed filtering, 7 had an FDR <0.05 and 2 were significant with a Bonferroni threshold (p-value <7e-05, Figure 1(a), Table S1). Of the 802 microRNA transcripts that passed filtering in RICHS, 11 had an FDR < 0.05 and 3 passed a strict Bonferroni threshold (p-value <7.2e-05, Figure 1(b), Table S2).

Figure 1. microRNA association with birthweight percentile. Volcano plots illustrate the results of the differential expression analyses for the New Hampshire Birth Cohort Study (round points; (a) and the Rhode Island child health study (square points; b), respectively. On the y and x axes, -log₁₀(p-values) in the association of each microRNA with BWP and effect estimates, or the log₂ fold change in each microRNA per one percent change in BWP, are shown respectively. The horizontal-dashed lines represent the Bonferroni threshold (log₁₀[0.05/M], where M is the number of microRNAs tested). (a) 7 microRNAs are significantly (FDR < 0.05) associated with BWP in NHBCS (grey filled shapes), with two microRNAs significant after Bonferroni correction. (b) 11 microRNAs are significantly (FDR < 0.05) associated with BWP in RICHS (grey filled shapes), with three microRNAs significant after Bonferroni correction. (c) Wald statistics for the association of each microRNA with BWP from the cohort-level differential expression analyses were plotted and labelled to indicate if they were: consistent in direction of effect (empty triangle), significant in RICHS (square), or significant in NHBCS (circle)

Of the microRNA transcripts that passed filtering in RICHS and NHBCS, 719 were common. Among the primary DEMirs, miR-532 had a significant positive association with BWP in both NHBCS (log₂ fold change = 2.4e-3; log₂ standard error: 6.7e-4; FDR = 0.04; Figure 1(a)) and RICHS (log₂ fold change = 2.4e-3; log₂ standard error: 5.7e-4; FDR = 6.6e-3; Figure 1(b)). Although there was little overlap between the microRNAs associated with BWP in the two cohorts, the Wald statistics from the cohort-level differential expression analysis were correlated (Pearson’s coefficient = 0.35, P = 1.3e-21), indicating that the NHBCS and RICHS differential expression analyses produced similar results despite sampling differences (Figure 1(c)).

In order to examine whether the placentas of babies born very large or very small exhibited different expression levels of miR-532, we also performed a differential expression analysis where microRNA transcript abundances were regressed on BWG (SGA/AGA/LGA). The positive correlation of miR-532 and birthweight was still evident in the extremes of birthweight groupings (Figure 2). When compared to the AGA group, there was a 0.12 and 0.24 log₂ fold change decrease in miR-532 in the SGA group, for RICHS and NHBCS, respectively. There is an apparent increase of miR-532 from the AGA to LGA groups (0.07 and 0.1 log₂ fold change, for RICHS and NHBCS), although the lower bounds of the confidence intervals crossed the null (Figure 2).

Figure 2. miR-532 associates with birthweight group. Estimates of log₂ fold change for SGA and LGA, both relative to AGA, illustrate a significant reduction of miR-532 in placentae from SGA newborns and an apparent increase in miR-532 in LGA newborns for NHBCS and RICHS. Birthweight groups, small (SGA) and large (LGA) for gestational age are represented by dashed and solid lines, respectively

miR-532 is associated with the expression of impactful targets

Bioinformatic targets of miR-532 were collected from the microRNA Data Integration Portal (mirDIP) [44]. Bioinformatic predictions are known to generate false positives. mirDIP predicted 1,883 targets of miR-532. Recent research has indicated that microRNA-guided transcriptional inhibition will result in transcript degradation, leading to less detectable target transcript [1]. In order to enrich our targets with true microRNA-target pairs, we utilized the parallel measures of total RNA and microRNA in the RICHS cohort. mRNA transcripts were considered putative targets of miR-532 if they were bioinformatically predicted targets of miR-532 and if the mRNA’s abundance was negatively correlated with miR-532 transcript abundance in RICHS (Table S3). Using these criteria, we found 56 putative miR-532 targets. We further curated our putative targets of miR-532 by requiring that they be nominally associated (p-value <0.05) with birthweight in the RICHS cohort. These 18 putative miR-532 target mRNAs used in pathway analysis (Figure 3, Table S3).

Figure 3. miR-532 target transcript abundance association with LGA or SGA, relative to AGA. Log₂ fold change of miR-532 putative target transcript abundance in placentae from SGA and LGA newborns, relative to AGA newborns

The 18 putative miR-532 target genes were tested for pathway overrepresentation with consensuspathdb, against all of the genes that passed QC in the RICHS whole transcriptome RNA-seq analysis [46]. Pathways related to adipogenesis, adipocytokine signalling and energy metabolism were over-represented among miR-532 putative targets. These enrichment results are primarily driven by putative miR-532 target mRNA transcripts: Leptin (LEP), Adiponectin Receptor 1 (ADIPOR1), Adiponectin Receptor 2 (ADIPOR2) and Solute Carrier Family 2 Member 1 (SLC2A1, GLUT1) (Figure 4).

Figure 4. Pathway analysis of miR-532 putative targets. Select consensus path DB results are presented in a bubble plot. Dot colour represents -log₁₀(FDR). Dot size indicates the percentage of pathway members that are represented by miR-532 targets. Pathway sources are indicated on the right y-axis

Sex bias in the associations of LEP with birthweight

We have previously shown a sexually dimorphic association of DNA methylation at the LEP promoter with birthweight in RICHS placentae [47]. Five of the 18 putative miR-532 target mRNAs, including LEP had significant interactions between sex and BWG (Table S4). Because of the enrichment of adipokine pathways in miR-532 target mRNAs and our previous findings regarding LEP, we used a stratified analysis to assess sex-specific differences. We found that LEP transcript abundance is negatively associated with BWP in males, but not associated with birthweight in females (Figure 5(b)). There is a one-fold decrease in LEP from AGA to LGA and a reciprocal (though not significant) one-fold increase in LEP from AGA to SGA group. For females in this analysis, there are apparent (though not significant) increases in LEP in SGA and LGA groups, relative to the AGA (Figure 5(a)).

Figure 5. Sex-specific association of birthweight percentile and group with placental leptin transcript abundance in RICHS. (a) log₂ fold change of LEP abundance in placentae from SGA and LGA newborns, relative to AGA newborns. (b) Estimates of log₂ fold change per one BWP. Colour represents sample stratification: black (all placentae), red (female placentae) and blue (male placentae)

Discussion

In this study, we have used an integrated, multi-omics approach to describe the microRNAs that associate with birthweight in human term placenta. We performed these analyses across two mother-infant cohorts, NHBCS and RICHS, to compare consistency in the overall miRNA relationships, which revealed an overall correlation in the relationships and one miRNA that achieved an FDR < 0.05 in both. We found evidence of differentially expressed microRNAs in both cohorts, with one microRNA found to be similarly and statistically significantly differentially expressed in both cohorts. We utilized parallel whole and small RNA-seq in RICHS to narrow microRNA-target mRNAs to only those that were negatively correlated with miR-532 transcript abundance. We further filtered miR-532 target mRNAs to only those displaying a nominally significant association with birthweight in the RICHS cohort. These miR-532 mRNA targets were used in a pathway analysis, revealing enrichments of adipogenesis, adipocytokine signalling and energy metabolism pathways.

The roles of microRNAs are dependent on the mRNAs they regulate, and these interactions are tissue and context-specific. Few non-targeted microRNA studies in placenta exist; fewer studying microRNA and birthweight. To our knowledge, miR-532 expression in placenta has not been previously associated with birthweight. In two studies of birthweight, miR-532 was chosen as an endogenous control, due to stable transcript levels [11,23]. However, in the study from Timofeeva et al., miR-532 was differentially expressed in the placenta of women with preeclampsia, suggesting it may not be a suitable stable control [11]. In this study, miR-532 was differentially expressed with birth weight percentiles and when comparing SGA and LGA to AGA in two independent cohorts, with similar estimates of effect in both cohorts, indicating that this association is robust and reproducible.

While we focus on miR-532, which had the strongest evidence for an association with birth weight, several other microRNAs yielded strong associations with BWP in either NHBCS or RICHS, some of which were consistent with the findings of other studies (Table S5). Among microRNA with an FDR < 0.1 in NHBCS, three were also found in similar studies (hsa-miR-193b-3p [21,28], hsa-miR-181a [23], hsa-miR-101-5p [21]). In RICHS, nine microRNAs were consistent with previous findings (hsa-miR-520a-3p [21], hsa-miR-193b-5p [21], hsa-miR-210-3p [21], hsa-miR-193b-3p [21,28], hsa-miR-335-3p [28], hsa-miR-515-5p [21,22], hsa-miR-519d-5p [21,22,26], hsa-miR-365b-3p [21], hsa-miR-365a-3p [21]). It is reassuring that despite different sampling schemes, study designs and nucleic acid technologies – we see an overlap in these findings even though they did not meet the statistical threshold of FDR < 0.05 in both cohorts for our study.

Bioinformatic gene target prediction and integration with mRNA target transcript abundances indicate that gene targets of miR-532 in our samples were enriched in adipogenesis and energy metabolism pathways. The adipokines leptin and adiponectin modulate insulin sensitivity, metabolism, and energy homoeostasis in adulthood [48,49]. While their roles in placenta are less clear, both of these hormones seem to be involved in placental development and function [50,51]. LEP has a placenta-specific enhancer and in placenta, leptin has autocrine functions related to angiogenesis, hypoxia response, and endocrine functions including mobilization of maternal fat [52]. Unlike leptin, adiponectin is not produced by the placenta [53]. However through the placental expression of adiponectin receptors, maternal adiponectin might limit foetal growth by inhibiting the production of placental lactogen, chorionic gonadotropin, and progesterone in trophoblast cells [53,54]. Both of these hormones have been linked to gestational diseases like preeclampsia [55–59], gestational diabetes [59–62] and foetal outcomes like intrauterine growth restriction [47,59,60,63–67]. While earlier work did find positive associations between placental leptin expression and birthweight [60,61,68], several studies since have found a negative correlation – which is in line with our findings [64,66,69–73].

In adulthood, leptin is known to be responsive to gonadal steroids, resulting in a sexually dimorphic response to leptin signalling [74]. We have previously reported sexually dimorphic methylation patterns in the LEP promoter [47]. Lesseur et al. found that RICHS male placentae had higher levels of promoter methylation than female placentae [47]. Assuming the canonical negative association of promoter DNA methylation and gene expression, those results are consistent with our own, in which the mean LEP transcript counts for male placentae were lower than for female placentae. To our knowledge, we are the first research study to find a dimorphic association of LEP expression with birthweight in placenta. Cagnacci et al. report that leptin in amniotic fluid at 16 weeks is inversely correlated with foetal growth in females and not correlated with any growth parameter in males [73]. This study makes a case for sexual dimorphic negative association between leptin and foetal growth, but in the opposite sex as in our study. We found that LEP expression inversely correlated to BWP in newborn males, and was not associated with BWP in females. The findings described in the Cagnacci study are at a much earlier time in gestation and in a different tissue. At 16 weeks of gestation, the placental production of leptin is low [75] and the amniotic fluid is supplied primarily by maternal serum [76], suggesting that leptin concentration could be a representation of the mother’s leptin status and unrelated to placental leptin [73].

To our knowledge, this is the largest study yet to examine the relationships between human placental microRNAs and birthweight and the only study in which results are replicated in an independent cohort. Our results suggest a role for miR-532 in the regulation of foetal growth through control of genes that control energy metabolism. However, our findings should be interpreted within the context of this study’s limitations.

This is an observational study in which RNA was assayed from term placentae of live births. We cannot conclude that our results are representative of microRNA associations throughout development, and we cannot determine the directionality of these observed relationships since BWP and microRNA were measured cross-sectionally. Although our results highlight a microRNA that may be involved in birthweight, we cannot claim that our findings represent a causal link between miR-532 and birthweight. There is a potential for selection bias in our study, under the assumption that both microRNA expression and intrauterine growth restriction can lead to pregnancy loss [77]. In this scenario, selection of only live births could lead to erroneous associations between microRNA expression and birthweight. While studies have found placental microRNAs that were differentially expressed in early pregnancy loss vs elective termination [78,79], we do not expect most microRNAs to have a strong enough effect, individually, to lead to foetal demise [80]. We are not aware of any evidence that placental miR-532 is essential to life. We adjusted for likely confounders in our study, but cannot rule out the possibility that unmeasured or residual confounding remains in our analysis. Two potential confounders are maternal BMI and gestational weight gain. These variables were not included in the analysis, because they were not available for the NHBCS cohort at the time of this analysis. Sensitivity analysis suggests that the addition of maternal BMI and gestational weight gain to our models would have little to no effect on our conclusions, though some of our top findings are sensitive to one or both (Tables S6-S11, Figures S1 and S2). To limit statistical confounding due to unknown confounding (e.g. differing cellular composition between individuals or population stratification), our models are adjusted using surrogate variable analysis. This is a data-driven approach that may incorrectly estimate confounding elements, adding variability to our model. Recent research indicates that SVA is one of the more robust and reliable methods in studies such as ours [81]. Lastly, both cohorts utilized in this study consisted predominantly of healthy white mothers from the New England region of the United States, potentially limiting the generalizability of our findings.

Conclusion

MicroRNAs are an important regulatory mechanism in gene expression. The activity of microRNAs has been linked to the development and function of placentae throughout gestation, making them essential to foetal development. Here we have shown that the transcript abundance of miR-532 is associated with birthweight in independent cohorts. Putative targets of miR-532 are involved in angiogenesis, energy metabolism and endocrine signalling. Although targets of miR-532 are already established as regulators of placental development and function, ours is the first study to report the association of miR-532 and birthweight.

Disclosure statement

The authors have declared no competing interest.

Supplementary Material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Institutes of Health under Grants NIH-NIEHS R24 ES028507, NIH-NIEHS R01 ES025145, NIH-NIEHS P30 ES019776, and NIH-NIMHD R01 MD011698. The study sponsors did not have any role in the study design, collection, analysis, and interpretation of the data, the writing of the report, or the decision to submit the paper for publication.