Placental microRNAs associate with early childhood growth characteristics

Poor placental function is a common cause of intrauterine growth restriction, which in turn is associated with increased risks of perinatal morbidity, mortality and long-term adverse health outcomes. Our prior work suggests that birthweight and childhood obesity-associated genetic variants functionally impact placental function and that placental microRNA are associated with birthweight. To address the influence of the placenta beyond birth, we assessed the relationship between placental microRNAs and early childhood growth. Using the SITAR package, we generated two parameters that describe individual weight trajectories of children (0-5 years) in the New Hampshire Birth Cohort Study (NHBCS). Using negative binomial generalized linear models, we identified placental microRNAs that associate with growth parameters (FDR<0.05), while accounting for sex, gestational age at birth, and maternal parity. Genes targeted by the six growth trajectory-associated microRNAs are enriched (FDR<0.02) in growth factor signaling (TGF/beta: miR-1290; EGF/R: miR-155, Let-7c; FGF/R: miR-155; IGF/R: Let-7c, miR-155, miR-1290), cyclic AMP signaling (miR-1246), calmodulin signaling (miR-216a, miR-1246), and NOTCH signaling (miR-629). These pathways function in placental proliferation, differentiation and function. Our results support the hypothesis that fetal environment, specifically placental cellular dynamics and function guided by microRNA expression, can have impacts beyond birth, into early childhood.


INTRODUCTION 1
The Developmental Origins of Health and Disease (DOHaD) posits that in utero exposures can 2 induce permanent maladaptive changes that "program" the fetus to have increased risk of later-3 life diseases 22 . The placenta is the master regulator of the intrauterine environment and an ideal 4 tissue in which to test the DOHaD paradigm. This ephemeral organ facilitates the exchange of 5 gases, nutrients and waste for the developing embryo 23,24 , in addition to aiding nutrient 6 metabolism, and acting as an endocrine organ critical for early development. Thus, its functions 7 are responsible for proper development and programming of the offspring. 8 MicroRNAs are 21-25 base pair non-coding RNAs that regulate gene expression post-9 transcriptionally via sequence complementarity to the 3' untranslated region of mRNA 10 transcripts. microRNA-mediated gene regulation is achieved through target mRNA translational 11 inhibition or degradation 25 . MicroRNAs are proposed to regulate more than 50% of human genes 12 given gene are over-represented among loci that associate with childhood obesity and BMI 29 . 24 Compared to eQTLs for seven adult tissues, placental eQTLs were most strongly enriched 25 among results from GWASs of birthweight, childhood obesity, and childhood BMI. These 26 results demonstrate that the placental transcriptional landscape can have a lasting impact on early 27 childhood growth and metabolism 29 . 28 In this study, we pose the question: Do microRNAs that regulate the placental mRNAs 29 have a lasting impact on childhood growth? To address this question, we model childhood 30 weight from birth to five years of age in the New Hampshire Birth Cohort Study (NHBCS). 31 Using parameters derived from individual growth curves, we relate placental microRNA 32 expression, from small-RNA sequencing, with early childhood growth. We use in silico methods 33 to predict mRNA targets of interesting microRNAs and pathway analysis to add biological 34 context to our findings. 35

MATERIALS AND METHODS 36
Cohort 37 The New Hampshire Birth Cohort Study. NHBCS was initiated in 2009 and is an ongoing 38 study comprised of a cohort of mother-infant pairs. Pregnant women between 18 and 45 years of 39 age were recruited from the study's participating prenatal care clinics in New Hampshire, USA. 40 Women were included in the cohort if their primary source of drinking water was from an 41 unregulated residential well, they had resided in the same household since their last menstrual 42 period and had no plans to move before delivery. All participants provided written informed 43 consent in accordance with the requirements of the Committee for the Protection of Human 44 Subjects, the Institutional Review Board (IRB) of Dartmouth College. In this study, NHBCS 45 participants were singleton pregnancies recruited between February 2012 and September 2013. 46 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Statistical analyses 93
Modeling childhood growth trajectory. The R package SITAR was used to model growth 94 trajectory for all study participants with at least three weight observations (n = 238) 36 . Individual 95 children's growth curves were visually assessed for outlier measures. Three observations were 96 removed (Fig. S1). The children's weights were natural log transformed to meet the normality 97 assumption of the SITAR method. Briefly, SITAR is a shape invariant model, with individual 98 random effects for each child, that estimates an average growth curve across samples. A set of 99 three parameters are generated for each individual that, through translation (vertical/horizontal 100 shift) and rotation (counter/clockwise), transform the average growth curve to match each 101 individual's growth. The parameters, random effects, have mean zero and standard deviations 102 estimated from the data. These parameters -size, intensity and tempo, are interpreted as 103 percentiles relative to the average because the weights were natural log transformed. The size 104 parameter represents average size for any child relative to the average child and is graphically 105 represented as a vertical shift of the weight curve (Fig. 2). Tempo, the age at peak weight 106 intensity, was not estimated for this analysis, since no measures were estimated around the peak 107 of infant growth (6 weeks); intensity is expressed as a percentage deviation from mean 108 intensity with higher values representing faster growth than average. Graphically, intensity 109 represents a rotation of the weight curve (Fig. 2). We modeled the transformed weights using 10 110 degrees of freedom, and adjusted for child sex in the models that generated the size and intensity 111 parameters. 112 Surrogate variable analysis. To adjust for batch effects, cell-type heterogeneity and other 113 unknown sources of technical variation, we estimated a surrogate variable from the normalized 114 transcript reads via the svaseq function in the sva package that incorporates the Combat 115 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) and gestational age (z-score). Three samples were missing information for parity and were 124 removed from the intensity analysis (n = 235). The size analysis included covariates for 125 primigravida and maternal educational (beyond high school). Five samples were missing 126 information for the maternal education variable and were removed from the size analysis (n = 127 230). One surrogate variable was also included as a covariate in both regression models 128 (discussed in previous section). We considered microRNAs with a false discovery rate less than 129 5% to be differentially expressed with either size or growth intensity. 130 Sensitivity analysis. To assess potential confounding of DEmiR effects by pre-eclampsia, 131 maternal pre-pregnancy BMI, gestational weight gain or birthweight, potential confounders were 132 individually added to the DEA model in DESeq2. Effect estimates and standard errors were 133 collected for intensity and size DEmiRs for comparison with original model statistics. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

151
This study analyzed data from 238 mother-infant pairs from the New Hampshire Birth 152 Cohort (NHBCS). Parallel placental microRNA transcript abundance was available for this 153 cohort, as well as up to 14 weight measures between birth and 5 years, collected from well-child-154 checks. Study participants were excluded from the analysis if they had less than three separate 155 weight observations, but were not excluded for intermediate missing observations (Fig. 1). 156 Observations as a fraction of participants in the study at each measurement remained high 157 throughout the study, despite attenuated participation with time, especially after one year (Fig.  158 1). 48% of participants had at least 12 of the 14 potential observations and only 2% had less than 159 5 observations (median observations across participants is 11, ranged 3-14). 160 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint SuperImposition by Translation And Rotation (SITAR) 36 . SITAR provided two parameters to 162 describe each child's growth. The size parameter is given as a percent, relative to the fitted 163 weight curve for all participants, for the average size of each individual child. Higher size is 164 graphically represented as an upward vertical shift in a child's weight growth curve. The size 165 parameter in this analysis ranged between -37% and 31% (Fig. 2, red and blue lines). The 166 intensity parameter is expressed as a percentage deviation from the mean intensity for all 167 participants, with higher values representing faster growth than average. Higher growth intensity 168 is graphically represented as a counterclockwise rotation in a child's weight growth curve. The 169 intensity parameter for this analysis ranged between -9% and 12% (Fig. 1, purple and cyan lines). 170 The demographics of the participants for this study are shown in Table 1 Covariates were included in the model if they were associated with the variable of interest 181 (intensity or size) in univariate regression analyses (Table 2). For the intensity analysis, models 182 also included maternal parity, gestational age (z-scores) and one surrogate variable (to account 183 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint for batch effects, cell type heterogeneity, and other sources of unknown variability). For the size 184 analysis, primigravida, maternal educational attainment and one surrogate variable were included 185 in the models. Of the 777 microRNA transcripts that passed quality control filtering for the 186 analyses, five and one had an FDR < 0.05 for the intensity and size analyses, respectively (p-187 value < 1.4x10 -4 , Fig.s 3, S2). For the intensity and size analyses, respectively, two and one of 188 those microRNAs had p-values less than the Bonferroni family-wise error rate threshold (p-value 189 < 6.4x10 -5 , Fig.s 3, S2). it was not included in the primary analysis. However, the potential for pre-eclampsia and other 203 potential confounders (pre-eclampsia, maternal pre-pregnancy BMI, gestational weight gain and 204 birthweight) to attenuate significant associations (Fig. 3) between placental microRNAs and 205 growth trajectory parameters was assessed in sensitivity analyses (Fig. S3). As expected for a 206 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. For each DEmiR, putative target genes were tested for pathway overrepresentation with 225 consensuspathdb (CPDB), against all of the genes that passed QC in the RICHS whole 226 transcriptome RNA-seq analysis 41 . The most significant pathways are listed in Table 3. CPDB 227 integrates pathways from 12 databases, meaning that similar results can be reported across 228 sources. To further prioritize significantly overrepresented pathways, we generated term 229 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint frequency matrices for all pathways for which microRNA targets made up more than 20% of 230 pathway genes and for which the enrichment had an FDR less than 0.02. To assess the 231 microRNA target genes leading the enrichments, a second term frequency matrix was made from 232 the microRNA target genes present in the significant pathways. From these analyses, we found 233 that FGFR1-4 and IGF/R signaling pathways were overrepresented among significant miR-155 234 pathways. The most frequently involved miR-155 target genes were PIK3CA and PIK3R1. For 235 miR-1290, BMP and TGF/beta signaling were the most frequent pathways, with MAPK1 and 236 SMAD2-4 the most frequent target genes. EGF and EGFR pathways were the most frequent 237 among let-7a significant pathways, with MAPK1, PIK3CA and PIK3R1 as the most frequently 238 involved gene targets. There was not repetition among the 13 significant miR-1246 pathways, 239 though the most frequent genes were ADCY1 and PRKACB. Among mir-629 pathways, PTEN 240 was the most frequent pathway, with TNRC6C and AGO1 and AGO4 were the most common 241 gene targets. Finally for the only size DEmiR, miR-216a, regulation of MECP2 was most 242 frequent among pathways, with miR-216a target CALM1-3 as the most frequent genes. When 243 targets of positively-associated growth intensity DEmiRs (Fig. 3) were pooled for pathway 244 analysis, we found that the most frequent pathway terms were FGFR1-4 and the most frequent 245 gene targets were PIK3CA and PIK3R1. 246

DISCUSSION 247
In this study, we have described the microRNAs from human placenta that associate with growth 248 trajectory from birth to five years old. We used a shape invariant model with random effects to 249 generate two parameters that describe children's growth intensity and average size during the 250 observation period. We found evidence of microRNAs that vary with both growth parameters 251 (DEmiRs). We narrowed bioinformatically predicted microRNA targets to only those in which 252 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. and IGFR (let-7c, miR-155, miR-1290) signaling pathways are influenced by growth trajectory 266 microRNAs (Fig. 4). 267 Multipotent villous cytotrophoblasts fuse to form the epithelial, multinucleated 268 syncytiotrophoblast, which facilitates exchange of gases, nutrients and waste for the growing 269 embryo, as well as acting as an endocrine organ to balance the needs of the mother and fetus 48 . 270 Syncytialization is mediated by cAMP and rising intracellular Ca+, as well as continued EGF 271 signaling and Activin A 46,49 . In this analysis, growth trajectory microRNAs are predicted to 272 target components of the cAMP (miR-1246), calmodulin (miR-216a, miR-1246), EGFR (miR-273 155, let-7c) pathways, as well as Activin A (miR-1290; Fig. 4). 274 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint form the invasive extravillous trophoblast that connects the placenta with the maternal decidua 276 and myometrium to establish maternal blood supply 48 . Signaling through NOTCH1/2 and TGFb 277 family members (BMP2, Activins) mediate differentiation to the invasive extravillous 278 cytotrophoblast 50-52 . The role of the TFG-beta superfamily in trophoblast differentiation is 279 complicated, as many of the pathway mechanisms overlap 53 , yet TGF-beta1/2/3 53 and NODAL 54 280 inhibit the invasive phenotype and favor syncytialization while BMP2 and Activin A/B/AB all 281 favor the invasive phenotype 50,51 . In this analysis, we found evidence that growth trajectory 282 associated microRNAs target components of the NOTCH1/2 (miR-629) and TGFb superfamily 283 (miR-1290) signaling cascades (Fig. 4). 284 Although these processes described above were primarily described in first trimester 285 placentae and cell culture, the placental cytotrophoblasts continue to maintain the 286 syncytiotrophoblast and extravillous trophoblasts to term 45 . The discovery of trophoblast 287 progenitor cells in term placenta suggests that some of the signaling pathways relevant in early 288 gestation, remain important later in gestation. Within this assumption, our findings suggest that 289 microRNAs influence the cellular dynamics of placenta at term, which in turn regulate the 290 plasticity and efficiency of the placenta to support fetal growth. This initial growth primes the 291 neonate's metabolism for early growth patterns. However, it is also likely that some of the 292 growth factors under microRNA influence have specific functions at term that are in addition to 293 or exclusive of their functions in early gestation. For instance, while EGF signaling guides 294 trophoblast maintenance, proliferation and differentiation in early gestation, it may stimulate the 295 release of hormones to the maternal and fetal circulations at term 44 . 296 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint been previously associated with birthweight in the NHBCS 1 including miR-155-5p, mir-629-3p, 298 let-7c and miR-1246. Birthweight was not included in our models, because it is an upstream 299 measure of weight, and included as the initial weight in the growth trajectory analysis. 300 Interestingly, the association of miR-1290 with growth trajectory is not sensitive to the inclusion 301 of birthweight to the model, suggesting that placental miR-1290 relates to early childhood 302 growth independently of birthweight. 303 Growth trajectory-associated microRNAs have also been associated with placental 304 characteristics in other studies. In cultured extravillous trophoblasts, miR-155 inhibits cell 305 proliferation by down-regulating cyclin D1/p27 55 . In cultured cytotrophoblasts, Let-7c is 306 associated with reduced proliferation potential and syncytialization 56,57 . Let-7c is associated with 307 the WNT/b-catenin signaling pathway in other progenitor cell types 58 . In cultured extravillous 308 trophoblasts, miR-1290 promotes rearrangement of maternal endometrium via placental 309 exosomes 59 . miR-1246 has been associated with syncytialization and targeting inhibitors of the 310 WNT/b-catenin signaling pathway 60 . 311 To our knowledge, this is the first study to examine the relationships between human 312 placental microRNAs and early childhood growth trajectory. Our results suggest that placental 313 microRNAs effect signaling cascades central to trophoblast proliferation, differentiation and 314 function. Most importantly, our results underscore the importance of placental function and the 315 intrauterine environment in establishing early growth trends. 316 Our findings should be interpreted within the context of this study's limitations. This is 317 an observational study in which RNA was assayed from term placentae. Thus, we cannot 318 conclude that our results are representative of microRNA associations throughout development. 319 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint where childhood growth was determined using child-specific estimates of size and intensity 321 derived using SITAR. In DESeq2, microRNA abundances (dependent variable) are regressed on 322 size or intensity (independent variable), allowing us to examine the association of early 323 childhood growth on placental microRNA expression. Although the models employed in 324 DESeq2 are temporally reversed in terms of what is being modeled as dependent and 325 independent variables, DESeq2 represents the best choice for modeling gene expression data as it 326 accounts for the over-dispersion of RNA-sequencing count data. While the estimates generated 327 by these models may be less intuitive and temporally reversed from our hypothesis, the 328 associations between placental microRNA expression with growth trajectory patterns ascertained 329 using such models are nevertheless valid. We adjusted for likely confounders in our study, but 330 cannot rule out the possibility that unmeasured or residual confounding remains in our analysis. 331 Potential confounders are maternal BMI, gestational weight gain and pre-eclampsia. Sensitivity 332 analysis suggests that the addition of these variables to our models would have little to no effect 333 on our conclusions, though some of our top findings are sensitive to one or both. To limit 334 unknown statistical confounding (e.g. differing cellular composition between individuals or 335 population stratification), our models are adjusted using surrogate variable analysis. This is a 336 data-driven approach that may incorrectly estimate confounding elements, adding variability to 337 our model. Recent research indicates that SVA is one of the more robust and reliable methods in 338 studies such as ours 61 . Lastly, the cohort utilized in this study consisted predominantly of 339 healthy white mothers from a rural New England region of the United States, potentially limiting 340 the generalizability of our findings. 341 342 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. missingness rate (y-axis) was low, meaning that observations at each measurement (obsage) 510 remained high despite missing data after 12 months (12m, Nage). 511 512 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)  The growth curves representing the minimum and maximum growth intensities are plotted in 518 cyan and purple, respectively. Deviations in growth intensity appear as rotations 519 (counter/clockwise) of the mean curve. 520 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)  proportional to -log10(P-value), such that larger points represent smaller p-values. 526 527 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint representing the involved placental cell types from which microRNAs in our study could arise: 531 fetal endothelial cells (FEC), cytotrophoblasts (CTB), syncytiotrophoblasts (STB) and 532 extravillous trophoblasts (EVT).Growth trajectory microRNAs are listed (in black) along with 533 their influential targets (gray) in the processes we predict they may influence (trophoblast stem 534 maintenance, proliferation and differentiation to the syncytiotrophoblast (STB) or extravillous 535 trophoblast (EVT) terminal lineage). microRNAs and their putative targets in the central blue or 536 outer red areas are predicted to encourage or inhibit the given process, respectively. Created with 537 BioRender.com. 538 539 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted February 4, 2022. ; https://doi.org/10.1101/2022.02.03.22270310 doi: medRxiv preprint