Regulation of birthweight by placenta-derived miRNAs: evidence from an arsenic-exposed birth cohort in Bangladesh

ABSTRACT

Altered expression of microRNAs (miRNAs) is implicated in fetal growth. However, the mechanisms by which placenta-derived miRNAs regulate birthweight are not well understood. In Phase 1, we compared the expression of 754 miRNAs in the placenta from two extreme birthweight groups (0.8–2.2 kg vs. 3.3–3.9 kg, n = 77 each) selected from an arsenic-exposed Bangladeshi birth cohort (n = 1,141). We identified 49 miRNAs associated with the extreme birthweight groups and/or gestational age in Phase 1, which were further analyzed in Phase 2 among 364 randomly selected mother-infant pairs. Gestational age was determined by ultrasound. Causal mediation analysis was used to estimate the effect of miRNAs on birthweight considering gestational age a mediator, adjusting for cord blood arsenic and other risk factors. miR-1290, miR-195, and let-7g showed significant inverse associations with gestational age, while miR-328 showed significant positive association [false discovery rate (FDR) <0.05]. Via changing gestational age, miR-1290, miR-195, and miR-27a showed significant inverse associations with birthweight, while miR-328 and miR-324-5p showed significant positive associations (FDR <0.05). The effect of miRNAs on birthweight varied by gestational age (for miR-1290, miR-195, miR-328) and in utero arsenic exposure (for miR-1290): stronger effect was observed among infants delivered early in gestation or exposed to higher concentrations of arsenic in cord blood. Gene enrichment analysis with in silico predicted targets identified cell proliferation, inflammation, apoptosis, insulin, and IGF family signaling cascades associated with these miRNAs. Future studies are warranted to replicate these findings and assess these miRNAs as early biomarkers of fetal growth.

KEYWORDS:

• microRNA
• placenta
• birth weight
• gestational age
• epigenetics
• arsenic
• miR-1290
• miR-195
• miR-328

Introduction

Low birthweight is an important population indicator of neonatal mortality and morbidity and an established risk factor for later-life chronic diseases, particularly cardiovascular disease and metabolic syndrome [1, 2]. Birthweight is governed by two main causal processes: shortening of gestation and intra-uterine growth restriction [3]. The placenta is thought to play a critical role in both processes, as it is involved in facilitating the supply of oxygen and nutrients to the fetus, secreting hormones and growth factors, and detoxifying potentially harmful maternal and environmental factors [4,5]. However, the molecular mechanisms by which the placenta regulates birthweight are not well understood, although epigenetic alterations – changes that regulate gene activity without the change in underlying DNA sequence – are potential players to influence fetal growth and the programming of health throughout the life-course [6].

One such epigenetic mechanism is the aberrant expression of microRNAs (miRNAs). microRNAs are 18–25 nucleotides long small noncoding RNA molecules that are involved in the regulation of gene expression at the posttranscriptional level [7]. Binding to the 3'-untranslated region of the messenger RNA (mRNA) target, miRNAs can effectively silence gene expression either by translational repression or mRNA degradation [8]. It is estimated that up to 60% of all human protein-coding genes are regulated by miRNAs [9], essentially controlling a variety of cellular processes, including metabolism, cell proliferation, differentiation, apoptosis, and development [10,11]. miRNAs are ubiquitously expressed in human tissues, although some are largely unique to the placenta [12]. Placenta-derived miRNAs have been implicated in normal placental development, function, and remodeling; trophoblast cell proliferation, migration, and invasion; apoptosis; and angiogenesis [13–15].

Several studies have investigated the link between miRNA expression profiles in maternal and fetal tissues and adverse maternal and perinatal outcomes. Most of these studies have focused on preeclampsia [12,14,16–20], although a few studies examined the associations with preterm birth (<37 weeks of gestation) [12,14,16–20], fetal growth restriction [18,21,22], and low birthweight (<2500 g) [23], reporting dysregulated expression of some miRNA species in relation to these outcomes. However, results from these studies are often conflicting or irreproducible, mainly because of small sample size, differences in miRNA source tissue and expression technology, inadequate adjustment for confounders, and targeted miRNA expression analysis approach adopted by most studies. Nonetheless, combined evidence suggests a role for miRNAs in determining birthweight, possibly involving both preterm birth and fetal growth restriction within the causal pathways. However, epidemiological studies investigating the causal relation between placenta-derived miRNAs and birthweight are lacking.

Furthermore, the expression of miRNAs is sensitive to environmental exposures, including inorganic arsenic [24], which is an independent risk factor for preterm birth [25] and lower birthweight [26,27]. Therefore, it is possible that the effect of placenta-derived miRNAs on birthweight can be modified by in utero arsenic exposure. In this study, we aimed to investigate the causal relation between placenta-derived miRNAs and birthweight, adjusting for in utero arsenic exposure and other risk factors, using data from mother-infant pairs enrolled in a prospective birth cohort in Bangladesh. We hypothesized that placenta-derived miRNAs will regulate birthweight involving pathways mediated through gestational age and pathways independent of gestational age (Figure 1). We also aimed to investigate whether the effect of placenta-derived miRNAs on birthweight was modified by in utero arsenic exposure.

Results

Phase 1 results

The characteristics of the study participants are summarized in Table 1. The median (range) birthweight of the two extreme birthweight groups selected in Phase 1 was 2.1 (0.8–2.2) kg and 3.5 (3.3–3.9) kg, respectively. The two extreme birthweight groups significantly (P value <0.05) differed by gestational age, maternal enrollment BMI, education, study site, and cord blood arsenic exposure.

Table 1. Distributions of selected characteristics of the study participants by study phases. Data were expressed as Mean ± SD or Median (range) except where indicated.

	Phase 1
Characteristics	Upper extreme birthweight group (n = 77)	Lower extreme birthweight group (n = 77)	P value	Phase 2 (n = 364)
Birthweight (kg)	3.5 (3.3 – 3.9)	2.1 (0.8 – 2.2)	<0.001	2.9 (1.0−4.3)
Maternal age (years)	23.5 ± 4.3	23.3 ± 4.2	0.7	22.7 ± 3.7
Gestational age (weeks)	38.5 ± 1.4	35.1 ± 3.0	<0.001	37.9 ± 2.0
Enrollment BMI, n (%)			0.03
≤ 18.5 kg/m^2	12 (15.6)	23 (29.9)		92 (25.3)
> 18.5 to ≤ 25.0 kg/m^2	53 (68.8)	50 (64.9)		239 (65.7)
> 25.0 kg/m^2	12 (15.6)	4 (5.2)		33 (9.0)
Maternal Education, n (%)			0.01
No formal education	5 (6.5)	18 (23.4)		57 (25.7)
Primary education	21 (27.3)	20 (26.0)		127 (34.9)
Secondary or higher	51 (66.2)	39 (50.6)		180 (49.4)
Infant sex, n (%)			0.4
Female	43 (55.8)	38 (49.4)		179 (49.2)
Male	34 (44.2)	39 (50.6)		185 (50.8)
Study site, n (%)			<0.001
Sirajdikhan	33 (42.9)	10 (13.0)		173 (47.5)
Pabna	44 (57.1)	67 (87.0)		191 (52.5)
No. of past pregnancies, n (%)			0.5
0	27 (35.1)	26 (33.8)		143 (39.5)
1	28 (36.4)	23 (29.9)		118 (32.6)
> 1	22 (28.6)	28 (36.4)		101 (27.9)
Cord blood arsenic (μg/dl)	0.54 (0.1–7.1)	0.70 (0.09–6.0)	0.003	0.60 (0.1–15.5)

Associations between miRNA expression and gestational age

In Phase 1, we selected 25 placenta-derived miRNAs for further analysis based on their associations with gestational age at the level of permutation P <0.01 (Table S1). Twenty-four of these miRNAs were associated with gestational age in at least 2 out of 4 normalization methods, while 16 miRNAs (let-7g, miR-16, miR-27a, miR-34b, miR-138, miR-181a, miR-181c, miR-195, miR-199a, miR-205, miR-296, miR-328, miR-331, miR-489, miR-542-3p, and miR-886-3p) were associated with gestational age in at least 3 normalization methods.

Associations between miRNA expression and the extreme birthweight groups

In Phase 1, we also selected 33 miRNAs for further analysis based on their associations with the extreme birthweight groups (Table S2). Twenty-three of these miRNAs were associated with the extreme birthweight groups in ≥2 normalization methods at the level of permutation P <0.01. Thirteen of these miRNAs (miR-124a, miR-625, miR-524, miR-518f, miR-376b, miR-34c, miR-337-3p, miR-324-5p, miR-30c, miR-205, miR-195, and miR-185) showed significant associations in ≥3 normalization methods. Among these 33 miRNAs, 9 miRNAs (let-7g, miR-138, miR-195, miR-205, miR-324-5p, miR-337-3p, miR-34b, miR-522, and miR-625) were also associated with gestational age in ≥2 normalization methods (Figure S1). We also selected 10 miRNAs that showed significant associations with the extreme birthweight groups in any normalization method but have known biological significance.

Phase 2 results

The median (range) birthweight of the participants selected in Phase 2 was 2.9 kg (1.0–4.3). Other characteristics of Phase 2 participants are presented in Table 1.

Associations between miRNA expression and gestational age

In Phase 2, we identified 9 miRNAs (let-7g, miR-1290, miR-195, miR-200c, miR-27a, miR-30c, miR-324-5p, miR-328, and miR-424) significantly associated with gestational age at the level of P <0.05 (Table 2). While 4 of these miRNAs (let-7g, miR-1290, miR-195, and miR-27a) were negatively associated with gestational age, 5 miRNAs (miR-200c, miR-30c, miR-324-5p, miR-328, and miR-424) showed significant positive associations at the level of P <0.05. After adjusting for multiple testing, the associations for let-7g, miR-1290, miR-195, and miR-328 remained statistically significant (FDR (false discovery rate) <0.05). Our results showed that a 1-standard deviation (SD) increment in placental expression of let-7g, miR-1290, miR195, and miR-328 in the placenta was associated with a change in gestation age of −0.46 weeks (95% CI: −0.78, −0.14), − 0.67 weeks (95% CI: −0.96, −0.38), −0.56 weeks (95% CI: −0.82, −0.30), and + 0.45 weeks (95% CI: + 0.26, + 0.64), respectively.

Table 2. The effects of a 1-SD increment in placental expression of miRNA on gestational age and birthweight in Phase 2 study (n = 364). The indirect, direct, and total effects of miRNA expressions on birthweight were estimated using causal mediation analysis considering gestational age a mediator.

				Associations with Birthweight
	Associations with Gestational Age			Indirect Effects			Direct Effects			Total Effects
miRNA	Coefficient (95% CI)	P value	FDR	Coefficient (95% CI)	P value	FDR	Coefficient (95% CI)	P value	FDR	Coefficient (95% CI)	P value	FDR	PM ^b %
let7g	−0.46 (−0.78, −0.14)	0.005	0.042	−37.5 (−77.4, −1.8)	0.044	0.132	18 (−86.5, 115.3)	0.798	0.582	−19.5 (−145.6, 79.8)	0.716	0.676
miR-124a	0.1 (−0.49, 0.69)	0.736	0.581	1.6 (−66, 74.4)	0.910	0.454	−46.3 (−212.2, 123.3)	0.650	0.547	−44.7 (−224.2, 137.2)	0.674	0.676
miR-127	−0.03 (−0.25, 0.18)	0.774	0.581	−1.7 (−23.5, 15.6)	0.832	0.453	27.3 (−25.7, 75.4)	0.300	0.453	25.7 (−32, 78.9)	0.360	0.539
miR-1290^a	−0.67 (−0.96, −0.38)	<0.001	<0.001	−52.9 (−95.4, −22.6)	0.000	0.003	−107.1 (−174.3, −31.8)	0.005	0.080	−160.1 (−235.6, −74)	<0.001	0.009	33.1
miR-130a	0.02 (−0.33, 0.37)	0.910	0.608	7.7 (−17.7, 41.4)	0.530	0.413	−59.8 (−150.2, 27)	0.194	0.443	−52.1 (−148.4, 37.1)	0.278	0.531
miR-135b	0.1 (−0.24, 0.43)	0.569	0.554	−6.4 (−44.5, 26.7)	0.764	0.444	−12.7 (−105, 71)	0.748	0.570	−19.1 (−111, 74.7)	0.686	0.676	33.6
miR-138	−0.2 (−0.51, 0.11)	0.210	0.385	−13.5 (−39.5, 2.6)	0.152	0.302	−50.1 (−102.2, −5.8)	0.042	0.363	−63.6 (−124.9, −13.9)	0.018	0.216	21.2
miR-149	−0.04 (−0.44, 0.35)	0.832	0.601	−9 (−42.6, 21.7)	0.586	0.413	−44.6 (−126.3, 46.2)	0.324	0.453	−53.6 (−136.6, 50.8)	0.270	0.531	16.8
miR-15a	−0.25 (−0.56, 0.07)	0.128	0.328	−15.1 (−43.4, 13.3)	0.296	0.327	3.2 (−79.3, 75.2)	0.984	0.642	−11.9 (−95.2, 65.3)	0.736	0.676
miR-16	−0.14 (−0.42, 0.15)	0.342	0.440	−11.1 (−34.5, 8.8)	0.300	0.327	42.2 (−13.6, 101.6)	0.140	0.413	31 (−24.9, 95.6)	0.306	0.531
miR-181a	−0.02 (−0.28, 0.23)	0.848	0.601	−3.8 (−23.7, 9.5)	0.562	0.413	−41.3 (−88, −0.4)	0.054	0.363	−45.1 (−98.3, −4.6)	0.040	0.240	8.3
miR-181c	−0.16 (−0.52, 0.19)	0.373	0.440	−7.5 (−33.2, 10.3)	0.450	0.387	−56.3 (−121, −0.6)	0.058	0.363	−63.8 (−135.8, −3.5)	0.038	0.240	11.7
miR-185	−0.42 (−0.93, 0.09)	0.106	0.300	−20.2 (−73.1, 26.1)	0.468	0.387	−8.4 (−130.2, 123.3)	0.878	0.611	−28.6 (−162.5, 107.7)	0.708	0.676	70.7
miR-195	−0.56 (−0.82, −0.3)	<0.001	<0.001	−40.9 (−78.3, −13.8)	0.002	0.024	−47.1 (−120.9, 34.2)	0.278	0.453	−87.9 (−170.4, 0.8)	0.053	0.271	46.5
miR-199a	0.27 (−0.11, 0.64)	0.164	0.366	35.1 (8, 70)	0.014	0.056	−70 (−167.3, 10.9)	0.112	0.413	−34.9 (−139.3, 49.8)	0.520	0.623
miR-200c	0.42 (0.03, 0.82)	0.038	0.131	19.8 (−8.7, 51.8)	0.230	0.327	−26.9 (−129, 82.7)	0.716	0.570	−7 (−114.9, 114.7)	0.958	0.702
miR-205	0.15 (−0.33, 0.63)	0.530	0.538	17.9 (−13.4, 53.2)	0.274	0.327	−37.9 (−127.5, 56.8)	0.448	0.541	−20.1 (−115.2, 77.6)	0.720	0.676
miR-20b	−0.23 (−0.52, 0.07)	0.137	0.328	−11.4 (−48.6, 17.9)	0.464	0.387	−31.6 (−108, 44.9)	0.374	0.498	−43 (−129.1, 43.8)	0.332	0.531	26.5
miR-21	−0.06 (−0.4, 0.28)	0.736	0.581	−2.2 (−47.1, 25)	0.902	0.454	31.4 (−54.6, 122.5)	0.550	0.541	29.2 (−68.2, 119.4)	0.596	0.669
miR-222	0.16 (−0.16, 0.49)	0.322	0.440	15.6 (−5.9, 48.4)	0.226	0.327	−17.8 (−101, 61.3)	0.612	0.541	−2.2 (−83.7, 79.1)	0.902	0.702
miR-27a	−0.4 (−0.74, −0.06)	0.023	0.089	−36.1 (−69.5, −11.4)	0.010	0.048	27.6 (−77.5, 100.1)	0.582	0.541	−8.5 (−123.9, 64)	0.774	0.678
miR-296	0.26 (−0.13, 0.64)	0.189	0.385	13.8 (−7.9, 44.2)	0.272	0.327	−46.4 (−105.1, 17.1)	0.140	0.413	−32.6 (−91.9, 37.2)	0.320	0.531
miR-30c	0.32 (0.05, 0.58)	0.019	0.086	15.1 (0.7, 40.3)	0.070	0.168	15.6 (−40.4, 89.2)	0.586	0.541	30.7 (−22.4, 96.8)	0.276	0.531	49.3
miR-320b	−0.05 (−0.36, 0.26)	0.750	0.581	3.4 (−19.6, 33.6)	0.802	0.447	33.6 (−30.9, 107.5)	0.312	0.453	37 (−31.2, 117.2)	0.328	0.531	9.1
miR-3245p	0.32 (0.08, 0.57)	0.009	0.051	28 (9.2, 54.9)	0.004	0.032	−5.9 (−72.4, 48)	0.800	0.582	22 (−41.6, 74.8)	0.500	0.620
miR-328	0.45 (0.26, 0.64)	<0.001	<0.001	20.4 (5.8, 45.4)	0.007	0.040	11.1 (−31.7, 63.9)	0.551	0.541	31.4 (−7.7, 86.5)	0.123	0.489	64.8
miR-331	0.13 (−0.15, 0.42)	0.363	0.440	5.2 (−16.8, 32.2)	0.660	0.439	44.3 (−31.2, 107.7)	0.226	0.453	49.5 (−25.9, 116.7)	0.184	0.531	10.5
miR-3373p	−0.12 (−0.4, 0.15)	0.381	0.440	−4.9 (−37.5, 20)	0.726	0.444	1.2 (−57.9, 71.7)	0.966	0.642	−3.7 (−78.3, 68.8)	0.928	0.702
miR-34b	0.22 (−0.12, 0.56)	0.201	0.385	7 (−26.4, 39.3)	0.704	0.444	27.7 (−60.3, 109.5)	0.538	0.541	34.7 (−56.6, 118.1)	0.492	0.620	20.2
miR-34c	−0.07 (−0.31, 0.18)	0.592	0.560	−8.1 (−35.1, 10.9)	0.464	0.387	−29.4 (−74.8, 24.1)	0.326	0.453	−37.5 (−90.1, 19.9)	0.220	0.531	21.5
miR-373	0.05 (−0.29, 0.39)	0.782	0.581	1 (−27.8, 30.5)	0.908	0.454	6.6 (−78.4, 83.1)	0.942	0.641	7.6 (−83.2, 84.8)	0.920	0.702	13.5
miR-376b	−0.23 (−0.65, 0.2)	0.303	0.440	−8.6 (−47.2, 24.3)	0.652	0.439	18.3 (−64.6, 101)	0.620	0.541	9.7 (−85.4, 101.5)	0.752	0.676
miR-424	0.33 (0.08, 0.57)	0.010	0.051	25.1 (5.5, 53.6)	0.024	0.082	−17.3 (−77.1, 45.2)	0.554	0.541	7.8 (−59.5, 69.3)	0.828	0.692
miR-454	−0.02 (−0.37, 0.32)	0.894	0.608	15.4 (−23.2, 56.8)	0.442	0.387	−16.8 (−128.2, 80.5)	0.674	0.553	−1.4 (−127.6, 106.1)	0.896	0.702
miR-455	0.13 (−0.29, 0.55)	0.534	0.538	19.1 (−15, 58.3)	0.264	0.327	45.1 (−33.7, 134.5)	0.274	0.453	64.3 (−31.3, 159.4)	0.148	0.531	29.7
miR-489	−0.21 (−0.45, 0.03)	0.084	0.262	−13.5 (−32.5, 4.9)	0.164	0.302	−57.5 (−109.4, 2.8)	0.068	0.363	−71 (−129.2, −3)	0.034	0.240	19.0
miR-518e	−0.21 (−0.56, 0.14)	0.241	0.417	−19.9 (−50.9, 0.9)	0.092	0.200	44.4 (−48.5, 136.4)	0.296	0.453	24.6 (−73.3, 118.8)	0.538	0.624
miR-518f	0 (−0.31, 0.31)	0.989	0.630	−4.9 (−37.7, 26.8)	0.778	0.444	−44.3 (−139.6, 34.2)	0.324	0.453	−49.1 (−146.8, 39.4)	0.330	0.531	10.0
miR-522	0.07 (−0.26, 0.4)	0.682	0.581	−2.8 (−25.1, 18.5)	0.754	0.444	52.3 (−16.1, 130.9)	0.142	0.413	49.5 (−32.1, 125.6)	0.198	0.531
miR-524	0.01 (−0.33, 0.36)	0.936	0.608	−5.4 (−31.8, 23.8)	0.728	0.444	17.5 (−71, 111.3)	0.740	0.570	12.1 (−84.1, 111.5)	0.828	0.692
miR-5423p	0.13 (−0.15, 0.41)	0.372	0.440	13.4 (−13.1, 39.9)	0.346	0.360	21 (−60.6, 85.9)	0.626	0.541	34.4 (−59.3, 105)	0.472	0.620	39.0
miR-625	−0.13 (−0.44, 0.19)	0.438	0.472	−15.8 (−50.9, 10.9)	0.274	0.327	46.1 (−29.3, 119.4)	0.174	0.443	30.3 (−53.6, 112.5)	0.384	0.552
miR-629	0.01 (−0.32, 0.35)	0.932	0.608	−5.7 (−27.8, 18.5)	0.580	0.413	45.5 (−22.7, 114.5)	0.194	0.443	39.8 (−37.7, 113.8)	0.308	0.531
miR-720	0.13 (−0.1, 0.35)	0.274	0.440	9.2 (−6.9, 27.3)	0.284	0.327	75.1 (26, 116.1)	0.002	0.070	84.2 (31.4, 128.1)	0.001	0.022	10.9
miR-8863p	0.07 (−0.31, 0.46)	0.718	0.581	−1.9 (−35.8, 22.2)	0.860	0.454	24.2 (−57.9, 106.2)	0.622	0.541	22.3 (−71.1, 108)	0.682	0.676
miR-9	0.11 (−0.17, 0.39)	0.439	0.472	6.7 (−12.6, 34.8)	0.542	0.413	−5.2 (−65.4, 56.4)	0.868	0.611	1.5 (−53.1, 62.7)	0.944	0.702
miR-99a	−0.1 (−0.53, 0.33)	0.648	0.577	−27 (−61.2, 0.9)	0.064	0.168	−27 (−119.4, 80.3)	0.604	0.541	−54 (−157.6, 64)	0.340	0.531	50.0
miR-99b	0.17 (−0.15, 0.48)	0.299	0.440	9.4 (−8, 37.2)	0.412	0.387	61.8 (−15, 117.7)	0.098	0.413	71.2 (−6.1, 128.7)	0.066	0.296	13.2
miR-99b*	−0.16 (−0.81, 0.48)	0.618	0.567	1.8 (−40.1, 45.7)	0.940	0.460	46.2 (−68.3, 154)	0.412	0.527	48 (−73.2, 159.6)	0.422	0.583	3.7

The models for gestational age (mediator model) were adjusted for maternal age, education, no. of past pregnancies, BMI, infant sex, study site, natural log of cord blood arsenic exposure, and surrogate variables. The models for birthweight were additionally adjusted for gestational age and included a miRNA-gestational age interaction term.

^a The outcome model for miR-1290 also included an interaction term between miR-1290 and natural log cord blood arsenic exposure.

^b PM = Proportion mediated. Calculated for miRNAs that showed same direction of associations via both indirect and direct effect pathways in relation to birthweight.

Associations between miRNA expression and birthweight

Causal mediation analysis revealed that placenta-derived miRNAs were associated with birthweight via changing gestational age and via pathways independent of gestational age. For example, let-7g, miR-1290, miR-195, miR-199a, miR-27a, miR-324-3p, miR-328, and miR-424 were significantly associated with birthweight via changing gestational age (indirect effect) at the level of P <0.05 (Table 2). After adjusting for multiple testing, the associations for 5 miRNAs – miR-1290, miR-195, miR-27a, miR-324-3p, and miR-328 – remained statistically significant (FDR <0.05). Three miRNAs – miR-1290, miR-138, and miR-720 – showed significant associations with birthweight (P <0.05) via pathways independent of gestational age (direct effect); none of the associations remained statistically significant after adjusting for multiple testing. Six miRNAs – miR-1290, miR-138, miR-181a, miR-181c, miR-489, and miR-720 – showed significant total effect associations with birthweight at the level of P <0.05. After adjusting for multiple testing, the total effect of associations for miR-1290 and miR-720 remained statistically significant (FDR <0.05). The indirect and direct effect associations for let-7g, miR-199a, miR-27a, miR-324-5p, and miR-424 were in the opposite direction, attenuating the total effect associations for these miRNAs on birthweight. miR-1290 was the only miRNA that showed statistically significant associations with birthweight via indirect, direct, and total effect pathways at the level of P <0.05. Specifically, a 1-SD increment in placental expression of miR-1290 was associated with an average of 160.1g (95% CI: −235.6, −74.0; FDR = 0.01) lower birthweight (total effect); 33.1% of the effect was mediated by a shortening of gestational age (β = −52.9g, 95% CI: −95.4, −22.6; FDR = 0.003), while the effect via pathways independent of gestational age was −107.1g (95% CI: −174.3, −31.8; FDR = 0.08). Sensitivity analysis indicated that the estimated mediation effects for our top-hit miRNAs were robust to the assumption of no unmeasured confounding effect (results not shown).

Gestational age also moderated the effect of placenta-derived miRNAs on birthweight. Specifically, the effects of miR-1290, miR-195, miR-328, and miR-489 on birthweight were stronger in magnitude and statistical significance among infants delivered early in gestation (all P-for-interaction <0.05) (Figure 2; Table S3). For example, higher expression of miR-1290 and miR-195 and lower expression of miR-328 in the placenta was associated with more than additively lower birthweight among infants born, roughly, before 38 weeks of gestation.

The effect of miR-1290 on birthweight was additionally moderated by in utero arsenic exposure. Specifically, the negative association between miR-1290 and birthweight was stronger when cord blood arsenic exposure was also higher (P-for-interaction <0.05) (Figure 2; Table S3). For instance, among infants delivered at 37th week of gestation, a 1-SD increment in placental expression of miR-1290 was associated with 86.8 g (95% CI: −164.1, −9.5; P = 0.028), 120.8 g (95% CI: −197.8, −43.8; P = 0.002), and 162.9 g (95% CI: −249.9, −75.9; P <0.001) lower birthweight when cord blood arsenic concentrations were at 25th, 50th, and 75th percentiles, respectively.

Figure 1. Directed acyclic graph (DAG) showing conceptual models for causal mediation analysis. We hypothesized that placenta-derived miRNAs will affect birthweight directly or indirectly via changing gestational age.

Figure 2. Controlled Direct Effects (CDEs) of selected placenta-derived miRNAs in relation to birthweight at gestational age ranging from 29 to 41 weeks. Panel A1, A2, and A3 illustrate the changes in mean birthweight [g] for a 1-SD increment in placental expression miR-1290 when cord blood arsenic exposure is at 25th, 50th, and 75th percentiles, respectively. Panel B, C, and D illustrate the changes in mean birthweight [g] for a 1-SD increment in placental expression of miR-195, miR-328, and miR-489, respectively. Solid black lines surrounded by shaded area represent the effect estimates and 95% confidence intervals. The horizontal blue dashed lines show the reference values and the vertical blue dashed lines represent the gestational age of 37 weeks. Models were adjusted for natural log cord blood arsenic exposure, maternal age, education, enrollment BMI, number of past pregnancies, infant sex, study site and surrogate variables.

We conducted a sensitivity analysis to determine the temporal relation between placenta-derived miRNAs and gestational age in relation to birthweight. Our results showed that the effect estimate of gestational age in relation to birthweight was consistent whether or not a miRNA was added to the model. This pattern was consistent for all miRNAs significantly associated with birthweight (Table 3). On the other hand, the effect estimate of a miRNA in relation to birthweight considerably changed whether or not gestational age was added to the model. These results suggest that gestational age is likely an intermediate or mediating variable for the associations between our top-hit miRNAs and birthweight.

In silico target mapping and miRNA-target gene enrichment analysis

Gene enrichment by process networks of candidate miRNAs with their in silico predicted targets revealed several biological processes, such as regulation of cell growth and proliferation, apoptosis, angiogenesis, embryonic development, erythropoiesis, reproduction and neurohormonal signaling, and inflammation that were associated with miR-1290, miR-195, and miR-328 (Table S4). These network processes involved several biological pathways, including growth factor signaling, insulin and IGF family signaling, WNT signaling, VEGF signaling, and TGF signaling cascades that were implicated in fetal growth (Figure 3; Table S4).

Figure 3. Gene enrichment pathways for selected miRNAs with in silico predicted targets. The vertical axis represents a description of the pathways, and the horizontal axis represents the – Log₁₀(P value) of the significant enrichment pathways.

Discussion

In our two-phase study involving a total of 454 mother-infant pairs, placental expression of miR-1290, miR-195, and let-7g was significantly and negatively associated with gestational age, while miR-328 was significantly and positively associated with gestational age (FDR <0.05). Causal mediation analysis revealed that placental expression of miR-1290, miR-195, and miR-27a was significantly and negatively associated with birthweight, while miR-328 and miR-324-5p was significantly and positively and associated with birthweight, primarily via gestational age. The effect of miR-1290, miR-195, and miR-328 on birthweight was modified by gestational age, whereas the effect of miR-1290 was additionally modified by in utero arsenic exposure. A stronger effect of these miRNAs on birthweight was observed among infants delivered early in gestation or exposed to higher concentration of arsenic in cord blood. Overall, our results indicated that a higher expression of miR-1290 and miR-195 and a lower expression of miR-328 in the placenta was not only associated with early termination of pregnancy but also associated with more than additively lower birthweight. Gene ontology enrichment analysis of these miRNAs with in silico predicted targets revealed several biological pathways, including growth factor signaling, insulin and IGF family signaling, WNT signaling, and VEGF and TGF signaling cascades implicated in fetal growth.

We identified several miRNAs (e.g., let-7g, miR-30c, miR-138, miR-181c, and miR-489) that were not previously implicated in fetal growth. We also identified miRNAs that were previously reported to be associated with fetal growth (miR-1290, miR-324, miR-424, miR199a) [21,28–30] or implicated in fetal growth based on their associations with preeclampsia (miR-181a [20,31,32], miR-424 [33,34], miR-27a [35], miR-199a [33], miR-200c [20,34], miR-195 [36,37], and miR-328 [14]). However, most of these studies employed a targeted miRNA expression analysis approach with limited sample size, often resulting in conflicting or irreproducible results. In this two-stage study, we identified several placenta-derived miRNAs associated with birthweight in Phase 1, which were further confirmed in Phase 2, minimizing the chance of false positive results.

We observed most robust associations for miR-1290, miR-195, and miR-328 in relation to birthweight, where, respectively, 33.1%, 46.5%, and 64.8% of the total effect of these miRNAs were mediated via gestational age. Our data showed that gestational age not only mediated but also moderated the effect of these miRNAs on birthweight. Specifically, a stronger effect of these miRNAs on birthweight was observed among infants who were delivered early in gestation (all P-for-interaction <0.05). These three miRNAs also showed significant associations with gestational age, where higher expression of miR-1290 and miR-195 and lower expression of miR-328 was associated with early termination of pregnancy. Future studies are warranted to validate our findings and assess miR-1290, miR-195, and miR-328 as potential biomarkers to predict fetal growth and the risk of early delivery. Studies have shown that placenta-derived miRNAs can be detected in maternal circulation, as they are released from placental tissues into the circulation encased in small vesicles, known as exosomes [15]. A nested case-control study in Mexico reported that the expression of miR-1290 in the extracellular vesicles collected from 2nd trimester maternal blood was negatively associated with birthweight-for-gestational age Z-score [21]. However, the study was exploratory in nature and the observed association was not confirmed in a validation cohort. We also found significant positive association for miR-720 in relation to birthweight; however, this miRNA has recently been identified as a fragment of a transfer RNA (tRNA) and subsequently removed from the miRbase database [38].

The mechanisms how placenta-derived miRNAs may regulate birthweight are not well understood. Evidence from observational and experimental studies indicates that fetal hypoxia plays a critical role in fetal growth [39–41]. Hypoxia triggers oxidative stress by generating oxygen-derived free radicals. Placental oxidative stress is considered to be one of the major factors responsible for placental dysfunction and remodeling [42]. Interestingly, of the 49 miRNAs we selected in Phase 1 for further analysis based on their associations with gestational age and/or the extreme birthweight groups, 13 miRNAs (let-7g, miR-15a, miR-16, miR-20b, miR-21, miR-27a, miR-181a, miR-181c, miR-195, miR-199a, miR-320, miR-373, and miR-424) were reported to be regulated by hypoxia [43,44]. Six of these miRNAs (miR-27a, miR-181a, miR-181c, miR-195, miR-199a, miR-424) showed significant associations with birthweight and/or gestational age at the level of P <0.05, while miR-27a and miR-195 showed significant indirect effect associations with birthweight at the level of FDR <0.05 in Phase 2 analysis.

miR-195 is a member of the miR-15 gene family, along with miR-424 and miR-15a, and has been reported to be dysregulated in response to cellular hypoxia in experimental studies [44,45]. For example, in ATDC 5 cell, a resident cartilage cell adapted to function in chronic hypoxic condition, miR-195 was downregulated, coupled with a reciprocal upregulation of its direct target, hypoxia initiation factor 1 alpha (HIF-1α)[13]. HIF1α is a transcription factor and considered to be the master regulator of cellular hypoxia. HIF1α has been linked with transcription of many genes, including genes that encode erythropoietin (EPO), vascular endothelial growth factor (VEGF), and glucose transporter 1 (GLUT1) [46]; the protein products of these genes are responsible for oxygen and nutrient transfer and glucose entry across cell membrane [47–49]. Therefore, downregulation of miR-195 in the placenta could potentially lead to increased oxygen and nutrient supply to the fetus via upregulation of HIF1α in the placenta and vice versa. Observational and experimental studies have also reported downregulation of FLT3, BCL-2, M-Ras, and A-Raf-1 in relation to higher cellular expression of miR-195 (Figure 4); the protein products of these genes are responsible for cell growth and differentiation, apoptosis, proliferation, and bone marrow hematopoietic cell differentiation [50–52].

Figure 4. Networks of miR-1290 and miR-195 with experimentally validated target genes. The green lines connecting the nodes represent activation, whereas the red lines connecting the nodes represent inhibition. The arrow indicates the direction of the association.

On the other hand, experimental studies have validated that miR-1290 has binding sites on MT1G, MT1H, CYP4F11, NAT-1 genes, and thereby downregulates the expression of these genes [53,54]. The protein product of these genes are responsible for arsenic and folate metabolism [55], metal transport and elimination [53], and cellular protection from free radical-induced oxidative damage [56]. Therefore, higher expression of miR-1290 in the placenta may result in increased arsenic toxicity. In this study, we observed that the effect of miR-1290 on birthweight was stronger when the arsenic concentration was higher in cord blood. In utero arsenic exposure has been associated with increased risk of preterm birth [57] and lower birthweight [27] in this population. Therefore, higher expression of miR-1290 in the placenta could potentially influence birthweight via shortened gestational age and increased arsenic toxicity. Future studies are needed to experimentally validate these postulations.

Our study has some limitations. Study participants selected in Phase 1 and Phase 2 were not completely independent: 52 participants were included in both study phases. In Phase 1, we implemented an extreme phenotype sampling design and selected two extreme birthweight groups from both tails of the birthweight distribution. In Phase 2, we selected a random sample of the full cohort without excluding samples that were selected in Phase 1. While this sampling strategy limited our ability to replicate findings in an independent cohort, it gained us predictive power in Phase 1 (because of the extreme phenotype design) and protected us from penalization for any loss of variance in birthweight distribution in Phase 2 analyses (because of not excluding participants with extreme birthweight). It also enabled our Phase 2 results to have a population-level inference. We conducted a sensitivity analysis excluding the 52 duplicate samples from our Phase 2 cohort and found consistent results in relation to gestational age (Table S5). However, the association between miRNA expression and birthweight was largely attenuated in terms of both the strength of association and statistical significance for most miRNAs. This could be due to the loss of variance in birthweight distribution in our Phase 2 cohort in addition to the loss of power because of a loss of sample size. We encourage validating our findings in external population to address this limitation. Other limitations of the study include the lack of epigenome-wide coverage of our selected miRNA expression platform in Phase 1, which likely have limited our ability to identify some miRNAs, especially the novel ones.

Our study has several strengths. The two-stage study design allowed us to validate our Phase 1 results in the Phase 2 cohort, reducing the risk of false discovery. In a two-stage design, however, technical variations may arise when different expression platforms are used, which potentially could introduce bias if not accounted for during the data analysis stage. In both phases of our study, we used miRNA expression platform based on qRT-PCR technology, which is considered as the ‘gold standard’ for small RNA quantification. Furthermore, we constructed surrogate variables and included those variables in all phases of data analysis as covariates to account for known or unknown, biological or technical variations that are common in high throughput gene expression studies [58]. Furthermore, epigenetic signatures are tissue-specific; it is important to measure these signatures in the tissue biologically most relevant to the exposure and outcome of interest. The placenta is likely the organ to express the strongest signal of any in utero perturbations related to fetal growth. Moreover, we examined the link between placenta-derived miRNAs and birthweight accounting for environmental exposure during the in utero period, when the critical window is narrow, and the fetus is most susceptible to environmental exposures.

Our results suggested that some placenta-derived miRNAs were significantly associated with birthweight, and that the effects of these miRNAs were primarily mediated via changing gestational age. The effect of placenta-derived miRNAs on birthweight was also modified by gestational age and in utero arsenic exposure, where a stronger effect was observed among infants delivered early in gestation or exposed to a higher level of arsenic in cord blood. Future studies are warranted to replicate these findings in external population and evaluate placenta-derived miRNAs as early biomarkers of fetal growth.

Materials and methods

All protocols were reviewed and approved by the Human Research Committees at Harvard T.H. Chan School of Public Health and Dhaka Community Hospital Trust. All participants provided written informed consent before participating in the study.

Study population and subject selection

The study population was derived from a prospective birth cohort established in Pabna and Sirajdikhan Upazilas of Bangladesh. During 2008–2011, we enrolled 1,613 pregnant women aged ≥ 18 years with an ultrasound-confirmed singleton pregnancy of ≤ 16 weeks of gestation and followed them throughout pregnancy. Details of this study, including recruitment and enrollment criteria were previously reported [25,27]. At the end of follow-up, 1,184 singleton live births were recorded after exclusions due to loss of contact before delivery (n = 99), participation withdrawal (n = 121), miscarriage (n = 132), stillbirth (n = 72), and multiple pregnancies (n = 5). Placental tissue samples were collected from 1,141 participants.

We applied a two-phase study design: Phase 1 and Phase 2, with a total of 466 mother-infant pairs without a history of maternal diabetes, gestational hypertension, or preeclampsia. More details of the study design and analytical approaches appear in Table 4. In Phase 1, we applied extreme-phenotype sampling design and selected participants from both tails of the birthweight distribution. We selected 77 participants with the lowest birthweights (birthweight range: 0.8–2.2 kg) and 77 participants with the highest birthweights (birthweight range: 3.3–3.9 kg) after excluding 24 participants with low RNA quality and 8 participants with birthweight >4,000 g. This design was chosen to maximize the predictive power in Phase 1. Top-hit miRNAs screened out in Phase 1 were further analyzed in Phase 2 in a cohort of 364 participants randomly selected from the full cohort, without removing participants already selected in Phase 1. Thus, 52 participants with extreme birthweights who were already selected in Phase 1 were also selected in Phase 2.

Table 3. Effect estimates of gestational age and selected miRNAs in relation to birthweight.

	Effect estimates of gestational age		Effect estimates of miRNA expression
miRNAs	Model 5^a	Model 7^c	Model 6^b	Model 7^c
let7g	78.7	78.9	−28.0	9.0
miR-1290	78.7	74.0	−121.1	−68.7
miR-138	78.7	77.1	−63.1	−50.0
miR-181a	78.7	78.1	−50.2	−46.1
miR-181c	78.7	77.8	−64.7	−57.5
miR-195	78.7	74.5	−98.3	−56.1
miR-199a	78.7	80.9	−35.3	−69.7
miR-27a	78.7	79.7	−8.9	27.9
miR-324-5p	78.7	78.6	32.6	1.1
miR-328	78.7	76.7	49.1	16.1
miR-424	78.7	79.7	7.8	−19.3
miR-489	78.7	77.1	−56.8	−43.0
miR-720	78.7	76.7	81.8	71.9

^a Model 5: $E\left[Birthweight\right]={\theta }_0+{\theta }_1\left(GestationalAge\right)+\theta {\mathrm{\prime }}_{2}Covariates$

^b Model 6: $E\left[Birthweight\right]={\beta }_0+{\beta }_1\left(miRNA\right)+\beta {\mathrm{\prime }}_{2}Covariates$

^c Model 7: $E\left[Birthweight\right]={\alpha }_0+{\alpha }_1\left(miRNA\right)+{\alpha }_2\left(GestationalAge\right)+\alpha {\mathrm{\prime }}_{4}Covariates$

Coefficients were scaled for a 1-SD change in miRNA expression in the placenta. Covariates included in the models were natural log cord blood arsenic exposure, maternal age, education, enrollment BMI, the number of past pregnancies, study site, infant sex and surrogate variables.

Table 4. Study design and analysis work flow.

Phase 1 (n = 154)
Study Design	Extreme Phenotype Sampling: 77 samples with birthweight ranging 0.8 – 2.2 kg 77 samples with birthweight ranging 3.3 – 3.9 kg
Expression platform	TaqMan RT-qPCR OpenArray for 754 preconfigured human miRNAs
Normalizations	Endogenous control; 2) Quantile; 3) Global; 4) Raw Ct
Association analyses	Associations between miRNA expression and gestational age Associations between miRNA expression and the extreme birthweight groups with and without adjusting for gestational age
MicroRNA selection	Total 49 miRNAs were screened out for further analysis in Phase 2
Phase 2 (n = 364)
Study Design	364 participants randomly selected from the full cohort (n = 1141) 312 new participants; 52 duplicates also selected in from Phase 1
Expression platform	TaqMan RT-qPCR Customized OpenArray Platform 310 samples in triplicate; 54 samples in duplicate
Normalization	Endogenous control
Association analyses	Associations between miRNA expression and gestational age Direct, indirect, and total effect associations between miRNA expression and birthweight considering gestational age a mediator Assess interactions between miRNA expression and gestational age as well as miRNA expression and arsenic exposure miRNA gene enrichment analysis with in silico predicted mRNA targets

Outcome and covariates

We collected participant’s socio-demographic and medical history, including maternal age, education, and the number of past pregnancies using a structured questionnaire. Trained healthcare workers measured maternal height and weight at the time of enrollment. Total inorganic arsenic was measured in cord blood using a dynamic reaction cell inductively coupled plasma mass spectrometer (DRC-ICP-MS, DRC II, Perkin Elmer) and subsequently transformed to their natural log for analysis. Gestational age was determined by ultrasonography using crown-rump length at median 12 (range: 4–16) weeks of pregnancy. Birthweight was measured on a pediatric scale calibrated before each measurement and rounded to the nearest 10 g.

MicroRNA expression in the placenta

Placental tissue was sampled from the fetal side of the placenta within 2 hours after delivery. Two placental tissue specimens, one from the inner region proximal to the umbilical cord insertion point and one from the outer region closer to the edge, were collected to control for within-placenta variability. Methods and rationale for the sampling scheme were described previously [59]. Tissue specimens were preserved in RNAlater (Ambion) to stabilize RNA and then stored at −80°C until analysis.

Total RNA, containing small RNA, was extracted from 50 mg of thawed and homogenized placenta using Qiagen miRNeasy kit (Qiagen) following manufacturer’s guidelines. RNA purity and concentration were evaluated by spectrophotometry using NanoDrop ND-2000 (ThermoFisher). RNA samples with an A260/230 ratio between 1.8–2.2 and a minimum concentration of >38 ng/μl (range: 38–687) were used to quantify miRNA expression. The quality and related size of total RNA, including small RNA was assessed by Agilent 2100 Bioanalyzer microfluidics-based platform (Agilent Technologies). RNA was transcribed into cDNA in 2 multiplex reactions, each containing 3 μL of the small RNA preparation and either TaqMan Megaplex RT primers, Human Pool A, or Pool B, and then preamplified using TaqMan Megaplex PreAmp Primers, Human Pool A, and Pool B. A total of 200 ng total RNA, containing the small RNA fraction, was used for each sample for reverse transcription with Megaplex^TM.

In Phase 1, we used TaqMan Human MicroRNA assay (Applied Biosystems) to quantify 754 preconfigured human miRNAs. This system can perform multiple simultaneous real-time polymerase chain reaction (qRT-PCR) running TaqMan assays. In Phase 2, we used TaqMan qRT-PCR based QuantStudio 12K Flex OpenArray platform (Applied Biosystems) to confirm the expression of top-hit miRNAs screened out in Phase 1. A total of 54 participants out of 364 Phase 2 participants were analyzed in duplicate, while 310 participants were analyzed in triplicate. Technical duplicates were used to assess confounding and bias from technical variations within and across the phases. The quality of qRT-PCR data between Phase 1 and Phase 2 appears in Table S6.

Statistical analysis

Phase 1 analyses

The demographic characteristics between the two extreme birthweight groups were assessed by ANOVA, Chi-square tests, or Fisher’s exact test, and a two-tailed P value was reported. Data were checked for outliers using box plots, and samples with low RNA yield or absorbance were excluded. One sample with very few detectable miRNAs was also excluded from analysis. microRNA expression data were quality controlled using the following criteria: cycle threshold (Ct) <30, amplification score > 1.1, Cq confidence >0.8, and detectable expression in ≥50% samples in the extreme birthweight groups. In Phase 1, 306 miRNAs passed the quality control criteria and were used for statistical analyses. To minimize the risk of false-positive results in the evaluation of 306 miRNAs due to multiple testing, we used an arbitrary conservative cut-off of permutation P value <0.01 to determine statistical significance. Because miRNA expressions are highly correlated, and Phase 1 study was designed for screening purpose, we did not use stringent criteria to define statistical significance in this phase. To estimate permutation-based P value, we conducted 1,000 permutation tests for each miRNA and randomly paired the outcome status with each miRNA within each study participant. Observed P values were estimated for each miRNA 1,000 times using permutation simulations. Permutation-based P value was then estimated for each miRNA given the distribution of 1,000 observed P values.

Imputation was used to assess missing value effect. Two imputation strategies were applied in Phase 1. Assuming that Ct values were not missing at random but missing due to no or low expression, we imputed the missing value based on the 95th percentile of that miRNA among the extreme birthweight groups without missing value. We also used multivariate imputation by chained equation (MICE) [60] under the assumption that missing values were at random as a result of a failed measurement. Test statistics for both imputed strategies were compared with that obtained in complete case analysis. Surrogate variables were constructed and were included in all association models to adjust for known/unknown biological or technical variations (e.g., batch effects) [58]. Endogene normalization was applied as primary data normalization technique. The geometric mean of 3 most stable internal control miRNAs (miR-364, miR-532, miR-532-3p) was used to normalize miRNA expression data following the method described by Vendesompele et al. [61]. To test the sensitivity of our results due to technical variations, we also analyzed the data on raw Ct values as well as after implementing global and quantile normalization in Phase 1 using R/Bioconductor.

In Phase 1, we screened out miRNAs associated with the extreme birthweight groups and/or gestational age using three sets of models (Model 1–3) based on our hypothesis that gestational age will mediate the effect of miRNA expressions on birthweight. Multivariate linear regression was applied to assess the association between miRNA expression and the extreme birthweight groups using miRNA expression as the dependent variable, adjusting for maternal age (continuous), infant sex (male, female), study site (Pabna vs. Sirajdikhan), and with (Model 1) and without gestational age (Model 2). Multivariate linear regression was also applied to assess the association between miRNA expression and gestational age using gestational age as the dependent variable, adjusting for maternal age, infant sex, study site, and the extreme birthweight groups (Model 3). All models were additionally adjusted for surrogate variables previously constructed. Statistical analyses were applied separately for each miRNA for each normalization method.

(1) $Y{\left(miRNA\right)}_{ij}={\theta }_0+{\theta }_1{\left(Birthweightgroups\right)}_i+{\theta }_2{\left(GestationalAge\right)}_i+{\theta }_3{\mathbf{C}}_i+e_i$(1)

(2) $Y{\left(miRNA\right)}_{ij}={\theta }_0+{\theta }_1{\left(Birthweightgroups\right)}_i+{\theta }_2{\mathbf{C}}_i+e_i$(2)

(3) $Y{\left(GestationalAge\right)}_i={\beta }_0+{\beta }_1{\left(miRNA\right)}_{ij}+{\beta }_2{\left(birthweightgroups\right)}_i+{\beta }_3{\mathbf{C}}_i+e_i$(3)

Top-hit microRNAs were selected for further analysis in Phase 2 if the following criteria were met: 1) a miRNA showed statistically significant association with either the extreme birthweight groups or gestational age in at least 2 out of 4 normalization methods at the level of permutation P <0.01; or 2) a miRNA showed statistically significant association with either the extreme birthweight groups or gestational age in any normalization method at the level of permutation P <0.01, but has known biological significance.

Phase 2 analyses

In Phase 2, we applied the same data quality control criteria as Phase 1. We calculated mean Ct value of each miRNA for each sample. Endogenous normalization was applied taking the geometric mean of the same three internal control miRNAs (miR-364, miR-532, miR-532-3p). Multivariate imputation by chained equation (MICE) was used to predict three missing values in miRNA expression [60]. Surrogate variables were estimated and included as covariates in all association models [58]. We used a significance level of false discovery rate (FDR) <0.05 to account for multiple testing.

The associations between miRNA expression and birthweight were analyzed using causal mediation analysis technique considering gestational age a mediator [62]. The directed acyclic graph (DAG) in Figure 1 illustrates our conceptual model [63]. We fit two sets of multivariate linear regressions models: the model for gestational age (mediator model, Model 4) and the model for birthweight (outcome model, Model 5).

(4) $E\left[GestationalAge\right]={\beta }_0+{\beta }_1\left(miRNA\right)+{\beta }_2\mathrm{\prime }C$(4)

(5) $E\left[Birthweight\right]={\theta }_0+{\theta }_1\left(miRNA\right)+{\theta }_2\left(GestationalAge\right)+{\theta }_3\left(miRNA\ast GestationalAge\right)+\theta {\mathrm{\prime }}_{4}c$(5)

We selected a priori the following list of covariates that were previously found to be associated with birthweight: natural log of cord blood arsenic exposure (continuous), maternal age (continuous), education (no formal education, primary, secondary or higher), enrollment BMI (<18.5, 18.5–24.9, >24.9), number of past pregnancies (0, 1, >1), infant sex, and study site in addition to the surrogate variables previously constructed to adjust for both the mediator and outcome models. The model for birthweight was additionally adjusted for gestational age. The natural direct, indirect effect, and total effects of mediation analysis were estimated combining the mediator and the outcome models (Supplemental file). Analyses were conducted with and without an interaction term between miRNA expression and gestational age. For miRNAs with statistically significant interaction with gestational age, controlled direct effect of that miRNA on birthweight was estimated at the gestational age ranging from 29 to 41 weeks. We also tested for interaction between miRNA expression and cord blood arsenic exposure in relation to birthweight. A significant interaction was present for miR-1290, but not for other miRNAs; hence, a miRNA-arsenic interaction term was added to the model for birthweight for miR-1290.

The direct and indirect effects were averaged across all study participants. Bias-corrected standard errors were computed from 10,000 bootstraps. Mediation analysis assumes that conditional on the covariates, there is no confounding of the exposure-outcome relation, exposure-mediator relation, or mediator-outcome relation and that there are no mediator-outcome confounders affected by the exposure [62]. The R-package ‘mediation’ was used to conduct causal mediation analysis [64]. The ‘medsens’ function of the ‘mediation’ package was used to conduct sensitivity analysis to test for unmeasured confounding assumptions for direct and indirect effects of mediation analysis.

The expression of miRNAs in the placenta was rescaled to their standard deviation (SD) to aid in interpretation and comparison of point estimates across miRNAs. Thus, effects estimate for each miRNA was presented as the mean change in the outcome of interest (gestational age or birthweight) for a 1-SD increment in the expression of respective miRNA in the placenta. The average expression of a miRNA in the placenta in full-term pregnancies (40 weeks) was considered as the baseline expression for that miRNA under the normal physiological condition, assuming that the expression of placenta-derived miRNAs temporarily precedes gestational age. To test whether our assumption of the temporal relation between miRNA expression and gestational age in relation to birthweight was valid, we conducted a sensitivity analysis using three sets of models (Model 5–7).

(6) $E\left[Birthweight\right]={\theta }_0+{\theta }_1\left(GestationalAge\right)+\theta {\mathrm{\prime }}_{2}C$(6)

(7) $E\left[Birthweight\right]={\beta }_0+{\beta }_1\left(miRNA\right)+\beta {\mathrm{\prime }}_{2}C$(7)

(8) $E\left[Birthweight\right]={\alpha }_0+{\alpha }_1\left(miRNA\right)+{\alpha }_2\left(GestationalAge\right)+\alpha {\mathrm{\prime }}_{4}C$(8)

We compared the effect estimate of gestational age in relation to birthweight between models adjusted for miRNA expression and not adjusted for miRNA expression (Model 7 vs. Model 5), in addition to previously selected confounding variables (C). Similarly, we compared the effect estimate of a miRNA in relation to birthweight between models adjusted for gestational and not adjusted for gestational age, in addition to the same set of covariates (Model 7 vs. Model 6). If miRNA expression were an intermediate variable for the association between gestational age and birthweight, we would see a change in the effect estimate of gestational age in relation to birthweight between Model 7 vs. Model 5. Similarly, if gestational age were an intermediate variable for the association between miRNA expression and birthweight, we would see a change in the effect estimate for a specific miRNA in relation to birthweight between Model 7 vs. Model 6. We conducted this analysis for each miRNA that was significantly associated with birthweight.

We used MetaCore (GeneGo) bioinformatics tool to identify experimentally validated mRNA targets of top-hit miRNAs and constructed miRNA-mRNA interaction networks. MetaCore was also used to conduct gene enrichment analysis by process networks and pathway maps of top-hit miRNAs with in silico predicted targets with a TargetScan total context score <-0.2 [TargetScan v7.1]). All analyses were conducted in R, version 3.2.2 (Austria, Vienna) and STATA, version 14.0 (College Station, TX).

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

This work was supported by the National Institute of Environmental Health Sciences (NIEHS) of the National Institutes of Health (NIH) grant R01ES015533.