Integrated Systems Biology Approach Identifies Novel Maternal and Placental Pathways of Preeclampsia

1Preeclampsia, one of the most severe obstetrical complications affecting 5–8% of pregnant women, is a leading cause of maternal and perinatal morbidity and mortality. In addition, pathologic changes in the affected mothers and fetuses lead to a higher risk of subsequent metabolic and cardiovascular diseases later in life, further increasing healthcare costs. In spite of the severity of the problem, there is yet no early diagnosis of all forms of preeclampsia, and the current therapy is still based on the delivery of the placenta, given the complexity of the disease and the lack of insight into the early perturbed molecular pathways.

2Indeed, preeclampsia is a syndrome with heterogeneous etiology and a spectrum of phenotypes (Figure 1). It may affect women at varying gestational ages with different degrees of severity and consequences for the fetus. The current classifications of preeclampsia are based upon its severity and the timing of clinical presentation, mostly dividing preeclampsia into preterm (<37 weeks) or term (≥37 weeks) and early-onset (<34 weeks) or late-onset (≥34 weeks) phenotypes. Preterm preeclampsia has a more severe clinical presentation and is often accompanied by fetal growth restriction compared to term preeclampsia. However, severe maternal disease and fetal growth restriction may be observed in both term and preterm preeclampsia, and their presentation may not be associated with each other, suggesting that the clinical phenotype is the result of an interplay between various factors and disease pathways, also supported by observations in animal models.

3A growing body of evidence offers support for the conclusion that the maldevelopment and/or dysfunction of distinct trophoblast lineages of the placenta have a central role in the pathogenesis of preeclampsia, and that the severity of the placental disease is subsequently reflected in the clinical phenotype of this syndrome. In the preclinical stage, extravillous trophoblast (EVT) development may be impaired, leading to EVT dysfunction, shallow trophoblast invasion, failure of the physiological transformation of the maternal spiral arteries, abnormal blood-flow to the placenta, and histological changes consistent with maternal vascular malperfusion. The frequency and severity of these lesions decrease from the preterm toward the term phenotype of preeclampsia (Figure 1), mirrored by the decreasing prevalence of fetal growth restriction or the delivery of small-for-gestational age (SGA) neonates. Thus, EVT development is less frequently and extensively affected in late-onset preeclampsia, which constitutes about 90% of all cases.

4In preterm preeclampsia, failure of the remodeling of the spiral arteries may lead to abnormal blood flow to the placenta and subsequently to placental structural damage and ischemic stress, villous trophoblast (VT) dysfunction and the release of detrimental placental substances (e.g. anti-angiogenic factors, pro-inflammatory cytokines, and syncytiotrophoblast debris) into the maternal circulation. As a consequence, the terminal pathway of preeclampsia, featuring an anti-angiogenic state and exaggerated maternal systemic inflammation, occurs in most cases, and its intensity correlates with the severity of preeclampsia, which may be coupled with damage to the maternal endothelium and to the kidneys, liver, and central nervous system during the clinical phase. VT development can also be impaired in preeclampsia, especially in the preterm form associated with SGA, where VT turnover is affected together with morphometric features.

5Term preeclampsia is characterized by a lesser magnitude of maternal systemic inflammatory and anti-angiogenic states. This phenotype may result from different stressors other than maternal vascular malperfusion and ischemia of the placenta, which include various preexisting maternal disorders, such as obesity, chronic hypertension, diabetes, and metabolic, kidney, and autoimmune diseases. These stressors may still trigger placental stress and VT dysfunction and induce a maternal pro-inflammatory milieu. Alternatively, maternal endothelial dysfunction may result from an exaggerated sensitivity to factors released from the placenta, which increases the risk of preeclampsia upon maternal genetic predisposition for cardiovascular disease. Indeed, preeclampsia has a genetic predisposition with high heritability of both phenotypes, and it shares common risk alleles with coronary artery disease.

6In spite of extensive research efforts, our understanding of the early pathologic pathways of preeclampsia has been limited given several obstacles. First, the complexity of the disease pathways and the heterogeneity of the syndrome have not been investigated in an integrative manner in both maternal and placental compartments throughout pregnancy. Second, it has been impossible to investigate the early placental disease pathways because of the invasive nature of placental biopsy and the limited information on placental functions obtained non-invasively. Consequently, an increasing number of high-dimensional studies aiming to detect molecular signatures of preeclampsia either in the placenta or in maternal blood have mostly targeted later stages of pregnancy, at a more advanced stage of placental development and pathology. Third, animal models of preeclampsia fail to mimic early placental pathways of preeclampsia due to the anatomical and physiological uniqueness of deep placentation in humans. Fourth, in vitro studies on human placental development and trophoblast functions are hindered by the lack of self-replicating trophoblast stem cells with the ability to differentiate into both VTs and EVTs.

7Here, we used a systems biology approach that integrated various omics and targeted methods to investigate both placental and maternal compartments and most aspects of preeclampsia, including placental disease; fetal and maternal outcomes; environmental, maternal, microenvironmental, and epigenetic factors; and trophoblastic functions. In the first placental study, we performed extensive investigations of the placenta at histologic, transcriptomic, epigenetic, and protein levels to target molecular pathways at the center of the disease. Molecular changes were correlated with maternal and neonatal morbidities associated with preeclampsia to uncover placental pathways affecting either maternal or fetal wellbeing. In the second, maternal study, we employed maternal blood proteomics and “virtual” liquid biopsy of the placenta to reveal blood factors of maternal or placental origin that can reflect disease conditions in early pregnancy. In the third, trophoblast study, we utilized in vitro functional assays on the trophoblast to investigate hub transcription factors at the center of placental disease gene modules and to model their in vivo involvement in placental pathways associated with maternal or placental/fetal disease (Figure 2).

8Given that the pathogenesis of preeclampsia has been implicated to originate from the placenta, we first aimed to investigate placental transcriptomic changes leading to placental dysfunction as well as regulatory networks involved in the pathologic pathways. The combined analysis of preterm preeclampsia cases and gestational age-matched controls (n = 17) in a Hungarian patient population (Table 1) in our placental microarray data revealed 1,409 differentially expressed (DE) genes (Data S1 in Supplementary Material), which are involved in fundamental cellular processes, including blood pressure (BP) regulation, apoptosis, development, hormone secretion, metabolism, homeostasis, and signaling (Figure 3A). DE genes included 137 transcription regulatory genes and 38 predominantly placenta-expressed genes. This latter set of genes (n = 164, Data S2 in Supplementary Material), which was defined by BioGPS microarray data (n = 153) or expression data from separate studies in the lack of BioGPS data (n = 11), was enriched among DE genes in preeclampsia [odds ratio (OR) = 3.4, p = 6.9 × 10−9]. This suggests that genes predominantly expressed by the placenta have pathologic and diagnostic significance in preeclampsia, a phenomenon indicated earlier by our targeted studies.

9Subsequently, we investigated the genomic links among DE genes by searching for genomic regions associated with the observed placental transcriptomic changes. We found that Chr6 (OR = 1.54, q = 1.6 × 10−3) and Chr7 (OR = 1.42, q = 0.02) were particularly affected by these gene expression changes (Data S3 in Supplementary Material), while Chr19 (OR = 2.6, q = 0.02) was enriched in dysregulated transcription regulatory genes (Data S4 in Supplementary Material). Of interest, predominantly placenta-expressed genes were also enriched on Chr19 (OR = 2.5, q = 0.002) (Data S5 in Supplementary Material). Figure 3B shows the non-random genomic localization of DE genes in preterm preeclampsia and the pronounced gene dysregulation associated with Chr19. These results are consistent with the fact that Chr19 harbors large transcription regulatory gene families, and probably reflect its regulatory role in placental/trophoblastic gene expression and their dysregulation in preeclampsia.

10Next, we aimed to identify functional links among DE genes by identifying gene co-expression network modules and hub transcription regulatory genes driving differential expression in the modules. Weighted co-expression network analysis (WGCNA) was conducted among DE genes resulting in the assignment of these into four major modules, labeled as M1 (green, n = 506), M2 (red, n = 442), M3 (blue, n = 381), and M4 (orange, n = 74) (Figure 3C; Data S1 in Supplementary Material). Most predominantly placenta-expressed genes belonged to modules M1 (n = 22 genes) and M2 (n = 12). Module M1 was enriched in down-regulated genes (OR = 1.88, p = 2.59 × 10−8), while module M2 was enriched in up-regulated genes (OR = 6.47, p = 2.2 × 10−16), suggesting the presence of distinct dysregulated gene-networks. Genes with predominant VT expression (14%) were mainly down-regulated, while genes with predominant EVT expression (9%) were mainly up-regulated even though both sets had the most members in module M1 (Figure 3D). These data suggested that the functions of both VT and EVT are strongly impacted in preterm preeclampsia, albeit in different ways.

11Predominantly placenta-expressed genes, down-regulated in module M1, are regulators of fetal growth (CSH1, HSD11B2), metabolism (ESRRG), estrogen synthesis (HSD17B1), and immune functions (LGALS14), some of which were reported to be down-regulated in preeclampsia. Within this module, ESRRG, POU5F1, and ZNF554 transcription regulatory genes had the highest number of significant correlations with predominantly placenta-expressed genes and, hence, deemed as hub factors (Figure S1A in Supplementary Material). Of note, these transcription factors have been implicated in the regulation of stemness and differentiation, pointing to the possible involvement of module M1 in the dysregulation of trophoblast differentiation in preterm preeclampsia. We selected ZNF554 for functional studies, since it belongs to the KRAB zinc finger family, crucial for early embryonic development and differentiation, and it may regulate genes in its co-expression network involved in biological processes affected by preeclampsia, such as development, chromatin assembly, signaling, adhesion, migration, and metabolism (Figure S1B in Supplementary Material).

12In module M2, FLT1, which expresses sFlt-1, the main driver of BP elevation in the terminal pathway of preeclampsia, was up-regulated. Moreover, module M2 genes were strongly overrepresented (OR = 29.9, p = 6.54 × 10−95) among genes that had correlated expression with mean arterial pressure (MAP) (Data S6 in Supplementary Material), suggesting a key role for this module in promoting hypertension. Within this module, BCL6, BHLHE40, and ARNT2 had the highest correlation in gene expression with predominantly placenta-expressed genes, including FLT1 (Figure S1A in Supplementary Material). Of note, these transcription regulatory hub genes are involved in hypoxia response. ARNT2 heterodimerizes with hypoxia-response regulator HIF-1α that is involved in trophoblast invasion and the pathogenesis of preeclampsia, and ARNT2 is a key regulator for adaption to hypoxic conditions at high altitudes, where the incidence of preeclampsia is much higher (16%) than at low altitude. BHLHE40 links immune and hypoxia-induced pathways. BCL6, a gene previously found to be up-regulated in the placenta in preterm preeclampsia, represses trophoblast differentiation and is regulated by stress-activated protein kinase signaling pathways. Our co-expression analysis revealed an enrichment of biological processes, such as differentiation, apoptosis, metabolism, signaling, and responses to stimuli including oxygen among genes co-expressed with BCL6 (Figure S1C in Supplementary Material). Since BCL6 is a key player in inflammation and oxygen-driven regulation of cell fate, we selected this gene for functional studies along with ARNT2.

13Next, to validate the microarray results on a larger patient population (n = 100) comprised mostly of subjects of African-American origin presenting various phenotypes of preeclampsia (Table 2; Figure S2 in Supplementary Material), we selected 47 genes for qRT-PCR profiling of 100 placentas (Table 3), provided these genes were: (1) dysregulated genes with predominant placental and syncytiotrophoblastic expression, potentially encoding biomarkers; (2) dysregulated hub transcription regulatory genes that had a high co-expression with dysregulated, predominantly placenta-expressed genes in modules M1 and M2; and (3) non-dysregulated genes with roles in trophoblast differentiation, trophoblast-specific gene expression, or the pathogenesis of preeclampsia.

14As depicted in the heatmap representing qRT-PCR data for the 100 placentas and 47 genes (Figure 4A), validation data supported microarray experiments in terms of (1) differential expression of selected genes, (2) high correlation of genes within modules M1 and M2, and (3) separate dysregulation of modules M1 and M2 from each other. Comparison between preterm preeclampsia and gestational age-matched controls revealed that qRT-PCR data validated microarray results (differential or non-differential gene expression) for 33 of 47 genes (70%), and further confirmed the differential expression of four genes that were non-dysregulated in the microarray cohort (Table 3). As expected, the extent of changes in gene expression was less pronounced in term preeclampsia (Table 3; Figure 4C; Figure S3 in Supplementary Material). As a further confirmation, there was a strong correlation between qRT-PCR data as well as tissue microarray (TMA) and immunostaining results for selected proteins in module M2 (Figures 4C–E). The extent of dysregulation was larger in preterm phenotypes of preeclampsia than in term phenotypes, in agreement with the more severe placental pathology. Among transcription regulatory genes of module M1, HLF had strong down-regulation only in preterm preeclampsia while ZNF554 was expressed at a lower level in all preeclampsia phenotypes (Figure S3 in Supplementary Material). Among the transcription regulatory genes in module M2, ARNT2, BCL6, and JUNB were highly expressed in preterm but not in term preeclampsia, suggesting that these might play a role only in the pathology of preterm cases (Figure 4C; Figure S3 in Supplementary Material). Overall, these data reflect the heterogeneous placental pathology and the more severely affected pathways in preterm preeclampsia.

15Because of the possible involvement of these two gene modules in distinct pathologic pathways, we correlated qRT-PCR data with maternal and fetal clinico-pathological indicators to further investigate their association (Figure 4B). This analysis showed that 7 of 9 genes related to the birthweight (BW) percentile were from module M1, while 10 of 15 genes related to BP were from module M2, confirming our observations with microarray data. Therefore, we decided to also refer to these as “M1-BW” and “M2-BP” modules. Of interest, the expression of most (9/14) M1 genes was negatively associated with the “maternal vascular malperfusion” score of the placenta, while the expression of most (10/14) M2 genes was positively correlated with this parameter (Data S7 in Supplementary Material). This observation fits well with the accepted concept that maternal vascular malperfusion of the placenta leads to oxidative stress and the increased placental expression of FLT1.

16Next, we investigated whether placental gene module dysregulation can be detected in maternal blood and whether liquid biopsy can be used to determine when the module dysregulation occurs. To study the relationship between placental gene expression and maternal biomarker concentrations, we selected predominantly placenta-expressed genes with the most differential expression in modules M1 (CSH1) or M2 (LEP), and correlated their placental expression levels with concentrations of their secreted protein products (hPL and leptin) in maternal blood collected in both first (n = 12) and third trimesters (n = 19) in a Hungarian patient population. Positive correlations were found between placental gene expression and maternal blood protein concentrations for both biomarkers in both trimesters (Figure S4 in Supplementary Material), suggesting that certain proteins in the maternal circulation reflect the placental expression of their encoding genes throughout pregnancy.

17To detect disease-associated protein signatures of placental dysfunction in maternal blood, similar to plasma DNA tissue mapping for noninvasive prenatal assessments, we performed “virtual” liquid biopsy of the placenta in preterm preeclampsia (Figure S5 in Supplementary Material; Figure 5). We identified five genes with predominant placental expression, which have products extensively investigated in maternal blood in preterm preeclampsia or in all cases of preeclampsia during all trimesters. From data of 61 reports, we built a database of 80,170 measurements in which preeclampsia data were expressed as a percentage of average levels in controls. Regarding biomarkers in module M1, we found hPL levels in preterm preeclampsia to be consistently below control levels throughout pregnancy. This is substantiated by the down-regulation of another M1 biomarker, HSD17B1, in first trimester maternal blood, confirming a generalized down-regulation of M1 biomarkers in early pregnancy. By contrast, the levels of biomarkers in module M2 had constant elevation in preterm preeclampsia compared to controls during gestation. When analyzing only data collected in the first trimester, prior to the onset of maternal circulation to the placenta, levels of M1 and M2 biomarkers were lower in women with preeclampsia compared to controls, while after 12 weeks of gestation, following the opening of the intervillous spaces to maternal blood flow, patients had increasing levels of M2 biomarkers in maternal blood compared to controls. These findings offer support for the conclusion that placental transcriptomic changes typical for preeclampsia in the third trimester are rooted in the first trimester.

18Next, we investigated the maternal serum proteome in early pregnancy in distinct phenotypes of preeclampsia and the potential effects of proteomic changes on the placental transcriptome. Comparing samples from women with preterm preeclampsia associated with SGA and their respective controls from an Israeli patient population (n = 10, Table 4), 19 DE protein spots were identified and investigated by mass spectrometry. According to public data, many of these proteins have a role in immune response, complement and coagulation cascades, lipid transport and metabolism, angiogenesis, BP regulation, and ion transport (Figure 6A, Data S8 in Supplementary Material). Comparing samples from women with term preeclampsia and their respective controls from this Israeli patient population (n = 10, Table 4), 14 DE protein spots could be identified (Figure S6A; Data S8 in Supplementary Material). Many of the proteins found in these spots are the same as those found in preterm preeclampsia or function in the same pathways. Of note, these pathways and the 26 differentially abundant proteins identified in term and preterm preeclampsia had mostly been implicated in a later stage of preeclampsia. We concluded that there is a common dysregulation of the maternal serum proteome in term and preterm preeclampsia; however, the extent of changes is larger in the latter, in agreement with the more fulminant and early pathogenesis.

19Supporting our findings, five studies found proteins with differential abundance in preterm preeclampsia in the same direction as our two-dimensional differential in-gel electrophoresis (2D-DIGE) study (Figure 6A). We also collected maternal blood specimens in the first trimester from a Hungarian patient population (n = 15, Table 5), and measured the concentrations of 10 immune proteins in patients with preterm preeclampsia associated with SGA and matched controls using liquid chromatography–mass spectrometry multiple reaction monitoring (MRM). In spite of the difference between the methods and ethnic background, MRM identified 4 of these 10 proteins as having differential abundance in the same direction as in the 2D-DIGE study, supporting the early pro-inflammatory changes in the maternal proteome in preterm preeclampsia (Figure 6A; Data S9 in Supplementary Material). MRM and proteomic evidence published to date supported 37% (7/19) of the proteomic changes detected by 2D-DIGE.

20To reveal whether early proteomic changes in maternal blood may affect the placental transcriptome, placental DE genes with documented connections to DE proteins in maternal serum were identified by Pathway Studio. The 121 DE placental genes with connections to DE serum proteins in preterm preeclampsia were marginally over-represented by those from module M2 (48/121, OR = 1.4, p = 0.057). Angiotensinogen had the largest number of connections to DE genes including FLT1 and LEP (n = 86, OR = 1.9, p = 4.9 × 10−5). This protein also had the largest number of connections to module M2 genes (n = 35), followed by plasminogen (n = 11) and kininogen-1 (n = 9), all involved in BP regulation (Figure 6B). These data were supported by the “renin-angiotensin signaling” as being a top pathway (p = 1.28 × 10−4) among the 35 angiotensinogen-connected DE genes (Data S10 in Supplementary Material). Similar results were obtained when analyzing connections between 116 DE placental genes and DE serum proteins in term preeclampsia (Figure S6B in Supplementary Material).

21Next, we aimed to study various factors implicated in the pathogenesis of preeclampsia in trophoblast models to determine whether these may drive the observed placental transcriptomic changes. Since the results noted above suggested that maternal serum proteins can influence the placental transcriptome, we first measured the effects of maternal blood from early pregnancy on trophoblastic gene expression. We treated primary VTs during differentiation with first trimester sera from women with preterm preeclampsia or normal pregnancy, and analyzed the expression of genes included in the placental validation study. Serum from preeclamptic women compared to controls induced up-regulation of seven placental DE genes on the first and third days of trophoblast differentiation, including FLT1 and LEP (Figure 6C). Among these, six genes were up-regulated in the placenta in term and three in preterm preeclampsia. These results support in silico findings and suggest that maternal serum factors can up-regulate FLT1 and may induce the terminal pathway.

22Next, we examined whether disturbance in VT differentiation may be reflected in the placental gene expression signature in women with different preeclampsia phenotypes. Since there were no (n = 46) or subtle (n = 2) gestational age-dependent differences in the expression of selected target or housekeeping genes in third trimester control placentas in our microarray study, the comparison of placentas and trophoblasts from various gestational ages was deemed to be valid. Thus, we performed qRT-PCR profiling of primary VTs during 7 days of differentiation to reveal the dynamics in the expression of genes (Figure 7B), which we similarly profiled in the placenta (Figure 7A). We determined whether maximum expression change of selected genes during VT differentiation compared to day 0 inversely correlated with their maximum expression change in different preeclampsia phenotypes. We found that this was the case for 17 genes (Figure 7C), suggesting the delay or inhibition of VT differentiation-related expression change of these genes in preeclampsia. Among these 17 genes, 15 showed this behavior in preterm and 9 in term preeclampsia, suggesting that VT differentiation problems are more pronounced in preterm than in term preeclampsia (p = 0.057).

23Subsequently, we tested how other factors implicated in preeclampsia pathogenesis, namely physiologic hypoxia (2% O2) or alternating hypoxic (1% O2) and normoxic (20% O2) conditions (ischemia), in combination with the overexpression of hub transcription regulatory factors, may affect trophoblastic gene expression in a widely used BeWo cell trophoblast model. Two percent O2 induced the dysregulation of only five genes, including LEP and FLT1, from the set of genes investigated in the placenta in preterm preeclampsia (Figure 8C), while alternating O2 concentrations induced the dysregulation of only three genes (Figure 8D). Since hypoxia or ischemia alone did not induce similar transcriptomic changes in BeWo cells as seen in the placenta in preeclampsia, we examined how the overexpression of hub transcription regulatory genes in module M2 may modify the effects of these two conditions. Of note, 2% O2 combined with ARNT2 or BCL6 overexpression led to the dysregulation of a large number of genes (Figure 8C). There were 9 genes (6 in the M2 and 3 in the M1 modules) dysregulated in BeWo cells, including FLT1, ARNT2, and ZNF554, similar to preeclampsia. Alternating O2 concentrations combined with ARNT2 or BCL6 overexpression led to the dysregulation of 11 genes (5 in the M2 and 3 in the M1 modules), including LEP, FLT1, and ENG, similar to preeclampsia (Figure 8D). However, ARNT2 or BCL6 overexpression at normoxic conditions did not lead to substantial gene expression changes (Figure 8B), underlining the importance of gene–environment interactions. A permutation test showed that BCL6 overexpression in ischemia mimicked overall expression changes of module M1 and M2 genes in preterm preeclampsia, while ARNT2 overexpression in ischemia (and also in hypoxia) mimicked the up-regulation of module M2 genes in term and preterm preeclampsia (Data S11 in Supplementary Material). Since BCL6 overexpression up-regulated ARNT2 both in ischemia and hypoxia but not vice versa, we propose that BCL6 is upstream of ARNT2. The up-regulation of these two transcription regulatory genes sensitize the trophoblast to ischemia, leading to the early dysregulation of modules M1 and M2, thus promoting preterm preeclampsia (Figure S7 in Supplementary Material).

24Since none of the investigated conditions could up-regulate BCL6, we wondered whether BCL6 overexpression might have an epigenetic background. Treating BeWo cells with 5-azacitidine down-regulated BCL6, supporting that its expression is regulated by DNA methylation in the trophoblast (Figure S8 in Supplementary Material). Of note, a recent study described the first intron of BCL6 to be key in its overexpression in Burkitt lymphoma via altered DNA methylation. The Human Reference Epigenome Mapping Project revealed a differentially methylated region (DMR) in this intron in H1 embryonic stem cells as well as trophoblastic and neuronal cells derived from H1 cells, suggesting that this DMR may be differentially methylated in the trophoblast compared to other cells. To address whether this intronic region may be affected in the trophoblast in preeclampsia, we investigated DNA methylation in this region in primary VTs compared to cord blood cells collected from the same normal pregnancies. Bisulfite sequencing showed that Chr3:187,458,083-187,458,651 and Chr3:187,460,304-187,460,374 regions contain 12 hypermethylated CpGs in the trophoblasts compared to cord blood cells (Figures S8 and S9 in Supplementary Material). Further, we tested DNA methylation in this region in 100 placentas with qRT-PCR data including patients with preeclampsia after laser capturing VTs. Three CpGs (Chr3:187,458,095, Chr3:187,458,163, and Chr3:187,458,327) were differentially methylated in preeclampsia (Figure S10 in Supplementary Material), of which CpG Chr3:187,458,163 was differentially methylated in preterm preeclampsia, suggesting that this CpG may have a role in BCL6 dysregulation in preeclampsia.

25Next, we were interested in how the dysregulation of module M1 genes may play a role in preeclampsia pathology. Among hub genes of this module, ZNF554 was of most interest due to the biological processes enriched in its co-expression network (Figure S1B in Supplementary Material), and also to its potential placenta- and preeclampsia-related regulation by transposable elements. This hypothesis was based on the fact that insertion of transposable elements into regulatory regions can lead to transcriptional changes, especially in the placenta, and that the 5′ flanking region of ZNF554 contains many LTR10A copies which had top enrichment among module M1 genes (OR = 17.4, p = 1.27 × 10−7; Data S12 in Supplementary Material). Of note, LTR10A drives placenta-specific expression of NOS3, and it may also have a similar effect on ZNF554. Indeed, ZNF554 had the highest expression in the placenta in comparison to 47 other human tissues, which mostly had negligible ZNF554 expression (Figure 9A).

26Subsequent in situ hybridization (Figure 9B) and immunostaining (Figure 9C) of first and third trimester placentas of women with a normal pregnancy showed dominant ZNF554 expression in the syncytiotrophoblast but not in the villous cytotrophoblast. Thus, we investigated ZNF554 expression during VT differentiation, in which it was up-regulated similar to CSH1 (Figure 9D), supporting that ZNF554 expression is developmentally regulated in VTs. Of interest, ZNF554 immunostaining was faint in the syncytiotrophoblast in preeclampsia compared to controls (Figures 9E,F). To characterize the loss of syncytiotrophoblastic ZNF554 function, we silenced ZNF554 in BeWo cells. At 74% ZNF554 knock-down (p = 5.24 × 10−6) (Figure 9G), decreased nuclear and cytoplasmic ZNF554 immunostaining was found (Figure 9H). Microarray analyses of ZNF554-silenced cells revealed 123 DE genes (Data S13 in Supplementary Material) including 9 DE placental genes in preeclampsia, and the dysregulation of the “glycolysis/gluconeogenesis” pathway (OR = 7.8, q = 0.06) as well as 18 molecular functions including “RNA binding” (down) and “activin binding” (up) (Figure 9I; Data S14 in Supplementary Material). The up-regulation of FSTL3 was confirmed by qRT-PCR (2.7-fold, p < 0.001) (Figure 9J). FSTL3 encodes a secreted glycoprotein that inactivates activin and other TGFβ ligands. It is involved in the regulation of EVT invasion and its placental up-regulation is associated with low BW in preeclampsia. This finding confirms that the dysregulation of ZNF554 may have key downstream effects on the pathogenesis of preeclampsia.

27We supposed that ZNF554 may also affect EVTs, since its expression in EVTs was detected in first and third trimester maternal decidua (Figures 10A–D), and ZNF554-positive intraluminal and endovascular trophoblasts were found in the wall of transformed spiral arteries (Figure 10B). In preterm preeclampsia, ZNF554 immunostaining of EVTs was weaker than in controls (Figures 10C,D). To characterize the loss of ZNF554 function in EVTs, we silenced ZNF554 in trophoblastic HTR8/SVneo cells. At 87% knock-down (p < 0.001) (Figure 10E), we observed decreased nuclear and cytoplasmic ZNF554 immunostaining (Figure 10F). Microarray analysis of ZNF554-silenced cells showed 185 DE genes (Data S15 in Supplementary Material) including 18 DE placental genes in preeclampsia. Gene ontology (GO) analysis revealed 16 molecular functions dysregulated, including “cyclin-dependent protein kinase regulator activity,” “metalloendopeptidase inhibitor activity,” and “insulin-like growth factor binding.” The 67 enriched biological processes included “regulation of growth,” “smooth muscle cell migration,” “smooth muscle cell-matrix adhesion,” and “response to oxygen levels,” all relevant to trophoblast invasion and placental pathology of preeclampsia (Data S16 in Supplementary Material).

28qRT-PCR confirmed the dysregulation of eight DE genes. Two genes (CDKN1A, STK40) are involved in the regulation of cell proliferation and differentiation, and proliferation assays showed that ZNF554 knock-down decreased cell proliferation in HTR8/SVneo cells slightly after 48 h (−14%, p = 0.02) (Figure S11 in Supplementary Material). Six genes (FSTL3, ITGB5, MYL12A, SDC1, SERPINE1, and TIMP3) encode proteins involved in cell adhesion, migration, invasion, and angiogenesis (Figure 10G). Since EVTs move through an environment with changing O2 levels, we used O2 concentrations for conditions relevant for endovascular (8%) and interstitial (2%) trophoblast invasion besides standard cell cultures (20%). The effect of ZNF554 knock-down was significant regardless of O2 levels on four genes (ITGB5, MYL12A, SERPINE1, and TIMP3), while there was an interaction between O2 levels and ZNF554 silencing on two genes (FSTL3 and SDC1) (Figure 10G).

29The up-regulation of SERPINE1 (PAI-1) and tissue inhibitor of metalloproteinases-3 (TIMP-3) was also confirmed at the protein level in supernatants of ZNF554-silenced cells (Figure 10H). Both proteins have an inhibitory function on trophoblast migration and invasion, and TIMP-3 is the major tissue metalloproteinase inhibitor at the maternal–fetal interface, which is up-regulated in preeclampsia. These results suggested that ZNF554-silenced cells have reduced migratory and invasive functions. Indeed, functional assays revealed that ZNF554 silencing had a strong inhibitory effect on trophoblast migration (p = 1.9 × 10−10) regardless of the O2 concentration, and also on invasion (p < 0.001), especially at 2% O2 concentration (Figure 10I). These data corroborated that ZNF554 supports trophoblast invasion via modulating a set of key genes that are involved in this process.

30We wondered whether placental ZNF554 down-regulation might have an epigenetic background, given that the ZNF554 flanking region contains several transposable elements, including several Alus, which are generally hypomethylated in germ cells and the placenta, while their hypermethylation leading to altered gene expression may be detected in preeclampsia. Thus, the down-regulation of ZNF554 expression in the placenta of patients with preeclampsia may also be reflected in the DNA methylation of the transposable elements in its 5′ flanking region.

31The treatment of BeWo cells with 5-azacitidine increased ZNF554 expression, showing the role of DNA methylation in trophoblastic ZNF554 regulation (Figure 11A). The subsequent search in the Human Reference Epigenome Mapping Project data revealed a DMR located in the AluY, which was hypomethylated in H1 embryonic stem cells and H1-derived trophoblasts compared to H1-derived neuronal cells (Figure 11B). This was of interest, since AluY is a retrotransposon evolved recently in primates, and its differential DNA methylation supports the expression of other gene transcripts in the placenta compared to somatic tissues. These data prompted us to investigate the DNA methylation in this genomic region in primary VTs and cord blood cells collected from the same normal pregnancies. In fact, bisulfite sequencing showed that the AluY, similarly to the AluSq2, is heavily methylated in cord blood cells compared to the hypomethylated trophoblast, suggesting its importance in the developmental regulation of ZNF554 expression (Figure 11B; Figure S12 in Supplementary Material).

32Further, we tested DNA methylation in this region in 100 placentas with qRT-PCR data, including patients with preeclampsia after laser capture of the VTs. Bisulfite sequencing of the trophoblastic DNA revealed four CpGs on AluY (Chr19:2,818,823, Chr19:2,818,864, Chr19:2,818,868, and Chr19:2,818,876) hypomethylated in controls but hypermethylated in preterm preeclampsia, with highest methylation in cases associated with SGA (Figure 11C; Figure S13 in Supplementary Material). Importantly, we found correlations between Chr19:2,818,823 CpG methylation and ZNF554 expression (R = −0.30, p = 0.04), maternal vascular underperfusion score of the placenta (R = 0.36, p = 0.03), and BW percentile (R = −0.41, p < 0.01) (Figure 11D). These data collectively provide evidence that the hypermethylation of certain CpGs in AluY in the trophoblast may result in low ZNF554 expression, impaired trophoblast invasion, preeclampsia, and fetal growth restriction.

33The placenta has a key role in the pathogenesis of the terminal pathway of preeclampsia, which may be triggered by discrete disease pathways at early stages of pregnancy, leading to the development of different preeclampsia phenotypes. In this study, an integrated systems biology approach was employed to gain insights into these complex pathways, given that this strategy offered the ultimate analytic solution to investigate and understand the complex disease pathways of the syndrome of preeclampsia. We incorporated “omics,” clinical, placental, and functional data from patients with distinct phenotypes of preeclampsia. We employed molecular network-based approaches to identify networks and modules of genes or proteins that are perturbed in the placenta and the maternal circulation of women with preeclampsia.

34Our placental transcriptomics study identified 1,409 DE genes involved in many biological processes (e.g. BP regulation, apoptosis, development, hormone secretion, metabolism, and signaling) that were previously implicated in the pathogenesis of preterm preeclampsia by other placental transcriptomics studies. Despite the differences in patient populations, design, or methodologies between these studies, many DE genes on our list have also been found by other groups. Indeed, from the 40-gene meta-signature that characterized the significant intersection of DE genes from independent placental gene signatures in preeclampsia in the meta-analysis of Kleinrouweler et al., our microarray and qRT-PCR studies found 26 (65%) to be DE in preterm preeclampsia. Of note, 16 of these 26 genes belong to module M2, while only six to the M1 module. This supports our observation that the dysregulation of module M2 is associated with BP elevation, the maternal disease condition required for patient inclusion into these studies. The weaker involvement of module M1 genes in the meta-signature may reflect the heterogeneity of preeclampsia transcriptomics studies regarding fetal (growth restriction) and placental disease conditions. Our microarray study was homogeneous for preterm preeclampsia cases with low BW and placental disease, while only a couple of other studies had this rigor. To overcome the inconsistency of smaller placental transcriptomics studies, Leavey et al. employed advanced bioinformatics methods to aggregate microarray datasets across multiple platforms to generate large datasets of patient samples. Unsupervised clustering of these datasets revealed three distinct molecular subclasses of preeclampsia. Among these, the “canonical” subclass, which is associated with the differential expression of our module M2 genes, was characteristic for preterm preeclampsia and consistent with stress response to poor oxygenation, further supporting our findings. However, our discovery on the two major dysregulated placental disease gene modules and their hub transcription regulatory genes, separately associated with maternal or fetal disease pathways, are novel (Figure 12).

35Functional assays on hub transcription factors of these two disease gene modules demonstrated that ZNF554 (M1) down-regulation leads to impaired trophoblast invasion, while BCL6 and ARNT2 (M2) overexpression sensitizes the trophoblast to ischemia, which are hallmarks in the pathogenesis of preterm preeclampsia. In the “ZNF554” pathway, hypermethylation of AluY in the ZNF554 5′ flanking region inhibits gene expression, leading to impaired trophoblast invasion, placental vascular malperfusion, and low BW. In the “BCL6-ARNT2” pathway, which is activated only in preterm preeclampsia, ischemic stress of the trophoblast coupled with BCL6 and ARNT2 overexpression increases FLT1 expression. This then eventually promotes the anti-angiogenic state, hypertension, and the early onset of this syndrome (Figure S7 in Supplementary Material).

36The perturbed placental disease gene modules in preterm preeclampsia can be detected by liquid biopsy in maternal blood in the early stages of pregnancy. Indeed, we detected the down-regulation of M1 and M2 disease gene module biomarkers in these patients during the first trimester. Of interest, the up-regulation of module M2 biomarkers can be detected after the establishment of maternal circulation in the intervillous space. These findings support the observation that placental transcriptomic changes, typical for preterm preeclampsia observed in the third trimester, are rooted in the first trimester. The positive correlation of gene expression in module M1 with BW and in module M2 with BP suggests that M1 genes may be biomarkers for placental and fetal growth and development while M2 genes can serve as biomarkers for placental stress.

37First trimester maternal blood proteomics uncovered the altered abundance of proteins of the renin–angiotensin and immune systems as well as complement and coagulation cascades in patients who subsequently developed both preterm and term preeclampsia. The same proteins and pathways were found to be dysregulated in maternal blood in later stages of preeclampsia by other proteomics studies; however, ours is the first revealing their dysregulation at an earlier stage of pregnancy when there is no or minimal direct connection between the placenta and the maternal circulation.

38From these dysregulated maternal serum proteins, in silico analysis pointed to candidates, which may drive trophoblastic transcriptomic changes, and corroborated earlier findings on angiotensinogen/angiotensin II in driving hypertension indirectly through FLT1 up-regulation in addition to its direct effects. Moreover, in vitro functional assays revealed that altered maternal serum proteome in the first trimester can affect the trophoblastic transcriptome and up-regulate FLT1. This is in agreement with reports indicating that maternal blood factors in preeclampsia can induce trophoblastic soluble endoglin overexpression and the development of preeclampsia-like symptoms in mice.

39Remarkably, most of the dysregulated maternal serum proteins in the first trimester in both preterm (11 of 19) and term (7 of 14) preeclampsia are implicated in immune functions (Figure 6; Figure S6 in Supplementary Material), suggesting a critical role for immune pathways and inflammation in the early pathogenesis of both phenotypes of preeclampsia. This is consistent with clinical, epidemiological, and immunological evidence showing that: (1) preeclampsia has multiple risk factors (e.g. dyslipidemia, hypertension, insulin resistance, and obesity) characterized by heightened inflammation; (2) the combination of inhibitory decidual NK (dNK) cell killer immunoglobulin-like receptor and the fetal HLA-C2 genotype increases the susceptibility to preeclampsia due to the loss of activating interactions between trophoblasts and dNK cells at early stages of placentation; (3) an altered local immune regulation and a shift toward the pro-inflammatory macrophage phenotype promotes a pro-inflammatory milieu in the maternal decidua in preeclampsia; (4) an imbalance between Th1/Th2/Th17/Treg cells in preeclampsia leads to failure of maternal–fetal tolerance mechanisms; and (5) complement system activation in preeclampsia leads to the activation of innate immune cells and placental damage. The role of inflammation in early preeclampsia disease pathways is also supported by in vivo studies showing that bacterial endotoxin administration to pregnant rats induces placentation defects and symptoms consistent with preeclampsia.

40Overall, our data show that there are distinct maternal and placental disease pathways, and their interaction influences the clinical presentation of preeclampsia. The activation of maternal disease pathways can be detected in both preterm and term preeclampsia earlier and upstream of placental dysfunction, not only downstream as described before, and distinct placental disease pathways are superimposed on these maternal pathways. This is a paradigm shift in our understanding of preeclampsia, which in agreement with epidemiological studies warrants for the central pathologic role of preexisting maternal diseases or perturbed maternal–fetal–placental immune interactions in preeclampsia.

41The superimposed placental disease pathways differ between preterm and term preeclampsia. For preterm preeclampsia, our functional data suggest that placental disease pathways are partly originated from altered trophoblast differentiation, which is followed by trophoblastic stress, induced by perturbed maternal blood proteome factors and/or ischemia after the onset of maternal circulation to the placenta. Our data are consistent with recent views indicating that defects in trophoblast proliferation, differentiation, invasion, and plugging are associated with defective decidualization, decidual inflammation, and the disturbance in endometrial-trophoblast dialog during the peri-conception period. Abnormal trophoblast invasion and plugging will subsequently lead to the aberrant onset of maternal circulation and malperfusion, causing placental oxidative stress in preterm preeclampsia. On the other hand, maternal disease pathways induce mainly placental dysfunction without maldevelopment in term preeclampsia. This is substantiated by the differences observed in the maternal proteome, placental transcriptome, and trophoblastic DNA methylation between term and preterm preeclampsia in our study. Moreover, this is also consistent with the major differences between these two preeclampsia phenotypes in etiology, placental histopathology, and stress levels as well as clinical presentation.

42Our findings are very timely in the light of recent clinical research showing that the administration of aspirin before 16 weeks of gestation to pregnant women at risk for preeclampsia prevents the preterm phenotype of this syndrome. Thus, the anti-inflammatory and anti-platelet actions of aspirin may ameliorate the early pro-inflammatory disease pathways leading to placental maldevelopment in preterm preeclampsia. Based on our discovery of these novel disease pathways and their hub molecules, we propose a “molecular phase” of preeclampsia (Figure 12), where early pathologic events can already be detected by maternal blood biomarkers, offering non-invasive diagnostics of maternal and placental disease pathways. Biomarkers of these disease pathways may open new venues for the molecular characterization of patients destined to develop preeclampsia, using multi-biomarker profiles that support preventive approaches for patients with distinct preeclampsia phenotypes.

43Placental tissue and maternal blood samples were collected from Caucasian women at the First Department of Obstetrics and Gynaecology, Semmelweis University in Budapest, Hungary as described previously. Pregnancies were dated according to ultrasound scans between 8 and 12 weeks of gestation. Patients with multiple pregnancies or fetuses having congenital or chromosomal abnormalities were excluded. The collection and investigation of human clinical samples were approved by the Health Science Board of Hungary. Written informed consent was obtained from women prior to sample collection, and the experiments conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report. Specimens and data were stored anonymously.

44Women were enrolled in the following groups: (1) preterm severe preeclampsia, with or without HELLP (hemolysis, elevated liver enzymes, low platelets) syndrome (n = 12) and (2) preterm controls (n = 5) (Table 1). Preeclampsia was defined according to the criteria set by the American College of Obstetricians and Gynecologists, and subdivided into preterm (<37 weeks) or term (≥37 weeks) groups. Severe preeclampsia was defined according to Sibai et al.. Preterm controls had no medical complications, clinical or histological signs of chorioamnionitis, and delivered neonates with a BW appropriate-for-gestational age (AGA). SGA was defined as neonatal BW below the 10th percentile for gestational age. Cesarean delivery was performed in all cases due to severe symptoms as well as in all controls due to previous Cesarean delivery or malpresentation.

45Placental tissue specimens were processed immediately after delivery as described previously. For the microarray study, 1 × 1 cm villous tissue samples were excised from central cotyledons close to the umbilical cord to reduce the possible bias due to regional differences in gene expression. These tissue blocks were then dissected from the choriodecidua on dry ice, snap-frozen, and stored at −80oC. For histopathologic evaluations, five representative tissue blocks were taken from each placenta to include central and peripheral cotyledons and the maternal side of the placenta with the fetal membranes. These blocks were embedded in paraffin after fixation in 10% neutral-buffered formalin (FFPE). Maternal blood samples were obtained at the time of admission into the delivery room; aliquots of maternal sera and plasma were stored at −80oC.

46Placental specimens were examined according to a standard protocol, describing the topography and size of macroscopic lesions. Four micrometer sections were cut from the five FFPE blocks and mounted on SuperFrost/Plus slides (Fisherbrand, UK). After deparaffinization, slides were rehydrated, stained with hematoxylin and eosin, and examined in 10 randomly chosen microscopic fields using bright-field light microscopy by a pathologist blinded to the clinical information. Macroscopic and microscopic lesions were defined according to published criteria.

47Tissues were homogenized using a ThermoSavant FastPrep FP120 Homogenizer (Thermo Scientific, Wilmington, DE, USA) with Lysing MatrixD (MP Biomedicals, Illkirch, France). Total RNA was isolated using RNeasy Fibrous Tissue Mini Kit (QIAGEN GmbH, Hilden, Germany), quantified with NanoDrop 1000 (Thermo Scientific), and assessed by Agilent 2100 Bioanalyzer (Matriks AS, Oslo, Norway). Total RNAs were labeled, and Cy3-RNAs were fragmented and hybridized to the Whole Human Genome Oligo Microarray G4112A (Agilent Technologies, Santa Clara, CA, USA) on an Agilent scanner, and processed with Agilent Feature Extraction software v9.5 according to the manufacturer’s guidelines.

48Demographics data were compared by the Fisher’s exact test and Mann–Whitney test using SPSS version 12.0 (SPSS Inc., Chicago, IL, USA). Microarray data analysis was performed using the R statistical language and environment1 following the MIAME guidelines and methodologies described previously. Microarray expression intensities were background corrected using the “minimum” method in the “backgroundCorrect” function of the “limma” package. After log2 transformation, data were quantile-normalized. From the 41,093 probesets on the array, 93 were removed before differential expression analysis, lacking annotation in the array definition file (Agilent Technologies). Subsequently, an expression filter was applied to retain probesets with intensity greater than log2(50) in at least two samples, yielding a final matrix of 30,027 probesets (15,939 unique genes). Differential gene expression was assessed using linear models including adjustment for batch effects, with coefficient evaluation via moderated t-tests. P-values were adjusted using the false discovery rate (FDR) method and Benjamini–Hochberg correction to compute q-values. Since there were no differences among patient groups in maternal age or ethnic background, we did not adjust for these parameters. Target gene Entrez IDs for the probesets were determined using the R package “hgu4112a.db.” For probesets without annotation in the package, Entrez IDs were taken from the array definition file (Agilent Technologies). Probesets remaining un-annotated (without Entrez ID and/or gene symbol) were removed from further analysis. Probesets were defined as DE (n = 1,409) if they had a q ≤ 0.2 and a fold-change ≥1.5 (Data S1 in Supplementary Material).

49From the DE genes in preeclampsia, those encoding for proteins with functions in transcription regulation (n = 137) were identified using the Metacore (GeneGo Inc., Saint Joseph, MI, USA) and GeneCards v32 databases.

50We downloaded the human U133A/GNF1H microarray data on 79 human tissues, cells, and cell lines from the Symatlas/BioGPS database to identify human genes with predominant placental expression. A probeset was defined as having predominant placental expression if its placental expression was (1) ≥1,000 fluorescence units; (2) six times higher than the median value in 78 other tissue and cell sources; and (3) two times higher than its expression in the tissue with the second highest expression. The resulting 215 probesets corresponded to 153 unique genes. Eleven additional genes that were not present on the microarray platform used by Symatlas/BioGPS (Affymetrix, Santa Clara, CA, USA) were added to this list based on previously published evidence of placenta-specific expression (Data S2 in Supplementary Material). Of the 164 predominantly placenta-expressed genes, 157 were present in our placental microarray data. These genes were tested for enrichment in DE genes compared to all genes on the array (1,409 of 15,939) using the Fisher’s exact test.

51The expression levels of DE genes in EVT compared to VT lineages were analyzed by retrieving published microarray datasets and reanalyzing expression data. Raw Affymetrix GeneChip Human Genome U133A 2.0 Array data from Bilban et al. and Tilburgs et al. was downloaded from GEO (GSE9773) and ArrayExpress (E-MATB-3217) respectively, and was processed with the “affy” and “limma” packages of Bioconductor.3 After RMA normalization, log-fold changes were calculated. Processed Illumina Human HT-12 V3 BeadArray data from Apps et al. was downloaded from ArrayExpress (E-MATB-429), and then log-fold changes were calculated.

52Chromosomal locations for all genes tested on the Agilent array were obtained from the R package “org.Hs.eg.db.” Of the 15,939 unique and 1,409 DE genes on the array, 15,935 and 1,408 could be assigned to chromosomes, respectively. Mapping the microarray probesets on the Affymetrix human U133A/GNF1H chips to ENTREZ identifiers was performed using the Bioconductor “hgu133a.db” and “hgfocus.db” packages. Chromosomal locations of the resulting list of genes were obtained from the package “org.Hs.eg.db” and from the National Center for Biotechnology Information for the 11 additional genes. Enrichment analyses for chromosomes among predominantly placenta-expressed genes, DE genes and DE genes encoding for transcriptional regulators (Data S3–S5 in Supplementary Material) were tested by the Fisher’s exact test. Chromosomal locations of these genes were visualized by Circos.

53Weighted gene co-expression network analysis was applied on the 1,409 DE genes across 17 samples to identify distinct regulatory modules and prioritize candidate genes for qRT-PCR validation. A gene pair-wise similarity (absolute Pearson correlation) matrix was first computed, then soft-thresholded by raising it to the power of 10 (chosen based on the scale-free topology criterion) to obtain an adjacency matrix. The topology overlap matrix (TOM) was then derived from the adjacency matrix. The topology overlap measures the node interconnectedness within a network and was generalized to WGCNA. This measure defines the similarity between the two genes based on both correlations within themselves and outside other genes. Gene distance matrix was defined as 1-TOM and used for average linkage hierarchical clustering. A hybrid dynamic tree-cutting method was applied to obtain modules (tree clusters).

54Gene modules identified with this approach were further tested for enrichment in DE genes using the Fisher’s exact test. Transcription regulatory genes expressed at high levels (average log2 intensity >9) and co-expressed (absolute Pearson coefficient >0.7) with the most genes among DE genes were treated as candidates for hub genes in the module. Hub genes then were selected based on the number and strength of their Pearson co-expression partners as well as their and their networks’ biological activities. The networks of biological processes enriched among genes co-expressed with hub factors modules were created by BINGO and visualized with Cytoscape.

55Enrichment analysis of transposable elements present in the 10,000 bp upstream region of DE genes was performed separately for the M1 (green) and the M2 (red) modules in preeclampsia versus all genes present on the microarray using the Fisher’s exact test. The locations of transposable elements and their families and classes were obtained from the “RepeatMasker” table of the UCSC Table Browser.4 P-values <0.05 were considered significant.

56Third trimester placentas (n = 100) collected predominantly from African-American women were retrieved from the Bank of Biological Specimens of the Perinatology Research Branch (Detroit, MI, USA). Pregnancies were dated according to ultrasound scans between 8 and 12 weeks. Patients with multiple pregnancies or fetuses having congenital or chromosomal abnormalities were excluded. The use of biological specimens and clinical data for research purposes was approved by the Wayne State University Human Investigation Committee and the Institutional Review Board of the Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD, NIH, DHHS). Written informed consent was obtained from women prior to sample collection, and the experiments conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report. Specimens and data were stored anonymously.

57The following homogenous patient groups were selected from a cohort: (1) preterm severe preeclampsia (n = 20); (2) preterm severe preeclampsia associated with SGA (n = 20); (3) preterm controls (n = 20); (4) term severe preeclampsia (n = 10); (5) term severe preeclampsia associated with SGA (n = 10); and (6) term controls (n = 20). Term controls, consisting of normal pregnant women with (n = 10) or without (n = 10) labor, and preterm controls with preterm labor and delivery had no medical complications or clinical or histological signs of chorioamnionitis, and delivered an AGA neonate (Table 2; Figure S2 in Supplementary Material).

58Labor was defined by the presence of regular uterine contractions at a frequency of at least two contractions every 10 min with cervical changes resulting in delivery. Preeclampsia was defined according to the criteria set by the American College of Obstetricians and Gynecologists and was subdivided into preterm (<37 weeks) or term (≥37 weeks) groups. Severe preeclampsia was defined according to Sibai et al.. SGA was defined as neonatal BW below the 10th percentile for gestational age. Cesarean delivery was performed in cases due to severe symptoms and in controls due to previous Cesarean delivery or malpresentation.

59Third trimester placental tissue specimens (n = 100) were processed immediately after delivery. For histopathologic evaluations, representative tissue blocks were taken from each placenta to include central and peripheral cotyledons, maternal side of the placenta, fetal membranes, and the umbilical cord. These blocks were embedded in paraffin after fixation in 10% neutral-buffered formalin (FFPE). Villous tissue blocks from central cotyledons were selected for TMA. For gene expression studies, systematic random sampling was used to obtain villous tissues to reduce the possible bias due to regional gene expression differences. Excised tissue blocks were homogenized and mixed in TRIzol reagent (Life Technologies), snap-frozen with liquid N2, and stored at −80oC.

60Placental specimens (n = 100) were examined according to a standard protocol, describing the topography and size of macroscopic lesions. Five µm sections were cut from the representative FFPE tissue blocks, stained with hematoxylin and eosin, and examined using bright-field light microscopy by two anatomic pathologists blinded to the clinical information. Histopathologic changes of the placenta were defined according to published criteria proposed by the Perinatal Section of the Society for Pediatric Pathology. “Maternal vascular malperfusion score” was calculated by summing the number of pathologic lesions consistent with this lesion category present in a given placenta.

61Total RNA was isolated from snap-frozen third trimester placental villous tissues (n = 100) with TRIzol reagent (Life Technologies) and RNeasy kit (QIAGEN, Valencia, CA, USA) according to the manufacturers’ recommendations. The 28S/18S ratios and the RNA integrity numbers were assessed using an Agilent Bioanalyzer 2100, and RNA concentrations were measured with NanoDrop 1000. Total RNA (500 ng) was reverse transcribed with the High Capacity cDNA Reverse Transcription Kit using random hexamers (Applied Biosystems, Foster City, CA, USA). TaqMan Assays (Applied Biosystems; Table 3, Data S17 in Supplementary Material) were used for high-throughput gene expression profiling on the Biomark qRT-PCR system (Fluidigm, San Francisco, CA, USA) according to the manufacturers’ instructions.

62Tissue microarrays were constructed from central tissue blocks of third trimester FFPE placentas (n = 100) as described earlier. Briefly, three 20 × 35 mm recipient blocks were made of Paraplast X-Tra tissue embedding media (Fisher Scientific, Pittsburgh, PA, USA). One mm diameter cores from tissue blocks were transferred in triplicate into recipient paraffin blocks using an automated tissue arrayer (Beecher Instruments, Inc., Silver Spring, MD, USA). Five µm sections cut from TMAs were placed on silanized slides and stained with antibodies and reagents (Data S18 in Supplementary Material) either on Ventana Discovery or Leica BOND-MAX autostainers.

63Tissue microarray immunostainings were semi-quantitatively scored by three examiners blinded to the clinical information with an immunoreactive score modified from a previously published one. Immunostaining intensity was graded as follows: 0 = negative, 1 = weak, 2 = intermediate, and 3 = strong. All villi in a random field of each of the three cores were evaluated by all examiners, and scores within each core were averaged to represent the target protein quantity of that core. Thus, each placenta had three scores corresponding to the three cores examined.

64Demographics data were compared by the Fisher’s exact test and Mann–Whitney test using SPSS version 12.0 (SPSS). All other data were analyzed in the R statistical environment (see text footnote 1).

65Data were analyzed using the ΔΔCt method. Data were first normalized to the reference gene (RPLP0), and the batch effect was adjusted through calibrator samples. Log2 mRNA relative concentrations were obtained for each sample as −ΔCt(gene) = Ct(RPLPO) − Ct(gene). The surrogate gene expression values (−ΔCtgene) were used to perform a hierarchical clustering with 1-Pearson correlation distance and average linkage. Between-group comparisons (in which groups were predefined based on the clinical characteristics of the patients) were performed by fitting a linear model on −ΔCt values, using the group variable indicator and the maturity status of the fetus (term vs. preterm) as covariates while allowing for an interaction effect (Figure 4; Figure S3 in Supplementary Material).

66The association between qRT-PCR gene expression and the “maternal vascular malperfusion” score was tested using a linear model. P-values of <0.05 were considered significant.

67Group comparisons using immunoscores were conducted in the same way as for the qRT-PCR data (Figure 4).

68To reveal whether the expression of any gene on the microarray was correlated with MAP while controlling for BW, a linear model (y~MAP + BW + Batch) was fit for every gene on the array, in which y represents gene expression, and the dependent variables represent MAP, BW, and batch, respectively. A moderated t-test was used to obtain p-values, which then were adjusted using the FDR method and Benjamini–Hochberg correction for multiple testing. Significance was determined using a q ≤ 0.2. Gene modules were also tested for the enrichment in genes with their expression correlated with MAP using the Fisher’s exact test (Data S6 in Supplementary Material).

69In the qRT-PCR validation study, we extended our analysis to include all 100 patients to test for the association between gene expression and mean arterial BP as well as BW percentile while adjusting for gestational age. All variables were treated as continuous. P-values of <0.05 were considered significant (Figure 4).

70Fifteen µm sections were cut from snap-frozen placentas that were also used for qRT-PCR expression profiling (n = 100) on 2 µm Glass Foiled PEN slides (Leica Microsystems). The trophoblast layer of 300–350 villi in each specimen was laser captured by a Leica DM6000B microscope (Leica Microsystems) into 0.5 ml microcentrifuge tubes. The captured material was digested with Proteinase K (PicoPure DNA Extraction Kit, Applied Biosystems) at 56°C by overnight incubation. Digestions were stopped at 95°C, and samples were stored at −70°C until DNA isolation.

71Genomic DNA was isolated from laser captured VTs (n = 100), from primary VTs (n = 3) collected for functional studies described below, and from respective umbilical cord blood leukocytes (n = 3) taken from the same pregnancies. The EZ1 Advanced Nucleic Acid Isolation System using EZ1 DNA Tissue and EZ1 DNA Blood Kits (QIAGEN) were utilized for DNA isolation, and DNA samples were quantified with Quantifiler Human DNA Quantification Kit (Applied Biosystems) according to the manufacturers’ instructions.

72The methylation status of CpGs in human H1 embryonic stem cells, H1-derived trophoblast cultured cells, and H1-derived neuronal progenitor cultured cells obtained by whole-genome shotgun bisulfite sequencing (University of California, San Diego, USA; UCSD Human Reference Epigenome Mapping Project) were visualized by the Epigenome Browser5 and used for the selection of regions of interest. Primer design, targeted amplification, and sequencing were conducted using the targeted sequencing service protocol of Zymo Research Corporation. For targeted bisulfite sequencing, 30 primer pairs were designed and validated. Primers were synthesized by Integrated DNA Technologies, Inc., (Coralville, IA, USA) and underwent quality control, which included duplicate testing for specific amplification of 1 ng bisulfite DNA using bisulfite converted human DNA. Quality control criteria included robust and specific amplification (Cp values <40 cycles and CV <10% for duplicates) of the bisulfite primers on a LightCycler 480 real-time qRT-PCR instrument (Roche Diagnostics Corp. Indianapolis, IN, USA).

73Genomic DNA samples from laser captured VTs, primary VTs, and umbilical cord blood cells as well as control samples were subjected to sodium bisulfite treatment using the EZ DNA Methylation-Direct Kit (Zymo Research Corporation). For non-methylated control, human DNA was extracted and purified with Quick-gDNA Miniprep Kit (Zymo Research Corporation) from the HCT116 cell line (American Type Culture Collection, Manassas, VA, USA), which is double knock-out for both DNA methyltransferases DNMT1 (−/−) and DNMT3b (−/−), and thus contains a low level (<5%) of DNA methylation. For methylated control, human DNA was purified similarly from the HCT116 cell line and was enzymatically methylated at all cytosine positions comprising CG dinucleotides by CpG Methylase (Zymo Research Corporation).

74Bisulfite-treated samples and 30 validated primer pairs were subjected to targeted amplification on the 48.48 Access Array System (Fluidigm), using the targeted sequencing service protocol of Zymo Research Corporation. Fluidigm’s protocols were used for sample loading, harvesting, and pooling, 1:100 dilution of amplicon pools for each sample, and for amplification using barcoded adapter-linkers (Fluidigm). Reactions were cleaned up using the DNA Clean and Concentrator-5 (Zymo Research Corporation), and products were normalized by concentration and pooled. The sequencing library was denatured, diluted, and sequenced with a 150-base paired-end run on the MiSeq Benchtop Sequencer (Illumina) according to Illumina’s protocols.

75Sequence reads from bisulfite-treated libraries were identified using standard Illumina base-calling software, and then analyzed using a Zymo Research Corporation proprietary analysis pipeline. Residual cytosines (Cs) in each read were first converted to thymines (Ts), with each such conversion noted for subsequent analysis. Reads were aligned by Bismark, a Bowtie-based alignment tool for bisulfite converted reads.6 The number of mismatches in the induced alignment was then counted between the unconverted read and reference, ignoring cases in which a T in the unconverted read was matched to a C in the unconverted reference. For a given read, only the best scored alignment was kept. If there were more than one best read, then only one was kept arbitrarily. The methylation level of each sampled cytosine was estimated as the number of reads reporting a C, divided by the total number of reads reporting a C or T. CpGs with coverage of less than four reads were removed from the analysis. The developed sensitive and robust bisulfite sequencing assays yielded a median total sequencing read of 533 (range: 30–1,725) per CpG in the trophoblast-fetal blood cell comparison and a median total sequencing read of 136 (range: 4–2,609) per CpG in the clinical sample comparison.

76Multiple sequencing counts (total and methylated) were summed for each sample at each CpG site, and samples with a total count <4 were dropped from the analysis. The mean methylation ratio in each group was computed for genomic visualization. In the comparison of methylation levels between trophoblasts and cord blood cells, a group sample size of two was considered as a minimum. In order to fit the count data, we used a generalized linear model of Poisson distribution with log link. When all samples in any of the two groups being compared had zero methylation counts, the maximum likelihood estimation of the Poisson model went to infinity. In such cases, the Student’s t-test was used alternatively. P-values and the group difference in methylation ratios were included above the bar plots in Figures S9 and S12 in the Supplementary Material. For comparisons of methylation levels between the clinical groups, only comparisons with a minimum group sample size of four were considered, and the Wilcoxon rank-sum test was used. P-values and the group difference in methylation ratios were included above the bar plots in Figure 11 as well as Figures S10 and S13 in the Supplementary Material. Differential methylation was considered to be mild, moderate, or strong when the p-value was <0.05 and the difference in methylation ratios was ≥0.125, ≥0.25, or ≥0.5, respectively. The correlation between methylation levels on each CpG in clinical samples and various demographical, clinical, or histopathological variables were evaluated by Kendall’s tau statistics. Correlation coefficients and p-values were plotted in the scatter plots. A Kendall’s tau p < 0.05 was considered significant.

77First trimester placentas and maternal blood samples were collected from healthy Caucasian women undergoing termination of pregnancy (n = 12), and processed at the Maternity Clinic and Semmelweis University in Budapest, Hungary. Villous tissues were dissected from the choriodecidua on dry ice and stored at −80oC. Aliquots of maternal sera and plasma were stored at −80oC.

78Pregnancies were dated according to ultrasound scans between 5 and 13 weeks of gestation. Patients with multiple pregnancies were excluded. The collection and investigation of human clinical samples were approved by the Health Science Board of Hungary. Written informed consent was obtained from women prior to sample collection, and the experiments conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report. Specimens and data were stored anonymously.

79Total RNA was isolated from snap-frozen first trimester placental villous tissues (n = 12) with Direct-zol RNA MiniPrep Kit (Zymo Research Corporation) according to the manufacturer’s recommendations. The 28S/18S ratios and the RNA integrity numbers were assessed using an Agilent Bioanalyzer 2100 (Agilent Technologies), and RNA concentrations were measured with NanoDrop 1000. Total RNA (500 ng) was reverse transcribed with the qScript cDNA Synthesis Kit (Quanta Biosciences, Gaithersburg, MD, USA). TaqMan Assays (Applied Biosystems; Data S17 in Supplementary Material) were used for expression profiling on the StepOnePlus Real-Time PCR System (Applied Biosystems).

80Concentrations of leptin and human placental lactogen (hPL) in first (n = 12) and third (n = 19) trimester maternal blood samples were measured with sensitive and specific immunoassays (Leptin ELISA Kit, Abnova, Taipei City, Taiwan; Human Placental Lactogen ELISA Kit, Alpco, Salem, NH, USA) according to the manufacturers’ instructions. Standard curves were generated, and sample assay values were extrapolated. The sensitivities of the assays were <4.4 ng/ml (leptin) and <550 ng/ml (hPL).

81Placental gene expression was measured with either microarray on third trimester samples or qRT-PCR on first trimester samples as described above. Maternal plasma protein concentrations were measured with the above-described immunoassays on respective blood samples taken from the same patients on the day of either delivery or termination of pregnancy. Correlations between placental gene expression and maternal plasma concentrations of gene product proteins were calculated with the Pearson method and visualized on scatter plots (Figure S4 in Supplementary Material).

82To detect disease-associated protein signatures of placental dysfunction in maternal blood, similar to plasma DNA tissue mapping for noninvasive prenatal, cancer, and transplantation assessments, we performed “virtual” liquid biopsy of the placenta. Briefly, human microarray data on 79 human tissues and cells were downloaded from the BioGPS database, which was used for the generation of placenta enrichment scores (placental expression/mean expression in 78 other tissues and cells). Five genes (CGB3, CSH1, ENG, FLT1, and LEP) with enrichment scores between 1.4 and 1,490 were selected based on a literature search due to the extensive investigations of their gene products in maternal blood in preeclampsia (Figure S5 in Supplementary Material). Next, an extensive PubMed search was conducted for first trimester maternal blood protein measurements of these five gene products in patients who developed preeclampsia later in pregnancy. Altogether, 61 scientific reports were identified that met the inclusion criteria, which contained data for 80,170 measurements. These reports were used to build a database for the “virtual” liquid biopsy. In this database, biomarker levels in preterm preeclampsia were expressed as the percentage of control levels as a function of gestational age. Then, the correlations of control percentage values with gestational age were evaluated using the Pearson method. Scatterplots were used to visualize data (Figure S5 in Supplementary Material; Figure 5).

83Women were enrolled in a prospective, longitudinal, and multicenter study in prenatal community clinics of the Maccabi Healthcare Services, Israel. Pregnancies were dated according to the last menstrual period and verified by first trimester ultrasound. Patients with multiple pregnancies or fetuses having congenital or chromosomal abnormalities were excluded. The collection and investigation of human clinical samples were approved by the Maccabi Institutional Review Board. Written informed consent was obtained from women prior to sample collection, and the experiments conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report. Specimens and data were stored anonymously.

84Preeclampsia was defined according to the criteria set by the International Society for the Study of Hypertension in Pregnancy and was subdivided into preterm (<37 weeks) or term (≥37 weeks) groups. Severe preeclampsia was defined by Sibai et al.. SGA was defined as neonatal BW below the 10th percentile for gestational age. Healthy controls had no medical or obstetric complications and delivered a neonate with a BW appropriate for gestational age.

85Peripheral blood samples were obtained by venipuncture in the first trimester from women who subsequently developed preterm severe preeclampsia (n = 5) and term severe preeclampsia (n = 5) as well as healthy controls (n = 10) matched for gestational age at blood draw (Table 4). Blood samples were allowed to clot and were then centrifuged at 10,000 × g for 10 min to separate and collect sera. Serum samples were stored in aliquots at −20oC in the Maccabi Central Laboratory until shipped on dry ice to Hungary.

86Sera were immunodepleted at Biosystems International Ltd., (Debrecen, Hungary) for 14 highly abundant serum proteins on an Agilent 1100 HPLC system using Multiple Affinity Removal LC Column-Human 14 (Agilent Technologies) according to the manufacturer’s protocol. To improve the resolution of 2D gels, immunodepleted serum samples were lyophilized, delipidated, and salt-depleted at the Proteomics Laboratory of the Eotvos Lorand University (Budapest, Hungary). The delipidated and salt-depleted plasma protein samples were dissolved in lysis buffer, and their pH was adjusted to 8.0.

87Protein concentrations of the immunodepleted, desalted, and delipidated serum samples were between 2 and 4 µg/µl as determined with PlusOne 2D Quant Kit (GE Healthcare, Little Chalfont, United Kingdom). Samples were equalized for protein content, and then 5 µg of each protein sample was labeled with a CyDye DIGE Fluor Labeling kit for Scarce Samples (saturation dye) (GE Healthcare) according to the manufacturer’s instructions. Individual samples from cases (n = 2 × 5) and controls (n = 2 × 5) were labeled with Cy5. An internal standard reference sample was pooled from equal amounts (2.5 µg) of all individual samples in this experimental set and was labeled with Cy3. Then, 5 µg of each Cy5-labeled individual sample was merged with 5 µg of the Cy3-labeled reference sample, and these 20 mixtures were run in 2 × 10 gels simultaneously. Briefly, labeled proteins were dissolved in isoelectric focusing (IEF) buffer and were rehydrated passively onto 24 cm immobilized non-linear pH gradient (IPG) strips (pH 3–10, GE Healthcare) for at least 14 h at room temperature. After rehydration, the IPG strips were subjected to the first dimension of IEF for 24 h to attain a total of 80 kVh. Focused proteins were reduced by equilibrating with a buffer containing 1% mercaptoethanol for 20 min. After reduction, IPG strips were loaded onto 10% polyacrylamide gels (24 × 20 cm) and SDS-PAGE was conducted at 12W/gel in the second dimension. Then, gels were scanned in a Typhoon TRIO + scanner (GE Healthcare) using appropriate lasers and filters with the photomultiplier tube biased at 600V. Images in different channels were overlaid using selected colors, and the differences were visualized using Image Quant software (GE Healthcare). Differential protein expression analysis was performed using the differential in-gel analysis and biological variance analysis (BVA) modules of the DeCyder 6.0 software package (GE Healthcare).

88The internal standard reference sample representative of every protein present in all experiments was loaded equally in all gels and thus provided an average image for the normalization of individual samples. The determination of the relative abundance of the fluorescent signal between internal standards across all gels provided standardization between the gels, removing experimental variations and reducing gel-to-gel variations. According to the standard proteomic protocol, the threshold for differential expression was set at 1.05-fold minimum fold-change. A p-value was determined for each protein spot using the Student’s t-test by the BVA module of the DeCyder software. A p-value of <0.05 was considered statistically significant.

89The density of spots in the case of Colloidal Coomassie Blue labeling depends only on the concentration of protein in the sample; however, the density of the spots in the case of saturation dye labeling also depends on the number of cysteines of the labeled proteins because the saturation dye labeling method labels all available cysteines on each protein. This results in the same pattern with different density among samples on the analytical and the preparative gels, rendering identification more difficult. To eliminate this problem for the exact identification of proteins in spots of interest, preparative 2D gel electrophoresis was performed using CyDye saturation fluorescent labeling and Colloidal Coomassie Blue labeling in the same gel. A total of 800 µg of proteins per each of the two gels was run. Following electrophoresis, gels were scanned in a Typhoon TRIO + scanner as described above, the significantly altered spots were matched among the “master” analytical and the fluorescent preparative gel images using the BVA module of the DeCyder 6.0 software package. The resolved protein spots were visualized by the Colloidal Coomassie Blue G-250 staining protocol. Individual spots of interest were excised from the gels based on the comparison of the matched images.

90The excised protein spots were analyzed at the Proteomics Research Group of the Biological Research Center of the Hungarian Academy of Sciences (Szeged, Hungary); the detailed protocol is described in http://ms-facility.ucsf.edu/ingel.html. Briefly, salts, SDS, and Coomassie brilliant blue dye were washed out, disulfide bridges were reduced with dithiothreitol, and then free sulfhydryls were alkylated with iodoacetamide. Digestion with side-chain protected porcine trypsin (Promega, Madison, WI, USA) proceeded at 37°C for 4 h. The resulting peptides were extracted from the gel using 1% formic acid in 50% acetonitrile; then the samples were dried down and dissolved in 0.1% formic acid.

91Samples were analyzed on an Eldex nanoHPLC system online coupled to a 3D ion trap tandem mass spectrometer (LCQ Fleet, Thermo Scientific) in “triple play” data-dependent acquisition mode, where MS acquisitions were followed by CID analyses on computer-selected multiply charged ions. HPLC conditions included in-line trapping onto a nanoACQUITY UPLC trapping column (Symmetry, C18 5 µm, 180 µm × 20 mm; and 15 µl/min with 3% solvent B) followed by a linear gradient of solvent B (5–60% in 35 min, flow rate: 300 nl/min, where solvent A was 0.1% formic acid in water and solvent B was 0.1% formic acid in acetonitrile) on a Waters Atlantis C18 Column (3 µm, 75 µm × 100 mm).

92Raw data files were converted into Mascot generic files (*.mgf) with the Mascot Distiller software v2.1.1.0 (Matrix Science Inc, London, UK). The resulting peak lists were searched against a human subdatabase of the non-redundant protein database of the NCBI (NCBInr, Bethesda, MD, USA) in MS/MS ion search mode on an in-house Mascot server v2.2.04 using Mascot Daemon software v2.2.2 (Matrix Science Inc). Monoisotopic masses with peptide mass tolerance of ±0.6 Da and fragment mass tolerance of ±1 Da were submitted. Trypsin with up to two missed cleavages was specified as an enzyme. Carbamidomethylation of cysteines was set as fixed modification, and acetylation of protein N-termini, methionine oxidation, and pyroglutamic acid formation from peptide N-terminal glutamine residues were permitted as variable modifications. Acceptance criterion was set to at least two significant (peptide score >40, p < 0.05) individual peptides per protein. Localization of identified peptides in the core protein sequences was visually analyzed in order to identify potential protein split products.

93Biological functions of the altered serum proteins were retrieved from Pathway Studio 9.0 software (Ariadne Genomics Inc., Rockville, MD, USA), and from the open-access GO database7 (Figure 6; Figure S6 and Data S8 in Supplementary Material). To elucidate possible interactions between the altered serum proteins and placental genes in the microarray data, bioinformatics analysis was performed using the same software. Molecular networks between the changed serum proteins in preterm (n = 19) or term preeclampsia (n = 14) and DE placental genes annotated in the GO database (n = 1,142) were built separately with a non-linear literature processing search engine, and the resulting connections were manually validated by reading full-text publications (Figure 6; Figure S6 in Supplementary Material). The Fisher’s exact test was used to test for the enrichment of the connections between the altered serum proteins and (1) DE genes in individual modules, taking the connections between the proteins and DE genes in all modules as a background; and (2) DE placental genes, taking the connections between the proteins and all genes tested on the array as a background. To reveal the pathways enriched among the DE genes connected to angiotensinogen, the Ingenuity Pathway Analysis software (QIAGEN, Redwood City, Redwood City, CA, USA) was used, which utilizes Fisher’s exact test and Benjamini–Hochberg correction for multiple testing for the analyses. Statistical significance was set at q < 0.01 (Data S10 in Supplementary Material).

94Caucasian women were enrolled in a prospective study at the Department of Obstetrics and Gynaecology of the University of Debrecen and at the Andras Josa Teaching Hospital in Nyiregyhaza, Hungary. Pregnancies were dated according to ultrasound scans between 8 and 14 weeks of gestation. Patients with multiple pregnancies or fetuses having congenital or chromosomal abnormalities were excluded. The collection and investigation of human clinical samples were approved by the Regional Ethics Committee of the University of Debrecen. Written informed consent was obtained from women prior to sample collection, and the experiments conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report. Specimens and data were stored anonymously.

95Women were included in the following groups: (1) preterm severe preeclampsia with SGA (n = 5) and (2) controls (n = 10). Women were matched for gestational age at blood draw (Table 5). Preeclampsia was defined according to the criteria set by the American College of Obstetricians and Gynecologists and was subdivided into preterm (<37 weeks) or term (≥37 weeks) groups. Severe preeclampsia was defined according to Sibai et al.. SGA was defined as neonatal BW below the 10th percentile for gestational age. Healthy controls had no medical or obstetric complications and delivered a neonate with a BW appropriate for gestational age.

96Peripheral blood samples were obtained by venipuncture in the first trimester. Plasma samples were separated by double centrifugation for 10 min. Samples were stored in aliquots at −80oC.

97All solvents were HPLC-grade from Sigma-Aldrich (St. Louis, MO, USA) and all chemicals, where not stated otherwise, were obtained from Sigma-Aldrich. Frozen plasma samples were thawed and denatured with denaturing buffer (Biognosys AG; Schlieren, Switzerland). Samples were alkylated using alkylation solution (Biognosys). Subsequently, samples were digested overnight with sequencing grade modified trypsin (Promega; Madison, WI, USA) at a protein:protease ratio of 50:1. C18 cleanup for mass spectrometry was carried out according to the manufacturer’s instructions using C18 Micro Spin columns (Nest Group Inc.; Southborough, MA, USA). Peptides were dried down to complete dryness using a SpeedVac system. Dried peptides were redissolved with LC solvent A (1% acetonitrile in water with 0.1% formic acid). Final peptide concentrations were determined for all samples by 280 nm measurement (SpectrostarNANO, BMG Labtech, Offenburg, Germany). Samples were spiked with PlasmaDive™ (Biognosys) reference peptides mix at known concentrations.

98Peptides (1 µg per sample, corresponding to injection of 0.0259 µl of initial plasma sample) were injected to a self-packed C18column [75 µm inner diameter and 10 cm column length, New Objective (Woburn, MA, USA); column material was Magic AQ, 3 µm particle size, and 200Å pore size from Michrom] on a ThermoScientific EASY-nLC1000 nano-liquid chromatography system. LC–MRM assays were measured on a ThermoScientific TSQ Vantage triple quadrupole mass spectrometer equipped with a standard nano-electrospray source. The LC gradient for LC–MRM was 6–40% solvent B (85% acetonitrile in water with 0.1% formic acid) for 30 min followed by 40–94% solvent B for 2 min and 94% solvent B for 8 min (total gradient length was 40 min). A subset of Biognosys’ PlasmaDive™ MRM Panel of 10 peptides representing 10 proteins (Data S9 in Supplementary Material) was used for the measurements of 10 altered proteins involved in immune responses. For the quantification of the peptides across samples, the TSQ Vantage was operated in scheduled MRM mode with an acquisition window length of 5 min. The LC eluent was electrosprayed at 1.9 kV and Q1/Q3 were operated at unit resolution (0.7 Da). Signal processing and data analysis were carried out using SpectroDive™ 8.0—Biognosys’ software for multiplexed MRM/PRM data analysis based on mProphet. A q-value filter of 1% was applied.

99Because of the small sample size of this pilot study, statistical simulations were carried out to predict power and significance in the extended validation study by multiple permutation testing. The probability of significant differences between the two groups was estimated by a paired t-test and the Mann–Whitney test. For both, p < 0.05 was taken as criteria to count successful simulations. The number of successful simulations was calculated for different sample sizes (n = 10 or 100) and repeats of simulations (n = 10, 100, or 1,000). The simulated significance level at p < 0.05 was accepted if the number of successful simulations was >25% (Figure 6; Data S9 in Supplementary Material).

100An extensive PubMed search was conducted for first trimester maternal blood protein measurements of the DE proteins found by 2D-DIGE. Altogether five scientific reports were identified that met the inclusion criteria. Figure 6 depicts biomarkers with the same direction differential abundance in preterm preeclampsia in published data as in 2D-DIGE assays.

101To build optimal in vitro cellular models of trophoblastic disease, an extensive PubMed search was first conducted for similar assays. Since no human trophoblastic stems cells were yet available that would enable a natural proliferative trophoblastic pool with differentiation potential into villous or extravillous lineages, choriocarcinoma-derived trophoblastic cell lines or immortalized EVTs were mostly used for such purposes. Our search revealed BeWo cells and HTR8/SVneo cells as the increasingly most accepted cell model systems based on 1,500 published articles.

102For in vitro trophoblast experiments, placentas (n = 6) were collected prospectively at the Perinatology Research Branch (NICHD, NIH, DHHS) from normal pregnant women at term who delivered an AGA neonate with Cesarean section. Cytotrophoblasts were isolated from these placentas by the modified method of Kliman et al.. Briefly, 100 g villous tissues were cut, rinsed in PBS, and sequentially digested with Trypsin (0.25%; Life Technologies, Grand Island, NY, USA) and DNAse I (60 U/ml; Sigma-Aldrich) for 90 min at 37oC. Dispersed cells were filtered through 100 µm Falcon nylon mesh cell strainers (BD Biosciences, San Jose, CA, USA), and then erythrocytes were lysed with 5 ml NH4Cl solution (Stemcell Technologies, Vancouver, BC, Canada). Washed and resuspended cells were layered over 20–50% Percoll gradients and centrifuged for 20 min at 1,200 × g. Trophoblast containing bands were collected and non-trophoblastic cells were excluded by negative selection using anti-CD9 (20 µg/ml) and anti-CD14 (20 µg/ml) mouse monoclonal antibodies (R&D Systems, Minneapolis, MN, USA), MACS anti-mouse IgG microbeads, and MS columns (Miltenyi Biotec, Auburn, CA, USA). Then, primary VTs were plated on a collagen-coated 12-well plate (BD Biosciences; 3 × 106 cells/well) in Iscove’s modified Dulbecco’s medium (IMDM; Life Technologies) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (P/S). To test the effect of trophoblast differentiation on selected genes’ expression, primary trophoblasts were kept in IMDM containing 5% non-pregnant human serum (SeraCare, Milford, MA, USA) and 1% P/S. The medium was replenished every 24 h, and cells were harvested for total RNA every 24 h between days 1 and 7.

103Villous trophoblasts were isolated from normal term placentas and plated as described above. To test the effect of preeclampsia serum on VTs at two time points of differentiation, VTs were kept in IMDM containing 1% P/S and 10% first trimester maternal sera from control or preeclamptic women. Medium was replenished after 48 h, and cells were harvested for total RNA either 24 or 72 h after the start of the human serum treatment. All experiments were run in triplicate.

104BeWo cells (American Type Culture Collection) were incubated in a T-25 flask or 6-well plate with Ham’s F12-K medium (Life Technologies) supplemented with 10% FBS and 1% P/S in a humidified incubator (5% CO2, 20% O2) at 37oC until reaching 50–80% confluence. To test the effect of ARNT2 or BCL6 overexpression on gene expression, cells were transiently transfected with ARNT2, BCL6, or control (GFP) vectors. Briefly, 4 µg expression plasmid (OriGene Technologies, Inc., Rockville, MD, USA) and 12 µl FuGENE HD transfection reagent (Promega) were mixed with 180 µl Ham’s F12-K medium (10% FBS, 1% P/S), incubated at RT for 15 min and added to cell cultures with 1.8 ml medium in each well of 6-well plate. Twenty-four hours after transfection, cells were split into three treatment groups and kept under either normoxic (20% O2), hypoxic (2% O2), or ischemic (1% and 20% O2 alternating for 6 h) conditions in an Oxycycler C42 (BioSpherix, Lacona, NY, USA) for 48 h before cell harvest. This study setup followed the generally accepted view in reproductive sciences that ~2% O2 concentration represents physiologic hypoxia at the implantation site and that placental development occurs under physiologic hypoxia in the first trimester.

105Total RNA was isolated from primary VTs on days 0–7 of differentiation and from BeWo cell cultures with TRIzol reagent (Life Technologies) and RNeasy kit (QIAGEN) according to the manufacturers’ recommendations. The 28S/18S ratios and the RNA integrity numbers were assessed using an Agilent Bioanalyzer 2100; RNA concentrations were measured with NanoDrop 1000. Total RNA (500 ng) was reverse transcribed with High Capacity cDNA Reverse Transcription Kit using random hexamers (Applied Biosystems). TaqMan Assays (Applied Biosystems; Data S17 in Supplementary Material) were used for high-throughput gene expression profiling on the Biomark qRT-PCR system (Fluidigm) according to the manufacturers’ instructions.

106Data were analyzed using the ΔΔCt method. Data were first normalized to the reference gene (RPLP0), and log2 mRNA relative concentrations were obtained for each sample as −ΔCt(gene) = Ct(RPLPO) − Ct(gene).

107The overall changes in gene expression during the 7 days of differentiation were analyzed by comparing the mean expressions on a given day versus Day 0. The highest fold change for a given gene was defined as the maximum of the daily expression differences during the 7-day time-period. Significant differences were defined by a paired t-test (p < 0.05) (Figure 7).

108Gene expression data were analyzed using the Student’s t-test to compare the effect of preeclampsia serum with the effect of control serum on gene expression at Days 1 and 3 of trophoblast differentiation. P-values of <0.05 were considered significant (Figure 6).

109Gene expression data were analyzed to compare the effect of ARNT2 or BCL6 overexpression with the effect of control vector overexpression on gene expression in normoxic conditions using a one-way ANOVA model. The same model was used to assess the differential effect of ARNT2, BCL6, or GFP overexpression on gene expression in hypoxic or ischemic conditions vs. normoxia. P-values of <0.05 were considered significant (Figure 8).

110A permutation test was used to measure the statistical significance of the matching between differential gene expression patterns in in vitro and in vivo conditions. Genes were discretized into three states, i.e. up-regulated (UP), down-regulated (DN), or unchanged (NS). For each gene in the two conditions, a score of 1 was assigned for a perfect match of UP/UP or DN/DN, 0 for a neutral match of NS/NS, −1 for a perfect mismatch of UP/DN or DN/UP, and −0.5 for all other patterns. The matching score for any pair of conditions was computed as the sum of all scores for each individual gene. The significance of the scores was assessed via a permutation on the class labels. Permutations were exhaustive when feasible, otherwise limited to a random sample of 5,000 (Data S11 in Supplementary Material).

111TaqMan assays for ZNF554 and RPLP0 (Data S17 in Supplementary Material) were run in triplicate for expression profiling of the Human Major Tissue qPCR Array (OriGene Technologies) that contains cDNAs from 48 different pooled tissues.

112In situ hybridization on third trimester FFPE placental tissues (n = 6) retrieved from the Bank of Biological Specimens of the Perinatology Research Branch was carried out using the RNAscope 2.0 FFPE Assay–Brown (Advanced Cell Diagnostics, Hayward, CA, USA) on a HybEZ Hybridization System (Figure 9). Briefly, tissue sections were incubated with ZNF554 target probe (Cat.No.: 423831, Advanced Cell Diagnostics) for 2 h at 40°C. After rinsing with 1× Wash Buffer, slides underwent a six-step amplification procedure at 40°C and were washed with 1× Wash Buffer between amplification steps. Chromogenic detection was performed using a 1:1 mixture of Brown-A and Brown-B solutions. Slides were counterstained with hematoxylin, dehydrated in graded ethanol, and mounted in xylene.

113Third trimester placentas were retrieved from the Bank of Biological Specimens of the Perinatology Research Branch. First trimester placentas were collected from healthy Caucasian women undergoing termination of pregnancy and processed at the Maternity Clinic and Semmelweis University in Budapest, Hungary as described above. Five µm sections of first and third trimester FFPE placental tissues (n = 15) were placed on silanized slides and stained using anti-ZNF554 or anti-cytokeratin-7 antibodies as well as reagents listed in Data S18 (Supplementary Material) either on Ventana Discovery (Ventana Medical Systems, Inc, Tucson, AZ, USA) or Leica BOND-MAX (Leica Microsystems, Wetzlar, Germany) autostainers (Figures 9 and 10).

114The experimental procedures were embedded in the study on VT differentiation as described above (Figure 9).

115BeWo cells were incubated in a T-25 flask or 6-well plate with Ham’s F12-K medium (Life Technologies) supplemented with 10% FBS and 1% P/S in a humidified incubator (5% CO2, 20% O2) at 37oC until reaching 50–80% confluence.

116To test the effect of ZNF554 knockdown on gene expression, cells were treated either with 100 nM ZNF554 siRNA (Ambion-Life Technologies, Foster City, CA, USA) or 100 nM scrambled (control) siRNA (Ambion) using X-tremeGENE siRNA transfection reagent (Roche, Mississauga, ON, Canada), and incubated at 37oC in 2 ml serum-free Opti-MEM (Gibco-Life Technologies) medium. After 6 h, the medium was replaced with 2 ml Ham’s F12-K medium supplemented with 10% FBS. After 48 h, cells were collected for RNA isolation, microarray, and qRT-PCR as well as confocal microscopy (Figure 9).

117To test the effect of DNA methylation on gene expression, BeWo cells were treated with 5 or 10 µM 5-azacitidine (Sigma-Aldrich), and control cells with DMSO. This experiment was also performed when both 5-azacitidine-treated and control cells received 25 µM forskolin (Sigma-Aldrich) to induce syncytialization. After 24 h incubation, cells were harvested for RNA isolation and qRT-PCR. The experiment was performed in six replicates (Figure 11; Figure S8 in Supplementary Material).

118HTR8/SVneo EVT cells (kindly provided by Dr. Charles H. Graham, Queen’s University, Kingston, Ontario, Canada) were incubated in a 6-well plate with RPMI-1640 medium (Gibco-Life Technologies) supplemented with 10% FBS and 1% P/S in a humidified incubator (5% CO2, 20% O2) at 37oC until reaching 50% confluence.

119To test the effect of ZNF554 knockdown on gene expression and functions, cells were treated either with 100 nM ZNF554 siRNA or 100 nM scrambled (control) siRNA as described for BeWo cells. After 6 h, the medium was replaced with 2 ml RPMI-1640 medium (Gibco-Life Technologies) supplemented with 10% FBS. On the following day, cells were kept in various O2 concentrations (2, 8, or 20%) in an Oxycycler C42. Cells were collected for functional assays after 24 h, while cells were collected for RNA isolation, microarray, qRT-PCR, or confocal microscopy and their supernatants for ELISA after 48 h. Cell cultures were used for cell proliferation assays after 0, 24, and 48 h.

120Total RNA was isolated from BeWo and HTR8/SVneo cell cultures with TRIzol reagent (Life Technologies) and RNeasy kit (QIAGEN, Valencia, CA, USA) according to the manufacturers’ recommendations. The 28S/18S ratios and the RNA integrity numbers were assessed using an Agilent Bioanalyzer 2100; RNA concentrations were measured with NanoDrop 1000. DNase-treated RNA from BeWo and HTR8/SVneo cells (500 ng) was amplified and biotin-labeled with the Illumina TotalPrep RNA Amplification Kit (Ambion-Life Technologies). Labeled cRNAs were hybridized to a HumanHT-12v4 Expression BeadChip (Illumina, Inc., San Diego, CA). BeadChips were imaged using a BeadArray Reader (Illumina, Inc.), and raw data were obtained with BeadStudio Software V.3.4.0 (Illumina, Inc.). Total RNA (500 ng) was also reverse transcribed with High Capacity cDNA Reverse Transcription Kit using random hexamers (Applied Biosystems, Foster City, CA, USA). TaqMan Assays (Applied Biosystems; Data S17 in Supplementary Material) were used for high-throughput gene expression profiling on the Biomark qRT-PCR system (Fluidigm, San Francisco, CA, USA) according to the manufacturers’ instructions.

121ZNF554 knock-down and control BeWo and HTR8/SVneo cells cultured in 6-well plate were detached with 0.05% Trypsin-EDTA (Life Technologies), washed, and resuspended in PBS, and then 4 × 104 cells were cytospined to Superfrost Plus slides (Fisher Scientific). Then, cells were fixed with 4% paraformaldehyde (Electron Microscopy Sciences, Hatfield, PA, USA), blocked with Protein Block (Dako North America, Inc., Carpinteria, CA, USA), and immunostained with anti-ZNF554 mouse polyclonal antibody (1:100 dilution, overnight; Abnova, Taipei City, Taiwan) and an AlexaFluor-488 goat anti-mouse antibody (1:1,000 dilution; Life Technologies). Cells were mounted with ProLong Gold antifade reagent with 4’,6-diamidino-2-phenylindole (DAPI; Life Technologies), followed by confocal microscopy using a Leica TCS SP5 MP spectral confocal system (Leica Microsystems) (Figures 9 and 10).

122Concentrations of human plasminogen activator inhibitor-1 (PAI-1) and TIMP-3 in HTR8/SVneo cell culture supernatants were measured with sensitive and specific immunoassays (Human PAI-1 ELISA Kit, Invitrogen; TIMP-3 Human ELISA Kit, Abcam Inc., Cambridge, MA, USA) according to the manufacturers’ instructions. Standard curves were generated, and sample assay values were extrapolated. The sensitivities of the assays were <30 pg/mL (PAI-1) and <2 pg/ml (TIMP-3) (Figure 10).

123Cell cultures of control and ZNF554 siRNA treated HTR8/SVneo cells were assayed using the CellTiter 96 Aqueous Non-radioactive Cell Proliferation Assay (Promega) according to the manufacturer’s instructions (Figure S11 in Supplementary Material).

124The migratory capacity of HTR8/SVneo cells was examined with 8μm-pore transwell inserts (Corning, NY, USA) inserted in a 12-well plate described previously. After transfection with ZNF554 or scrambled siRNAs for 24 h, 5 × 105 HTR8/SVneo cells were plated in the upper chambers in a serum-free RPMI-1640 medium, whereas the lower chambers contained an RPMI-1640 medium supplemented with 10% FBS. After incubation for 36 h in 2, 8, or 20% O2 concentrations, cells on the upper side of the membranes were removed by cotton swab, and the inserts were fixed in methanol for 10 min at RT and washed once with PBS. Then, the membranes were cut out and mounted on Superfrost Plus slides (Fisher Scientific) with ProLong Gold antifade reagent with DAPI. Comprehensive images of each membrane were taken using a Leica TCS SP5 MP spectral confocal system. The number of invaded cells was quantified using Image-Pro Premier v9.0.2 (Media Cybernetics, Inc., Rockville, MD, USA). The experiment was performed in six replicates (Figure 10).

125The invasiveness of HTR8/SVneo cells was examined with a Matrigel invasion assay using 8 μm-pore cell culture inserts (BD Biosciences) pre-coated with Matrigel (125 ng/ml; BD Biosciences) and inserted in a 24-well plate described previously. After transfection with ZNF554 or scrambled siRNAs for 24 h, 2 × 105 cells were plated in the upper chambers in a serum-free RPMI-1640 medium, whereas an RPMI-1640 medium supplemented with 10% FBS was added to the lower chambers. After incubation for 48 h in 2, 8, or 20% O2 concentrations, cells on the Matrigel side of the membranes were removed by cotton swab, and the membranes were processed as in the migration assay. Comprehensive images taken using the Leica TC5 SP5 spectral confocal system were quantified using Image-Pro Plus 6.2 (Media Cybernetics, Inc.). The experiment was performed in triplicate (Figure 10).

126The expression of ZNF554 relative to RPLP0 in the placenta was compared to 47 other human tissues using the Student’s t-test. P-values of <0.05 were considered significant (Figure 9).

127The Student’s t-test was used to evaluate ZNF554 knock-down efficiency and the effect of ZNF554 knock-down on gene expression in BeWo and HTR8/SVneo cells (Figures 9 and 10).

128Data were analyzed using the Bioconductor packages in R following the MIAME guidelines and methodologies described previously. Raw microarray gene expression data were normalized by a quantile normalization approach. A moderated t-test was used to select DE genes using a cutoff of >1.5 fold-change and q < 0.1. GO analysis and pathway analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database were also performed.

129qRT-PCR data were analyzed using the ΔΔCt method relative to RPLP0 expression. The Student’s t-test was used to evaluate ZNF554 knock-down efficiency in HTR8/SVneo cells and the effect of ZNF554 knock-down on gene expression and cell proliferation. A linear model was built to quantify the effects of ZNF554 knock-down and various O2 concentrations on the gene expression and protein secretion of HTR8/SVneo cells as well as on their migratory and invasive capacity. O2 concentration was treated as a continuous variable, and the interaction between ZNF554 knock-down and O2 concentrations was retained in the model when the coefficient was significant. P-values of <0.05 were considered significant (Figure 10).

130All relevant data are within the paper and its Supplementary Material files. MIAME compliant microarray data are available from the Gene Expression Omnibus (GEO) under accession numbers GSE65866, GSE65940, and GSE66273.

131The collection and use of human biological specimens and clinical data for research purposes were approved by the Health Science Board of Hungary, the Institutional Review Board of the Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD, NIH, DHHS), the Wayne State University Human Investigation Committee, the Maccabi Institutional Review Board, and the Regional Ethics Committee of the University of Debrecen. Written informed consent was obtained from women prior to sample collection, and the experiments conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report. Specimens and data were stored anonymously.

132NGT conceptualized study and designed research. NGT, KAK, YX, KJ, RJL, EH-G, ZsD, AS, KE, SzSz, VT, HE-A, CL, AB, GSz, SL, and ZD performed research. NGT, RR, ALT, ZX, LO, OT, HM, SD, SSH, THC, CJK, and ZP contributed new reagents/analytic tools/clinical specimens. NGT, RR, ALT, KAK, YX, ZX, KJ, GB, ZsG, JP, THC, BAGy, AD, ASz, ZsD, GSz, IK, AF, MKr, MKn, OE, GJB, CJK, GJ, and ZP analyzed and interpreted data. All authors contributed to manuscript writing and approved the paper.

133Part of the data was submitted by NGT, RR, ALT, KAK, RJL, THC, HM, MK, GJ, and ZP as patent applications in 2012 and 2014 to describe biomarkers for preeclampsia. Other authors declare no conflict of interest.