The study utilized samples from Estonian REPROMETA (full study name “REPROgrammed fetal and/or maternal METAbolism”; recruitment 2006–2011) and HAPPY PREGNANCY (full: “Development of novel non-invasive biomarkers for fertility and healthy pregnancy”; 2013–2015) data sets. Both studies were approved by the Ethics Review Committee of Human Research of the University of Tartu, Estonia (Permissions No 146/18, 27.02.2006; 150/33, 18.06.2006; 158/80, 26.03.2007; 221/T-6, 17.12.2012; 286/M-18, 15.10.2018). All study participants were recruited, and the study material was collected at the Women’s Clinic of Tartu University Hospital, Estonia. Written informed consent to participate in the study was obtained from each individual prior to recruitment. The study was carried out in compliance with the Helsinki Declaration and all methods were carried out in accordance with approved guidelines. All participants were of white European ancestry and living in Estonia.

Genome-wide profiling of placental miRNome was performed for 52 placental samples representing first (n = 5) or second (n = 7) trimester of gestation, and term pregnancy (n = 40) (Table 1). Term placental samples were drawn before or shortly after delivery during the REPROMETA study (Supplementary Methods). The analyzed 40 term pregnancy cases (delivery after 37th gestational week, g.week) represented normal gestation, preeclampsia (PE), gestational diabetes (GD), and delivery of a small- or large-for-gestational-age newborn (SGA and LGA, newborns <10th or >90th birth weight centile, respectively). Each clinical subgroup included eight cases, matched for the gestational age, delivery mode and proportions of male/female newborns. First and second trimester placental samples had been collected from women who underwent elective surgical termination of pregnancy or medically induced abortion due to maternal medical risks. In all analyzed cases, fetal anomalies and gross chromosomal aberrations were excluded. The definition of clinical subgroups, details of collection and processing of placental samples, and DNA and RNA extraction protocols are provided in Supplementary Methods and recent publications (Sõber et al., 2015; Reiman et al., 2017; Kikas et al., 2019, 2020).

Initial small-RNA libraries were prepared from 1 μg total RNA (TruSeq Small RNA kit, Illumina), followed by miRNA enrichment (Caliper LabChipXT, PerkinElmer) according to manufacturer’s protocols. Small RNA-Seq libraries were sequenced on Illumina HiSeq 2000. Library preparation and sequencing were conducted in FIMM Sequencing Laboratory, University of Helsinki, Finland. Quality control of the raw reads was performed using FastQC (ver. 0.11.7) and MultiQC (ver. 1.7) (Ewels et al., 2016). Trimmomatic (ver. 0.38) was implemented to remove adapters and trim the quality of reads with the following settings - ILLUMINACLIP:2:30:9, LEADING:3, CROP:50, TRAILING:3, SLIDINGWINDOW:4:20, MINLEN:16. Reads were aligned to human genome reference GRCh38 using bowtie (ver. 1.2.2, settings: -n 1 -l 20 -q -m 40 -k 1 -t –best –strata) (Langmead et al., 2009). miRNA quantification was performed using featureCounts from the Rsubread package (ver. 1.20.6) (Liao et al., 2019) for R with miRNA annotations from miRBase 22.1 as reference (Kozomara et al., 2019).

From 2,652 quantified placental miRNAs only those with median raw read counts over 50 across all analyzed samples (n = 417; 15.7%) were included in statistical analyses. All subsequent computational profiling and differential expression analyses were implemented using read counts normalized with DESeq2 (ver. 1.22.2) package for R with default settings (Love et al., 2014; Supplementary Table 1). In addition, counts per million reads mapped (CPM) were quantified for the graphical presentation of miRNA expression in subgroups. False-discovery rate (FDR) in differential expression tests was applied according to Benjamini and Hochberg (1995). Test results with FDR P < 0.05 were considered as statistically significant. Placental miRNome was compared for the following subgroups: first (n = 5) vs. second trimester (n = 7), second trimester vs. normal term (n = 8) pregnancy; late gestational complications (PE, GD, SGA, LGA) vs. normal term pregnancy (each group n = 8); term pregnancy 46, XX (n = 21) vs. 46, XY (n = 19). Comparisons of miRNA expression between the three trimesters of gestation were carried out without and with adjustment for sex as cofactor. miRTarBase database (Huang et al., 2020) was used to assemble the lists of experimentally validated target genes for differentially expressed miRNAs (DEmiRs) across normal pregnancy, and for miRNAs encoded by the C19MC and C14MC clusters. Only target genes with high confidence were considered for the gene enrichment analysis (see below). Expression correlations of miR-143-3p, miR-92a-3p, miR-26a-5p, as well as C19MC and C14MC microRNAs with high confidence target genes were calculated using Kendall correlation coefficient (parameter tau).

mRNA/lincRNA expression data was derived from the RNA-Seq datasets that had been previously generated for the same placental samples as utilized for miR-Seq. Placental RNA-seq library preparation, sequencing and raw data processing are detailed in Supplementary Methods and in previous studies (Sõber et al., 2015, 2016; Reiman et al., 2017; Pilvar et al., 2019).

Analysis of inter-relatedness between the expression of miRNAs and mRNA/lincRNA genes in 40 term placentas also utilized the above-mentioned published RNA-Seq data. Expressional correlation of miRNA/mRNA transcripts was evaluated using Spearman’s correlation coefficient (parameter rho). Correlation analysis included 66 miRNAs showing differential expression in PE in the miR-Seq dataset and 16,567 genes with raw median read counts >50 in the RNA-Seq dataset. Spearman’s rho values for 1,093,422 miRNA-gene pairs were estimated in R and visualized as a heatmap, using R package heatmap.2 (Gregory et al., 2015). Lists of genes showing confident expressional correlation with miRNA hierarchical cluster groups G1-G5 were formed using the following criteria: median Spearman’s rho across 40 term placentas <−0.3 and for individual samples <−0.1 (negatively correlated genes); or median rho > 0.3 and for individual samples higher than rho > 0.1 (positively correlated genes). These gene lists were used as input for the gene enrichment analysis for in silico functional profiling.

All gene enrichment analyses based on miRNA hierarchical cluster groups were implemented in g:Profiler (ver. 1760) with default settings (Reimand et al., 2016). The recommendation for a more conservative analysis to compute functional enrichment in a custom gene list instead of all human genes in Ensembl database was followed to avoid overestimating statistically significant results (Reimand et al., 2016; Raudvere et al., 2019).

To avoid potential confounding effect of gestational expression dynamics, the discovery analysis of placental miR-eQTLs included only term placental samples (n = 40). SNV genotypes were derived from the previously published genome-wide genotyping dataset of the same placental samples [Illumina HumanOmniExpress-12-v1 BeadChip (>733,000 SNVs; median spacing 2.1 kb)] (Kasak et al., 2015; Pilvar et al., 2019). The analysis was targeted to ±100 kb window extending to both directions from the start and end of miRNA genes, annotated based on miRBase (ver. 22.1). The genomic regions flanking the analyzed 417 miRNAs included 6,274 common SNVs (MAF > 0.1). In total, 17,302 linear regression association tests were carried out between SNV genotypes and miRNA expression levels, quantified as normalized miRNA read counts. All tests with miR-eQTLs were implemented in PLINK v1.07 using fetal sex and gestational age as cofactors (Purcell et al., 2007). The results were corrected for multiple testing using the Benjamini–Hochberg method, with cut-off FDR < 0.05. All of the miR-eQTLs were tested for Hardy–Weinberg equilibrium (Supplementary Table 2).

The REPROMETA study recruited 366 pregnant couples before or shortly after delivery of a singleton newborn at term (Supplementary Table 3). The cases represented pregnancies with uncomplicated gestation, PE, GD, SGA or LGA. Maternal and newborn clinical and epidemiologic data were documented retrospectively from self-reported questionnaires and medical records, collected biological materials included placental tissue (available for 326 cases) and parental blood samples.

The HAPPY PREGNANCY study recruited prospectively 2,334 pregnant women during their first antenatal visit. Longitudinal clinical and epidemiological data covers reproductive history, parental lifestyle, the course and outcome of pregnancy. The collected biological material included placental tissue (available for 1,772 cases), maternal blood and urine samples.

In both studies, the diagnosis of PE and GD followed the international guidelines at the time of recruitment (Metzger, 2010; American College of Obstetricians and Task Force on Hypertension in Pregnancy, 2013). Fetal growth was evaluated using the gestational age and sex adjusted weight centiles based on the Estonian Medical Birth Registry data (Sildver et al., 2015). Details of the REPROMETA and HAPPY PREGNANCY studies are provided in Supplementary Methods and in recent publications (Kikas et al., 2019, 2020).

Genetic association testing between the identified miR-eQTLs and pregnancy traits (placental weight; newborn’s weight and height, head and chest circumference; diagnosis of PE or GD) was carried out using either linear or logistic regression (additive model) adjusted for fetal sex and gestational age as cofactors. Results were corrected for multiple testing using the Benjamini-Hochberg method, with cut-off FDR < 0.05. Initially, all miR-eQTLs were tested in the discovery dataset (n = 40, Table 1). Three SNVs (rs12985296, rs7046565, and rs12420868) were further analyzed in the REPROMETA (n = 326) and HAPPY PREGNANCY (n = 1,772) samples, including all cases with available placental tissue for genotyping (Supplementary Table 3).

Genotypes were generated with the TaqMan SNP Genotyping Assays (Supplementary Methods and Supplementary Table 4) using recommended experimental conditions (Applied Biosystems, Life Technologies). All association tests were carried out in PLINK v1.07 using (Purcell et al., 2007). Meta-analysis of REPROMETA and HAPPY PREGNANCY datasets was implemented in R package meta (ver. 4.15-1) (Balduzzi et al., 2019), under fixed-effect model.

The dataset of 52 placental miRNomes was generated for the samples collected from first and second trimester and term pregnancy cases (n = 5, 7 and 40, respectively; Table 1). All subsequent analyses included 417 of total 2,652 identified miRNAs (15.7%, Table 2A), filtered for transcript levels that allow confident statistical testing (median raw read counts >50 across analyzed samples; Supplementary Table 1).

Placental miRNome in uncomplicated pregnancies was assessed in 20 samples: five cases representing first trimester [median 60 (51–81) gestational days, g.days], seven second trimester [121 (108–140) g.days], and eight term pregnancy [284 (260–291) g.days] cases. Overall, a broad variability in expression levels of individual placental miRNAs were measured in all trimesters (Figure 1A and Supplementary Figure 1). Median expression level of 417 analyzed miRNAs did not differ across the three trimesters of pregnancy (Kruskal–Wallis test, p = 0.24). However, a non-significant decreasing trend of median values was observed from 140 CPM (range 1–79,604) in first trimester, to 132 CPM (4–123,631) in second trimester and 103 CPM (7–172,159) in term placental samples.

A handful of miRNAs were identified with extremely high (CPM > 25,000) transcript levels throughout gestation (Figures 1A,B). Among these, miR-143-3p showing gradually increasing transcript levels from early pregnancy until term appeared as a potential marker for pregnancy progression and placental maturation. In total 47 high-confidence target genes for miR-143-3p were identified in the miRTarBase database (Supplementary Table 5). Majority, 93% of them clustered to the Gene Ontology (GO) pathway ‘cellular response to stimulus’ (GO:0051716; FDR = 7.2 × 10^–10; Supplementary Table 6). More specific enriched (FDR < 0.05) functional categories relevant to pregnancy included, e.g., ‘PI3K-Akt signaling pathway’ (KEGG:04151; 34% of target genes), ‘Endocrine resistance’ (KEGG:01522; 31%), ‘collagen metabolic process’ (GO:0032963; 12%), ‘regulation of glucose transmembrane transport’ (GO:0010827; 12%), ‘growth factor binding’ (GO:0019838; 12%). Other examples of gestational age-dependent highly expressed miRNAs were miR-92a-3p, miR-30d-5p and miR-512-3p (specifically increased in first trimester), miR-26a-5p (second trimester), miR-125a-5p (at term). Some major miRNAs exhibited constant expression throughout gestation (e.g., miR-181a-5p, miR-22-3p).

Comparative assessment of the transcript levels of miR-143-3p, miR-92a-3p and miR-26a-5p with the expression of predicted target genes in the corresponding RNA-Seq dataset revealed a substantial group of mRNA/lincRNAs with reverse expression compared to the miRNA (Figure 1C and Supplementary Table 7). However, significant expressional correlations identified for miR-143-3p and miR-92a-3p (Kendall rank correlation coefficient, p < 0.05) included nearly equal proportions of reversely and directly correlated target genes (9/7 and 8/5, respectively; binominal test, p > 0.58). There was an under-representation of genes showing tight negative vs. positive correlation with the miR-26a-5p transcript levels (8/21; p = 0.024). Interestingly, positively correlated loci represented cancer driver genes (MYC, PRKCD, RB1, CDK6), molecules involved in immune response (MALT1, PIK3CG), as well as cellular and hormonal signaling (HGF, PTPN13, IGF1, SMAD1, ESR1, CTGF). These functional groups are well-known to be involved in placental development and role in supporting the pregnancy.

The majority, 319 of 417 (76.5%) of tested miRNAs, exhibited significant gestational expression dynamics (Table 2B and Supplementary Tables 8, 9). In total, 227 (54.4%) miRNAs were differentially expressed between first and second trimester [FDR < 0.05; log₂(FC) from −4.91 to 2.84; 125 down- and 102 upregulated], and 211 miRNAs (50.1%) between second trimester and term pregnancy placental samples [FDR < 0.05; log₂(FC) from −2.41 to 2.52; 110 down- and 101 upregulated]. More than a quarter of tested miRNAs (n = 119/417; 28.5%) represented DEmiRs in both comparisons, indicating their potential critical contribution in fine-tuning placental transcriptome profile in gestational age-dependent manner until term (Figure 2A and Supplementary Table 10).

The 319 placental miRNAs exhibiting trimester-dependent differential expression were explored for their predicted target genes in the miRTarBase database (Supplementary Table 11). An overrepresentation of target genes for upregulated compared to downregulated miRNAs was observed – 1305 vs. 773 genes in first vs. second trimester and 1207 vs. 888 in second trimester vs. term comparisons (χ²-test, p = 6.2 × 10^–4). Only a small fraction of genes (14.1–17.7%; Supplementary Table 11) represented trimester-specific targets of dynamic miRNAs, despite nearly two thirds of DEmiRs (200 of 319) being detected significant increase or decrease only in one trimester (Figure 2A). The shared list of target genes of dynamic miRNAs supported the coordinated action of the placental miRNome in modulating the expression of key placental genes during gestation.

In silico functional profiling of genes targeted by placental miRNAs with progressively increasing transcript levels from early pregnancy to term revealed a significant enrichment of GO pathways representing broad basic cellular and tissue functions (FDR < 0.05; Table 3 and Supplementary Tables 12, 13). Examples of modulated processes were ‘protein binding’ (GO:0005515; >90% of target genes) and ‘cytosol’ (GO:0005829; ∼48%), ‘cytoplasmic stress granule (GO:0010494; ∼31% of the pathway genes) and ‘extracellular matrix structural constituent conferring tensile strength’ GO:0030020; 43.5%). Target genes for downregulated miRNAs in second compared to first trimester placentas represent more focused functional groups, possibly supporting the placental role in fine-tuning the fetal development in mid-gestation. These include, e.g., ‘membrane-enclosed lumen’ (GO:0031974; 54.4% of target genes) and ‘cellular response to chemical stimulus’ (GO:0070887; 52.9%); as well as ‘thymocyte apoptotic process’ (GO:0070242; 41.2% of pathway genes), ‘DNA alkylation’ (GO:0006305; 24.4%), ‘PML body’ (GO:0016605; 22.6%), ‘regulation of cell size’ (GO:0008361; 16.5%) and ‘regulation of calcium ion transport’ (GO:0051924; 15.6%). Loci that represent targets for downregulated miRNAs at term cluster in biological pathways relevant to the preparation for the delivery, such as ‘positive regulation of smooth muscle contraction’ (GO:0045987; 42.1% of pathway genes), ‘phosphatidylinositol-3-phosphate biosynthetic process’ (GO:0036092; 39.4%) and ‘regulatory RNA binding’ (GO:0061980; 38.7%).

The tested 417 placental miRNAs were assigned to one of nine subgroups representing their temporal expression dynamics pattern across three trimesters of pregnancy (Table 2B and Supplementary Figure 2). The most prevalent expression dynamics pattern represented miRNAs exhibiting specifically high transcript levels in early pregnancy (n = 67 miRNAs, ∼16%). The second frequent pattern reflected miRNAs that were downregulated only at term (n = 54, ∼13%). High miRNA expression restricted to second trimester was the rarest observed expressional pattern (n = 26, ∼6%). A stable expressional window from early pregnancy to term was identified for 98 miRNAs (23.5%) with 900 predicted target genes (Supplementary Tables 14, 15). Among these, the most significantly enriched biological pathways (FDR < 1.0 × 10^–17) are also implicated in basic cellular and tissue functions, ‘protein binding’ (GO:0005515; 89.5 of target genes), ‘cellular response to organic substance’ (GO:0071310; 43.8%) and ‘nuclear lumen’ (GO:0031981; 43.5%) (Supplementary Table 16).

When incorporating sex as a cofactor (46, XX; 46, XY) in differential expression testing between the three gestational trimesters, the list of DEmiRs remained largely the same (Figure 1D and Supplementary Tables 8, 9). In second trimester [(M)ale/(F)emale, n = 4/3) compared to normal term placentas (M/F, n = 5/3), only six out of 211 (2.8%) DEmiRs were ‘lost’ after correction for fetal sex. The analysis of first (M/F, n = 2/3) compared to second trimester pregnancy samples with or without fetal sex as a cofactor resulted in 200 shared miRNAs, whereas four and 27 miRNAs were additionally detected in one sub-analysis. Only a single miRNA, miR-193-3p, was highlighted in both tests comparing trimester-specific miRNomes and may represent a locus with sex-specific expression level. The nature of other minor differences in the identified miRNAs lists cannot be speculated as testing small sample sets may result in spurious outcomes.

To evaluate further the role of fetal sex on placental miRNome, term pregnancy samples representing 46, XY (n = 19, 260–291 g.days) vs. 46, XX (n = 21, 262–288 g.days) cases were compared (Supplementary Table 17). Statistically significant differential expression was identified for three miRNAs: X-linked miR-361-5p, autosomal miR-378a-3p and miR-130b-3p. Additionally, miR-101-5p showed placental sex-modulated expression when the analysis was adjusted for pregnancy complications. These miRNAs did not overlap with the candidates identified in trimester-specific analysis. Overall, the analysis outcome suggested that fetal sex is not a major modulator of placental miRNome.

A notably high fraction, 125 of 417 (∼30%) expressed miRNAs belonged to the primate-specific miRNA cluster C19MC (detected mature placental miRNAs, n = 65; 15.6%) or to the eutherian-specific clusters C14MC (n = 58; 13.9%) and miR-371–373 (n = 2, 0.5%) (Table 2A and Supplementary Table 1). These clusters showed markedly different patterns of gestational expression dynamics. About ∼2/3 of C19MC and the miR-371–373 clusters are specifically highly transcribed in early pregnancy with a significant drop in second trimester and a slight increase at term (Figure 3A, Table 2B, and Supplementary Tables 8, 9). The C14MC cluster showed diverse expression in first trimester, but more coordinated transcript levels in later gestational ages. The majority of C14MC miRNAs showed high expression in second trimester and significant downregulation before term. Only five C14MC, but 22 C19MC miRNAs exhibited stable expression levels from early pregnancy until delivery.

Among confident target genes in miRTarBase database reported for the C19MC and C14MC miRNAs, 63/76 and 215/262 loci were expressed in the current RNA-Seq dataset (Supplementary Tables 18, 19). Indicating the distinct roles of C19MC and C14MC in shaping human placental transcriptome across gestation, only 21 of these target genes overlap. Most target genes are modulated explicitly by either C19MC or C14MC miRNAs. Assessment of gestational expression dynamics of C19MC and C14MC miRNAs and their predicted target genes supports their mutual functional relationship (Figures 3A,B). The major cluster of C19MC miRNAs have the lowest transcript levels at ∼180–200 g.days and a substantial proportion of their predicted target genes show increased expression during mid-gestation. Interestingly, the minor cluster of C19MC miRNAs is strongly positively correlated with over 1/3 of the analyzed genes. The majority of C14MC target genes belong to distinct co-expression groups that form clusters based on either negative or positive transcriptional correlation with miRNAs. Large notable clusters of miRNA–gene pairs showing directly and inversely correlated expression dynamics suggested potential functional relationships in both scenarios.

The top GO terms with significant enrichment of C19MC target genes (FDR < 1.0 × 10^–5, Supplementary Table 20) were indicative of processes related to transcriptional activation and cellular signaling that are critical in early pregnancy, such as ‘response to organic substance’ (GO:0010033), ‘transcription regulatory region DNA binding’ (GO:0044212), ‘kinase binding’ (GO:0019900). Notably, most significantly enriched KEGG and REAC categories comprise of genes implicated in various cancers or are responsible for signaling pathways of cancer driver genes. The C14MC miRNAs were implicated in more basic cellular functions that are essential throughout pregnancy, such as ‘cellular response to chemical stimulus’ (GO:0070887), ‘protein binding’ (GO:0005515), ‘intracellular organelle lumen’ (GO:0070013). However, some target genes of C14MC reflect specific placental features, such as ‘DNA-methyltransferase activity’ (GO:0009008) potentially linked to epigenetic programming or ‘regulation of cell size’ (GO:0008361) possibly referring to the development of large multinucleated placental cell types.

Placental miRNomes representing term cases of late-onset preeclampsia (LO-PE), gestational diabetes (GD), small- and large-for-gestational-age newborns (SGA, LGA) were tested for differential expression in comparison to uncomplicated pregnancies (n = 8 in each group; all cases after 37th g.week; Table 1). Only PE placentas demonstrated a major shift in their miRNome profile that affected 66 of 417 (15.8%) miRNAs (FDR < 0.05; Figure 2 and Supplementary Tables 21, 22). None of these miRNAs was detected to be modulated by fetal sex (Supplementary Table 17). Seven significantly upregulated miRNAs overlapped with the placental DEmiRs reported in early-onset PE cases (EO-PE, before 34th g.week) (Awamleh et al., 2019; Figure 2B). Several of these showed large changes in their expression level: miRNAs miR-210-3p (FC = 2.64), miR-193b-3p (2.53), miR-193b-5p (2.29), miR-365b-3p (1.93), miR-365a-3p (1.92), miR-520a-3p (1.82) (Figure 4 and Supplementary Table 21). Notably, 13 of 38 (34%) upregulated DEmiRs were transcribed from the C19MC cluster and only one from C14MC, whereas 11 of 28 (39%) downregulated miRNAs were transcribed from the C14MC and none from the C19MC cluster (Figure 2C).

Differentially expressed miRNome in PE was comprised of both dynamic and stable miRNAs. No specific pattern of normal gestational dynamics was preferentially altered (Figures 2A,C). Several miRNAs normally downregulated at term were characterized by increased transcript levels in PE placentas corresponding to their typical mid-gestation values (e.g., miR-210-3p, miR-31-5p, miR-96-5p, miR-193a-3p, miR-519a/b-5p; Figure 4). On other occasions, miRNA expression in PE placentas was significantly downregulated compared to other analyzed samples (e.g., miR-196b-5p, miR-411-3p, miR-1247-5p). PE miRNome also showed aberrant upregulation of several miRNAs that are normally stably expressed across gestation (e.g., miR-365a/b-3p).

Six miRNAs with altered transcript levels in LO-PE represented also DEmiRs in term SGA placentas – upregulated miR-210-3p, miR-512-5p, miR-32-5p, miR-19a-3p and miR-590-3p, and downregulated miR-379-5p (Supplementary Table 23). No miRNAs exhibited significant expressional changes in the GD and LGA cases compared to normal term placenta (Supplementary Tables 24,25).

Correlation analysis between the expression levels of 66 placental DEmiRs identified in PE and placental transcriptome was performed using the corresponding miR-Seq and RNA-Seq datasets of 40 term pregnancy samples. Hierarchical clustering based on the expressional correlation with the transcript levels of 16,567 genes assigned the tested miRNAs into five groups G1-G5, each containing 6–22 miRNAs (Figure 5 and Supplementary Data 1). In these groups, there was a highly non-random distribution of miRNAs from C19MC (G1:10 miRNAs, G5:3) and C14MC (G4:9, G3:2, G5:1) clusters (χ²-test, p = 1.5 × 10^–5), providing further support to their distinct roles in modulating placental transcriptome. The C14MC cluster outlier miRNA that didn’t properly fit in either groups G4 or G5 was miR-376a-5p. Furthermore, this miRNA already stood out in the differential expression analysis with an opposite behavior compared to the rest of the C14MC miRNAs (Figure 2C).

Six of 10 miRNAs with the highest number of negatively correlated placental transcripts belonged to miRNA group G2, e.g., miR-490-3p (median Spearman’s rho across 40 placentas <–0.3, 4375 genes vs. rho > 0.3, 55 genes), miR-27a-3p (2957 vs. 25), miR-362-5p (2823 vs. 427) (Supplementary Tables 22, 26). Among the C19MC miRNAs that are significantly upregulated in PE placentas, the largest proportion of negatively correlated genes was detected for miR-522-5p (1780 negative vs. 148 positive correlations) and miR-518a-5p (1731 vs. 154) clustering to group G1. In general, groups G1 and G2 showed confident negative correlation with a large number of genes (median rho <–0.3 and for all individual samples <–0.1: 342 and 239 genes, respectively; Supplementary Table 27). In contrast, groups G3-G5 were positively correlated with a substantial proportion of the placental transcriptome (median rho > 0.3 and for all individual samples >0.1; 1077, 484, and 2537, respectively). This included also C14MC cluster miRNAs that are significantly downregulated in PE.

Gene enrichment analysis was performed to clarify the potential functional link between the PE-linked miRNAs clustered into groups G1–G5 and their most significantly correlated placental transcripts (Figure 5 and Supplementary Tables 28–32). Group G2 miRNAs appeared as potential key modulators of genes implicated in basic cellular processes, such as organelle function and ribosomal biogenesis, DNA and RNA metabolism, cell cycle and senescence (Table 4). However, expression of miRNAs in groups G1, G3, and G4 was correlated with functionally more specific gene categories, encoding proteins involved in extracellular matrix formation (G1, G4), DNA replication and immunoglobulin complex (G3), synaptic and ion channel activity, adrenergic signaling, renin and insulin secretion (G4).

Positive expressional correlation between miRNAs in group G5 and thousands of genes may reflect either true functional relationship or alternatively, co-correlated expression of miRNAs and some housekeeping genes through shared upstream regulatory modulators. Notably, ∼50% of them belonged to the functional category ‘nucleus’ (GO:0005634; enrichment P = 5.0 × 10^–3), being responsible for various DNA, chromatin and RNA-related processes. Another large category represented ‘immune system process’ (GO:0002376, 19.5% of query genes, P = 1.8 × 10^–4).

Finally, the effect of SNVs on the placental miRNome was investigated. Placental Expression Quantitative Trait Loci (eQTLs) for 417 miRNAs were mapped by genetic association testing between their transcript levels in 40 term placental samples and genotypes of 6,274 common SNVs located ±100 kb from the miRNA genes. In total 66 miR-eQTLs for 16 miRNAs were detected (FDR < 0.05; 3.8% of tested miRNAs; Figure 6, Supplementary Table 33, and Supplementary Data 2). Four of 16 placental miRNAs modulated by eQTLs had also been identified as DEmiRs in PE (miR-30a-5p, miR-210-3p, miR-490-3p, miR-518-5p). Despite the limited sample size, the effect of miR-eQTL on some miRNAs was observed in all three trimesters of pregnancy (e.g., pairs rs447001 and miR-130b-3p/5p, rs2427554 and miR-941, rs12642661 and miR-1269a). The most extreme identified SNV-miRNA pair was rs7046565 (A/G) and miR-3927-3p. The major allele AA-homozygosity completely suppressed the expression of miR-3927-3p. This effect was also detected in second trimester placental samples (Figure 6).

Among 66 identified placental miR-eQTLs, 18 eQTLs were unique to placental miRNAs and 48 have also been reported in the GTEx database as expressional modulators of 53 coding genes (The GTEx Consortium, 2020). In our placental RNA-Seq dataset, 32 of them were expressed. Statistically significant associations (FDR < 0.05) were detected with transcript levels of the KLHL3 (rs10515496), SNX11 (rs11651097), ANO9 and PTDSS2 (rs12420868) genes (Supplementary Table 34). However, each of these statistical associations were weaker than the originally detected effect on the adjacent miR-874-3p, miR-152-3p, and miR-210-3p, respectively.

Three miR-eQTLs (SNVs: rs12420868, rs12985296, rs7046565) showing some nominal associations in the discovery dataset (n = 40, Supplementary Table 35) were targeted to replication testing with pregnancy traits in the REPROMETA (n = 326) and HAPPY PREGNANCY (n = 1,772) pregnancy-related cohorts (Supplementary Table 3 and Supplementary Methods). No statistically significant associations were identified with the height, weight, head and chest circumference of newborns, placental weight and PE or GD incidence in independent cohorts or their meta-analyses (all tests, FDR > 0.05; Supplementary Table 36). A non-significant trend was detected between rs12420868 (eQTL for miR-210-3p) and newborns’ head circumference (meta-analysis: nominal P < 0.05; Supplementary Figure 3).