Preeclampsia is a hypertensive disorder of pregnancy. It is associated with a higher risk of hypertension and cardiovascular disease later in life. Women who had preeclampsia have a twofold increased risk of death from cardiovascular diseases (Bokslag et al., 2016; Neerukonda et al., 2020). There is evidence that preeclampsia originates in part from genetic causes that include contributions from the maternal, paternal and fetal genome (Zusterzeel et al., 2002; Cnattingius et al., 2004; Nilsson et al., 2004; Kobayashi, 2015; Than et al., 2018). The heritability of preeclampsia is up to 52% (Salonen Ros et al., 2000; Chappell and Morgan, 2006). The role of genetics in preeclampsia is also supported by family-based observations (Chappell and Morgan, 2006; Nejatizadeh et al., 2008) with more than 100 studies showing a 2- to 5-fold increased risk among family members of affected women (Chesley et al., 1968; Sutherland et al., 1981; Arngrimsson et al., 1990; Cincotta and Brennecke, 1998; Mutze et al., 2008; Ward, 2008). The recurrence risk for preeclampsia in the daughters of either eclamptic or preeclamptic patients is 20–40% (Serrano, 2006; Genc and Schantz-Dunn, 2007). Nonetheless, there is no current consensus among the published results in regards to associated genes and the pathogenesis of the disease.

The genetic basis for complex diseases involves the interaction of multiple genes in discrete networks and pathways (Loscalzo et al., 2007). Although complex diseases show increased recurrence risk in families, they do not follow a simple Mendelian pattern of inheritance (Smith and Ebrahim, 2003). Computational methods have been used to analyze networks of genes and to find biological subnetworks due to the genetic heterogeneity of e diseases that are linked to a variety of disorders (Wall et al., 2009). There are several studies employing computational methods to identify important genes associated with hypertension. Ran et al. analyzed the protein protein interaction (PPI) network topology and molecular connectivity between protein pathways to identify associations with hypertension (Ran et al., 2013). Researchers have also developed a machine-learning algorithm to predict novel hypertension associated genes (Li et al., 2017).

We hypothesize that the genetic architecture of complex diseases like preeclampsia is described by clusters of patients with variants in genes in shared PPI networks. We sought to test this hypothesis using whole exome sequencing in carefully selected patients with early onset preeclampsia with severe clinical features. We compared variants identified in women with early onset, idiopathic severe preeclampsia with term controls without personal or family history of pregnancy related hypertensive disorders. We built and implemented Proteinarium, a novel multi-sample PPI tool, to identify clusters of patients with shared PPI networks.

The clinical characteristics and the race/ethnicity distribution of the patients are shown in Table 1. As can be seen, the gestational age at delivery, systolic blood pressure, frequency of proteinuria, impaired liver function, thrombocytopenia, cerebral visual symptoms and fetal growth retardation were all significantly different between the groups, which was expected by our definition of severe preeclampsia.

High quality sequence data with a Phred score ≥30 from well-balanced pools with over 19, 000, 000 reads/patient, 40X average depth of coverage, with more than 80% of sequence reads with at least 20X coverage were observed. We identified a total of 528,630 variants including 187,915 exonic variants. The work flow for the univariate analysis is shown in Figure 1. After application of the initial filters for coverage and variant pathogenicity, there were 8,868 predicted deleterious variants (available at Online_Supplementary Table S1 ). Among these, 21 variants were nominally associated with preeclampsia by genotype testing. All were non-synonymous, exonic variants (Table 2). Nonetheless, none of these variants met genome-wide significance after correction for multiple comparison testing.

Proteinarium was used to carry out the multivariate protein protein interaction analysis to identify clusters of patients with shared networks associated with severe preeclampsia. The resulting, circularized dendrogram is shown in Figure 2A (Ciccarelli et al., 2006). Out of the 143 patients sequenced, 129 patients were assigned to two statistically significant clusters. (p < 0.0001). The inset in Figure 2A shows the number of cases and controls in each cluster. Cluster A had significantly more cases than controls, containing 47 of the 61 case subjects (Bastian et al., 2009). The layered network for the case-dominated Cluster A is shown in Figure 2B. There are 13 genes which are unique to Cluster A, highlighted in red in the layered network graph. Most have defined functional roles or implications for preeclampsia, Table 3. Cluster B had significantly more controls than cases, including 61 of the 82 subjects. The layered network for Cluster B is shown in Figure 2C. The unique genes from the layered network graph of Cluster B, shown in blue, are listed in Supplementary Table S1.

We compared the sequence data of the samples not assigned to clusters with those that were assigned and we did not find significant differences in the average depth of coverage. Likewise, there were no significant differences in clinical/phenotypic characteristics when comparing the subjects in the significant clusters with the subjects that were not in these clusters (data not shown).

The comparison of the unique genes from the case and the control dominated clusters revealed a positive separation score, confirming that the layered PPI networks of these two patient subgroups indeed exist in distinct areas of the interactome. We ran GO term analysis using DAVID software on all genes of the network from Cluster A and from Cluster B, Supplementary Table S2 (Huang da et al., 2009b; a). We found significantly enriched biological processes, molecular functions and cellular components based on Bonferroni corrected p-value for the case and control dominated networks. Prominent among the biological processes and molecular functions associated with preeclampsia were antigen processing and presentation, cellular movement (axon guidance and microtubules) and T cell receptor signaling.

We previously reported the database for Preeclampsia (dbPEC) which archives a curation-based collection of genes associated with preeclampsia and their association with clinical features and concurrent conditions (Triche et al., 2014). We compared the genes from our univariate analysis and the genes from both case and control dominated layered networks to those in the database. We found two overlapping genes from the univariate gene list (TTN and CCL14) that were included in dbPEC. We also found three overlapping genes from the layered network of Cluster A (FN1, KIF2A, VCAM1). By over representation analysis, Cluster A is significantly enriched for genes previously shown to be associated with preeclampsia in dbPEC (p < 0.0033).

Preeclampsia is a life-threatening, multi-system hypertensive disorder of pregnancy, which complicates up to 15% of US deliveries (Chappell and Morgan, 2006; Valenzuela et al., 2012; Bornstein, 2020; Steinthorsdottir et al., 2020). The incidence is increasing (Bornstein, 2020). It is recognized as a leading cause of maternal and fetal morbidity and mortality worldwide (Valenzuela et al., 2012). Preeclampsia is characterized by varying degrees of maternal symptoms including elevated blood pressure, proteinuria and fetal growth restriction (Jebbink et al., 2012). Many clinicians believe that preeclampsia, severe preeclampsia, and early vs late preeclampsia are different disorders (Carreiras et al., 2002; Raymond and Peterson, 2011; American College of et al., 2013). Previously, using bioinformatic methods, we showed that there are discrete gene sets associated with these different phenotypes of preeclampsia (Triche et al., 2014).

We performed whole exome sequencing on women with idiopathic early-onset preeclampsia with severe features and singleton births <34 weeks’ gestation and compared them to term controls with no family history of preeclampsia. We developed Proteinarium, a novel multi-sample PPI analysis and visualization tool, to identify clusters of patients with shared protein protein interaction networks (Armanious et al., 2020). Unlike other tools which use the differentially expressed genes or the phenotype associated variants to build a single network for all samples, or clustering models based solely on gene expression, Proteinarium builds individual networks for each patient based on the STRING database. The similarities between individual PPI networks are then evaluated using a distance metric for clustering the samples. We identified a single, significant, large cluster of patients with a predominance of cases with early-onset, severe features of preeclampsia. We also identified a single control-dominated cluster. The separation test of the unique genes from case and control dominated clusters confirmed that the two subnetworks forming clusters A and B exist in the different regions of the interactome. These results support our hypothesis that the genetic architecture of complex diseases is characterized by clusters of patients that have variants in shared gene networks and provide insights into the genetics of severe preeclampsia.

Several of the unique genes from the case dominated network have very plausible mechanistic connections to preeclampsia. Laminin β2 (LAMB2) is a glomerular basement membrane (GBM) component, required for proper functioning of the glomerular filtration barrier. It has a role in proteinuria (Zhang and Huang, 2012) and serum laminin levels in preeclamptic patients are significantly higher than those in normal pregnancy (Furuhashi et al., 1993). Hypoxia-induced upregulation of Quiescin Sulfhydryl oxidase 1 (QSOX1) and an elevation in intracellular H₂O₂ leads to increased apoptosis in the placentae of pregnancies complicated by preeclampsia (Li et al., 2019). QSOX1 protein is found in circulating extracellular vesicles of both preeclampsia and healthy pregnant women (Tan et al., 2014). Fibronectin 1 (FN1) might promote the development of preeclampsia by modulating differentiation of human extravillous trophoblasts, as well as formation of focal adhesions (Brubaker et al., 1992; Auer et al., 2010; Zhao et al., 2017). Vascular Cell Adhesion Molecule 1 (VCAM1) is involved in cellular adhesion and serum concentrations of sVCAM-1 are significantly elevated in both mild and severe preeclampsia (Kim et al., 2004). Invasion of maternal decidua and uterine spiral arteries by extravillous trophoblasts is required for establishment of normal placenta. Human trophoblast migration requires Rac Family Small gtpase 1 (RAC1) and Cell Division Cycle 42 (CDC42) (Grewal et al., 2008). Lower levels were found in preeclampsia samples than in normal term pregnancy samples, and decline significantly in severe preeclampsia (Fan et al., 2016). Protein tyrosine kinase 2 (PTK2) (focal adhesion kinase) is differentially expressed in preeclampsia and reported as among the promising biomarkers for preeclampsia (Sado et al., 2011). In the case-dominated subnetwork we observed Kinesin Family Member 2A (KIF2A) which is upregulated in the preeclamptic placenta (Kobayashi, 2015). Up-regulated genes in the preeclampsia placenta have been shown to be associated with the regulation of diverse cellular processes, including matrix degradation, trophoblast cell invasion, migration and proliferation (Kobayashi, 2015).

There have been several sequencing efforts, including whole genome, whole exome and targeted sequencing, on an array of preeclampsia phenotypes from diverse populations (Johnson et al., 2012; Emmery et al., 2016; Kaartokallio et al., 2016; Thomsen et al., 2017; Gammill et al., 2018; Hansen et al., 2018; Soellner et al., 2018; Glotov et al., 2019; Melton et al., 2019; Steinthorsdottir et al., 2020; Zhang et al., 2020). There is no consensus among the published results in regards to associated genes and variants. Since preeclampsia is a complex, polygenic disease, the lack of a consensus among these univariate comparisons might be expected in these early-stage studies. Among the 20 genes identified in our univariate analysis, only Titin (TTN) was identified in prior studies (Gammill et al., 2018; Zhang et al., 2020). Protein-altering mutations in TTN have been identified in patients with cardiomyopathy and women with preeclampsia are more likely to carry TTN mutations associated with idiopathic cardiomyopathy and peripartum cardiomyopathy (Gammill et al., 2018). Additionally, we found 2 genes, Major Histocompatibility Complex, Class II, DQ Alpha 1 (HLA-DQA1) and Inositol 1,4,5-Trisphosphate Receptor Type 1 (ITPR1) that were reported in previous studies of preeclampsia (Emmery et al., 2016; Hansen et al., 2018). None of these overlapping genes were among the unique genes identified in the shared layered network from the case-dominated cluster. Likewise, no overlapping variants or genes were found in a recent genome-wide association meta-analysis investigating genetic predispositions associated with preeclampsia (Steinthorsdottir et al., 2020).

Our analysis allowed us to identify a cluster of patients with shared PPI networks associated with severe preeclampsia. Within the significant cluster, there were unique imputed genes (RAC1, KIF5B, PTK2, KIF5A, FN1, QSOX1, ARF4, VCAM1, CDC42, KIF2A) that were not among the top 60 seed genes selected by genotype testing. Nonetheless, our approach allowed us to identify these influential genes in the mechanism(s) underlying preeclampsia that would not otherwise have been identified by whole genome univariate variant analysis.

We also examined the unique genes in the network of the control dominated cluster. Proteins in this network are associated with the ubiquitination process. They may serve a role that confers resilience against preeclampsia (Fredrickson and Gardner, 2012; Berryman et al., 2019). Although there are studies showing a relationship with hypertension - ubiquitination process and pregnancy, this still needs further investigation (Fredrickson and Gardner, 2012).

Whole exome sequencing, combined with a novel, multi-sample network analysis, and careful phenotyping contributed to our discovery despite the relatively modest size of our study. Concepts developed from network theory suggest that related diseases involve proteins in similar neighborhoods of the interactome (Menche et al., 2015). Based on these concepts, we hypothesized that the genetic architecture of preeclampsia is described by subgroups of patients with variants in shared genes in specific networks and pathways. We identified a significant subgroup of cases with shared PPI networks associated with severe preeclampsia. We believe that the careful phenotyping resulted in the high percentage of subjects being successfully assigned to significant clusters and the ability to observe distinct separation between the case and control dominated clusters. We acknowledge that genetic variation may not be the sole mechanism for preeclampsia, but rather epigenetic changes or protein conformational disturbances may also play a significant role.

While we were not expecting each patient to appear in a significant cluster and our study included only a modest sample size, we identified a significant subgroup of patients with shared PPI networks associated with severe preeclampsia. In order to leverage the likelihood of genetics discovery, we focused exclusively on women with severe, early-onset preeclampsia. Our analysis was restricted to evaluation of genetic variants in the maternal genome only. Future studies including fetal and/or paternal data will enhance the likelihood of genetic discovery.

Using our unique network analysis, we were able to identify subsets of patients with shared networks, thus confirming our hypothesis about the genetic architecture of preeclampsia. Strict phenotyping of both cases and controls improved the likelihood of identifying these otherwise difficult to find genetic associations. Our network analysis identified genes which were imputed from the interactome and these imputed genes provide insights for severe preeclampsia that may otherwise have not been identified. As such, these are important candidates to include in meta-analyses of genetic associations with preeclampsia. These results provide promise to further our understanding the mechanism underlying complex diseases like preeclampsia.