Identification of underlying mechanisms and hub gene-miRNA networks of the genomic subgroups in preeclampsia development.

Preeclampsia is a hypertensive disorder of pregnancy that can lead to multiorgan complications in the mother and fetus. Our study aims to uncover the underlying mechanisms and hub genes between genomic subgroups of preeclampsia. A total of 180 preeclampsia cases from 4 gene profiles were classified into 3 subgroups. Weighted gene coexpression analysis was performed to uncover the genomic characteristics associated with different clinical features. Functional annotation was executed within the significant modules and hub genes were predicted using Cytoscape software. Subsequently, miRNet analysis was performed to identify potential miRNA-mRNA networks. Three key subgroup-specific modules were identified. Patients in subgroup II were found to develop more severe preeclampsia symptoms. Subgroup II, characterized by classical markers, was considered representative of typical preeclampsia patients. Subgroup I was considered as an early stage of preeclampsia with normal-like gene expression patterns. Moreover, subgroup III was a proinflammatory subgroup, which presented immune-related genomic characteristics. Subsequently, miR-34a-5p and miR-106a-5p were found to be correlated with all 3 significant gene modules. This study revealed the transcriptome classification of preeclampsia cases with unique gene expression patterns. Potential hub genes and miRNAs may facilitate the identification of therapeutic targets for preeclampsia in future.

1. Introduction

Preeclampsia (PE) is characterized by hypertension and multiorgan dysfunction (mainly heart, brain, and kidney) after 20 weeks of gestation. The condition may lead to long-term maternal and neonatal complications.[ The global incidence rate of PE among pregnant women is 5% to 7%, and this pregnancy-specific disorder is responsible for over 70,000 maternal deaths and 500,000 fetal deaths per year worldwide.[ Several risk factors for progression of PE have been identified, such as previous medical history, environmental factors, oxidative stress, genetic factors, and immunologic dyshomeostasis.[ There are limited clinical interventions for PE, such as antihypertensive drugs and magnesium sulfate therapy; these are used to relieve symptoms and prolong the gestational age. However, these therapies have a minimal effect in preventing the development of PE.[ A subset of patients experiences severe PE symptoms, leading to termination of pregnancy.[

Studies have shown that PE progression can be divided into 2 stages: starting with abnormal placentation and uteroplacental ischemia early in the first trimester and progressing toward excess production of antiangiogenic factors resulting in maternal syndrome.[ In the first stage, inadequate spiral arteriolar remodeling may cause intermittent hypoxia and reoxygenation contributing to oxidative stress.[ Rajakumar et al[ found overexpression of HIF-1α and HIF-2α (markers of cellular oxygen deprivation) in proliferative trophoblasts and placental tissue of patients with PE. AP39, a mitochondrial targeting hydrogen sulfide donor, was found to decrease the expression of sFLT1 (a classic antiangiogenic factor) in placental tissues of PE patients to facilitate instability of HIF-1α.[ In the later stage, there is accumulating evidence to suggest that PE is a proinflammatory disease with insufficient trophoblast invasion, which is characterized by imbalance of interleukin-10 and proinflammatory cytokines, leading to failure of T-cell differentiation.[ The advent of microarray and high-throughput sequencing technologies have provided a newfangled method to further understand the genetic and molecular mechanisms of PE. Gong et al[ conducted high-quality RNA-seq study and found abundant expression of several circ-RNAs, such as hsa-circ-0007444 in PE cases. Moreover, miR-210 was shown to inhibit the proliferation, invasion, and angiogenesis of trophoblasts through targeting mitogen-activated protein kinase (MAPK),[ Notch homolog 1,[ and potassium channel modulatory factor 1.[ However, the precise etiopathogenetic mechanisms of PE are not clear and there is an immense need for identifying high risk groups and discovering novel therapeutic targets.

Here, we acquired data from 4 datasets (GSE10588, GSE25906, GSE60438 and GSE75010) from the Gene Expression Omnibus (GEO). We then divided PE patients into 3 subgroups and constructed a weighted coexpression network analysis of subgroup-specific differentially expressed genes to investigate the genes associated with maternal and fetal complications. Significant gene coexpression modules and genes related to clinical characteristics were identified. Furthermore, the brown, turquoise, and blue module gene network was constructed and hub genes were corroborated using Cytoscape. Finally, miRNAs related to these module hub genes were detected through miRNet, which may serve as biomarkers and therapeutic targets of PE.

2. Materials and Methods

2.1. Data preprocessing

Figure 1 illustrates the flowchart of the data preparation, processing and analysis. As shown in Table 1, 4 datasets (GSE10588, GSE25906, GSE60438, and GSE75010) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Batch effects are subgroups of measurements that exhibit qualitatively different behavior across conditions and are unrelated to the biological or scientific variables in a study.[ The R package “sva” was used to remove the batch effect among these platforms and the principal component analysis was conducted to evaluate the efficiency of batch-effect removal. Given that this study was based on the public database without human or animal experiments, ethical approval and patient consent were not needed.

Figure 1.

The flowchart of the data preparation, processing, and analysis. GO = Gene Ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, WGCNA = weighted gene coexpression network analysis.

Table 1

Summary data of GEO database.

Dataset ID	Platform	No preeclampsia	Preeclampsia
GSE10588	GPL2986	50	21
GSE25906	GPL6102	37	23
GSE60438	GPL6884	65	60
	GPL10558
GSE85307	GPL6244	101	47
	GPL19162
GSE75010	GPL6244	173	157

GEO = Gene Expression Omnibus.

2.2. Verification of subgroup-specific genes and clinical characteristics among 3 subgroups

Samples were classified into different subgroups using consensus clustering, which was performed by K-means algorithm with the spearman distance. The maximum cluster number was 10 and the final cluster number was decided by the consensus matrix and cluster consensus score. Then, Wilcoxon sum-rank text was used to verify the differential expression of subgroup-specific genes among cases of each subgroup with the thresholds of Benjamini–Hochberg adjusted P < .05. The R package “ggpubr” was used to assess differences of maternal age, body mass index (BMI), gestational age, and infant z-score[ between the subgroups.

2.3. Weighted gene coexpression network analysis construction

We used R package termed “weighted gene coexpression network analysis” (“WGCNA”)[ to construct the coexpression network of the filtered genes, which could characterize the biological function of each group. In this study, the soft-thresholding power was set to β = 5 (scale free R2 = 0.9). The topological overlap matrix[ was used to cluster the genes by analyzing the average linkage and Pearson correlation. Next, we constructed a cluster tree with the following major parameters (cutHeight = 0.4; minimum module size = 55; deepSplit value = 2) and identified functional modules using hierarchical clustering analysis. Finally, the correlation between genes in functional modules and clinical characteristics in the WGCNA package were assessed using Spearman correlation coefficient; P values of < .05 were considered indicative of statistical significance.

2.4. GO and KEGG analysis of Module genes

Gene set enrichment analysis (GSEA) was performed using GSEA software (version 4.1.0). The Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed with the module genes. The R package “clusterProfiler” was used to analyze GO and KEGG enrichment.[ Terms with false discovery rate of <0.05 and P value of <.05 were considered for subsequent analysis and the dotplot was processed using the R package “ggplot2.” The R package “pheatmap” was performed to create the heatmap of subgroups with KEGG significant terms.

2.5. Hub gene identification and module miRNA–mRNA network construction

Using Maximal Clique Centrality methods of cytoHubba package, Cytoscape software (version 3.8.0) was utilized to rank the candidate genes and the top 10 nodes were selected as the hub genes. By potentially affecting many biological processes (BPs), miRNAs, with potent posttranscriptional regulatory functions, are closely related to the occurrence and development of diseases. Each miRNA would target hundreds of different genes and multiple miRNAs could regulate their target genes cooperatively. Through miRNet analysis (www.mirnet.ca) in human tissue PLACENTA, based on 11 well-annotated miRNA databases, we could build miRNA-target interaction network and rank the miRNAs by number of interactions, which are most relevant to the hub genes studied.[ Networks were visualized using Cytoscape software.

2.6. Statistical analysis

R software (v4.0.3) was used for statistical analysis. The correlation between mutant genes and TMB was determined by the Mann–Whitney U test. The enriched functions associated with the differentially expressed genes were determined by the R package “clusterProfiler.” Two-tailed P < .05 was considered statistically significant for all comparisons.

3. Results

3.1. Clinic characteristics of PE subgroups

A total of 385 samples (180 PE subjects and 205 no-PE cases) from 4 independent GEO datasets (GSE10588,[ GSE25906,[ GSE60438,[ and GSE75010[) were analyzed to predict the PE subgroups by evaluation of gene signatures. Principal component analysis method was performed to remove the batch effect among different platforms and datasets (Fig. 2A). The scatterplot revealed successful removal of the batch effect (Fig. 2B). Then, 180 PE cases were classified into 3 subgroups by consensus clustering analysis of gene expression profiles. The minimal cluster consensus score was >0.6 and the cluster count ranged from 2 to 10 (Fig. 2C). There were 77, 88, and 15 cases in subgroups I, II, and III, respectively. The consensus matrix revealed a high similarity of gene expression (Fig. 2D).

Figure 2.

PCA and consensus clustering analysis of the gene expression profiles of PE cases. The colors indicate the samples from 4 different datasets. Each dot in the scatter plots represents the samples on the top 2 principal components (PC1 and PC2) before (A) and after (B) batch-effect removal of gene expression profiles. (C) The bar-plots visualize the consensus scores for subgroup with cluster count ranging from 2 to 10. (D) The heatmap shows the consensus matrix with cluster count of 3, which was determined by the minimal consensus scores for subgroup (>0.6). PC = principal component, PCA = principal component analysis, PE = preeclampsia.

Next, we analyzed the differences between the 3 subgroups with respect to maternal age, BMI, gestational age, and infant z-score. We discovered that the patients in subgroup II were older and their gestational age was shorter than those in subgroups I and III (Fig. 3A, B); this indicated that subgroup II may develop more severe PE. However, no significant differences were observed with respect to maternal BMI or infant z-score among the 3 subgroups (Fig. 3C, D). We further clarified the relationship of subgroup with maternal age and gestational age using analysis of variance; the results indicated that the PE subgroups were a maternal age-independent indicator of gestational age (Table 2; P = .021).

Figure 3.

Box-plots representing the pairwise comparison of clinical features, including maternal age (A), gestational age (B), maternal body mass index (C), and infant z-score (D) between the subgroups. *P < .05, **P < .01, ns: no significance. BMI = body mass index.

Table 2

Analysis of variance of gestational age, fetal z-score, and their interaction in the subgroups.

	Df	Sum square	Mean square	F value	Pr. (>F)
Subgroup	2	109.7	54.83	4.001	0.0205*
Maternal age	1	8.1	8.14	0.594	0.4421
Subgroup: maternal age	2	39.8	19.90	1.453	0.2376
Residuals	134	1836.2	13.70

Df = degree of freedom.

Significance codes: *P < 0.05.

0.05.

3.2. Construction of weighted gene coexpression network and modules detection for each subgroup

A total of 2845 subgroup-specific genes (991, 376, and 1478 genes in subgroup I, II, and III, respectively) were obtained by pairwise differential expression analysis across subgroups (Benjamini–Hochberg adjusted P < .05; mean difference threshold = 0.2; Table 3; Table S1, Supplemental Digital Content, http://links.lww.com/MD/G938). We further performed differential expression analysis to identify the upregulated genes in each subgroup. GSEA indicated that specific subgroup genes were significantly expressed when compared with normal samples (false discovery rate <0.05; Fig. 4A–C). In addition, subgroup III contained the largest number of subgroup-specific genes and differentially expressed genes compared with normal samples; considering that these patients had moderate PE symptoms, these findings indicated that specific categories of functional genes in subgroup III may have a protective effect against PE.

Table 3

Data of differentially expressed genes in case–control and case–case comparative analysis and the subgroup-specific WGCNA modules.

Subgroup	Genes in case–control comparison	Gene in case–case comparison	WGCNA modules
I	318	991	Green, turquoise
II	1124	376	Brown
III	3476	1478	Blue, gray

WGCNA = weighted gene coexpression network analysis.

Figure 4.

Expression profiles of subgroup-specific genes. Enrichment plots of (A), (B), and (C) indicate that the subgroup-specific genes have different expression levels than the normal controls.

The R package “WGCNA,” a well-known data mining tool, was utilized to construct a weighted coexpression network of these subgroup-specific genes. No outliers or strong clusters were observed; therefore, the selected samples were well clustered in hierarchical cluster analysis. The soft-thresholding power and soft threshold 5 were considered as the optimal power with a scale-free topology fix index 0.90 (Fig. 5A, B). Under specific conditions, the “blue” and “yellow” modules and the “brown” and “red” modules were merged (Fig. 5C, D). Thus, 5 gene coexpression modules were identified in the WGCNA analysis. The correlation of these 5 modules with PE clinical characteristics was assessed to identify the module that showed the most significant level of correlation (Fig. 5E). The module “turquoise” showed the strongest positive correlation with infant z-score and gestational age, while module “brown” showed a significant negative correlation with these characteristics. Module “blue” showed a positive correlation with gestational age and infant z-score, but showed no relationship with maternal age. Other modules, such as green and gray, showed no significant association with maternal age, gestational age, or infant z-score.

Figure 5.

Module detection, gene significance, and module membership of PE samples. (A) Network topology analysis for different soft-thresholding powers. (B) The scale-free fit index and the mean connectivity as a function of the soft-thresholding power. (C) The clusters of the modules, the eigengenes, and the threshold (red line) of modules were merged. (D) Clustering dendrogram of genes with dissimilarity based on the topological overlap, together with assigned module colors. The unmerged colored modules in first row below the dendrogram and the merged colored modules are in the second. (E) The association of the modules and the traits was constructed. Each row refers to a module eigengene and each column refers to a trait. Numbers in each cell represent the correlation and P value. The right side of the heatmap reflects the intensity and direction of the correlation (the color depth indicates the strength of the correlation; red, positively correlated; green, positively correlated). BMI = body mass index, PE = preeclampsia.

3.3. Module function enrichment analysis

We used GO analysis (including cell component, BP, and molecular function) to perform enrichment analysis of the biological significance on the 5 WGCNA modules (Fig. 6A–D). In BP group (Fig. 6A, B; Table S2, Supplemental Digital Content, http://links.lww.com/MD/G939), module blue was mainly enriched in neutrophil activation involved in immune response, including neutrophil activation involved in immune response, neutrophil degranulation, T-cell activation, lymphocyte proliferation, and mononuclear cell proliferation. Module brown showed a significant association with multi–multicellular organism process, covering cell–substrate adhesion, response to metal ion, tissue migration, multi–multicellular organism process, female pregnancy, reproductive structure development, and reproductive system development. Module turquoise was mainly associated with extracellular structure organization, extracellular matrix organization and cell-substrate adhesion. In cell component group (Fig. 6A, C; Table S3, Supplemental Digital Content, http://links.lww.com/MD/G940), module blue was enriched in specific granule, secretory granule membrane, tertiary granule, specific granule membrane, and external side of plasma membrane. Module brown was mainly associated with membrane microdomain, membrane region, and membrane raft. Module turquoise showed a significant association with focal adhesion, cell–substrate junction, collagen-containing extracellular matrix, and endoplasmic reticulum lumen. In molecular function group (Fig 6A, D; Table S4, Supplemental Digital Content, http://links.lww.com/MD/G941), genes of module blue were mainly associated with immune receptor activity, carbohydrate binding, protein serine/threonine kinase activity, and signaling receptor activator activity. Module brown was enriched in protein serine/threonine kinase activity and signaling receptor activator activity. Module turquoise was mainly related to extracellular matrix structural constituent.

Figure 6.

The expression patterns of subgroup-specific genes and functional characterization of WGCNA modules. (A) The scaled expression values of genes that comprise of the 5 gene coexpression network analysis modules are displayed in the heatmap. The bubble plots show the top significantly enriched genes in GO analysis including biological process (B), cellular component (C), and molecular function (D), as well as the KEGG pathways for each WGCNA module (E). (F) The heatmap demonstrates the most significant KEGG pathway in each WGCNA module associated with each subgroup. GO = Gene Ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, WGCNA = weighted gene coexpression network analysis.

The results of KEGG analysis of the 5 WGCNA modules are shown in Figure 6E and F and Table S5 (Supplemental Digital Content, http://links.lww.com/MD/G942). Module blue was enriched in chemokine signaling pathway, osteoclast differentiation, hematopoietic cell lineage, natural killer cell-mediated cytotoxicity, viral protein interaction with cytokine and cytokine receptor, MAPK signaling pathway, and Rap1 signaling pathway. Module brown was enriched in MAPK signaling pathway, Rap1 signaling pathway, proteoglycans in cancer, focal adhesion, cAMP signaling pathway, and breast cancer. Module turquoise was enriched in focal adhesion, MAPK signaling pathway, Rap1 signaling pathway, proteoglycans in cancer, vascular smooth muscle contraction, and cGMP-PKG signaling pathway. Notably, MAPK signaling pathway and Rap1 signaling pathway were both enriched in these 3 key modules.

3.4. Hub gene analysis and prediction of potential miRNA–mRNA networks

According to the relatedness, as assessed by Cytohubb package, Figure 7A to C shows the specific top 10 hub genes in modules blue, brown, and turquoise from strong to weak. The 10 hub genes in module blue were C3AR1, CCR2, CCL5, CCL4, CXCL9, CXCR4, CCR7, CX3CR1, CXCL16, and CXCL11. In module brown, the top 10 hub genes were EGFR, MAPK3, FGF7, FGF10, FGF9, FLT1, FLT4, FN1, PIK3CB, and PPARG. The 10 hub genes in module turquoise were COL6A1, COL4A1, COL4A2, COL6A2, COL1A2, COL3A1, COL5A2, COL6A3, COL15A1, and COL1A1.

Figure 7.

Interaction network of hub genes and potential miRNA–mRNA regulatory networks. (A–C) Interaction network of the top 10 respective hub genes in the module turquoise, brown, and blue. Using Cytoscape package to obtain a network from DEGs based on WGCNA method and cytoHubba package to search the hub genes in R Studio. The color depth indicates the strength of hub genes. Low to high significance ranked by colors from yellow to red. (D–F) Interaction networks of predicted miRNAs and their target genes; blue squares represent miRNAs, yellow squares represent hub gene mRNAs. Edges between squares denote interactions between miRNAs and hub genes mRNAs. DEG = differentially expressed genes, WGCNA = weighted gene coexpression network analysis.

Moreover, miRNet database was utilized to predict the networks of hub genes and potential miRNA in the 3 key modules. To the hub genes in module brown, the top 5 miRNAs (ranked by degrees) of 36 miRNAs predicted were hsa-mir-34a-5p (degree 7), hsa-mir-19b-3p (degree 4), hsa-mir-522-5p (degree 3), hsa-mir-210-3p (degree 3), and hsa-mir-145-5p (degree 3; Fig. 7D). In module turquoise, the top 6 miRNAs of the 52 miRNAs predicted by the hub genes were hsa-mir-145-5p (degree 6), hsa-let-7a-5p (degree 6), hsa-mir-34a-5p (degree 4), hsa-mir-106a-5p (degree 4), hsa-mir-19b-3p (degree 3), and hsa-mir-498 (degree 3; Fig. 7E). In module blue, the top 3 miRNAs of the 20 miRNAs predicted by the hub genes were hsa-mir 34a-5p (degree 3), hsa-mir-210-3p (degree 2), and hsa-mir-106a-5p (degree 2; Fig. 7F).

4. Discussion

In the present study, based on the analysis of the gene expression profiles of 205 normal controls and 180 PE cases from 4 independent GEO datasets, we classified the PE cases into 3 subgroups after eliminating the batch effect and revealed the subgroup-specific gene modules associated with clinical features using WGCNA. By analyzing the correlation between module genes and clinical characteristics, we identified 3 key gene modules and their hub genes. Moreover, we successfully constructed the hub gene-miRNA ceRNA networks.

As a placental disease, the pathogenesis of PE is not well characterized. The suggested pathogenetic mechanisms include reduced invasive capacity of placental trophoblasts, recast dysfunction of the spiral uterine artery caused by an excess of antiangiogenic factors, and inappropriate immune responses to the allogenic fetus.[ Previous microarray studies that have sought to identify early detection markers and candidates for therapeutic intervention have indicated the involvement of oxidative stress and cytokine-associated genes.[ However, there is considerable heterogeneity in the clinical manifestations and therapeutic effect, which cannot be explained by classical separation of early-onset and late-onset (diagnosis after 34 weeks) PE patients. Leavey et al[ divided PE patients into 5 clusters. They found high levels of PE markers (such as FLT1) in cluster 2, while cluster 3 showed a maternal-fetal incompatibility presentation; this may represent the considerable interindividual variability observed in clinical settings.

Likewise, in our study, we divided PE patients into 3 subgroups with consensus clustering and the more severe subgroup II had greater maternal age and shorter gestational age. Genes in module brown were significantly upregulated in subgroup II and functional enrichment analysis showed that pathogenesis of this subgroup was mainly related to MAPK, Cyclic adenosine monophosphate (cAMP), and Ras-associated protein 1 (rap1) signaling pathways. Studies have demonstrated the participation of these 3 pathways in the regulation of trophoblast invasion and migration in normal pregnancy.[ Padmini et al[ revealed that aberrant expression of extracellular signal-regulated kinase 1/2 (ERK1/2, which is also known as MAPK3 and MAPK1), a member of MAPK family, plays an important role in PE progression. Studies have also shown that the downregulation of ERK signaling can aggravate oxidative injury in trophoblasts.[ Moreover, PE cases showed oxidative stress-induced downregulation of p38 MAPK, which may be caused by elevated phosphorylation of p38 MAPK, leading to defects in growth of placenta.[ In pregnancy, expression of epidermal growth factor receptor activated MAPK/ERK1/2 signaling pathway in placental tissues to regulate blastocyst implantation may result in occurrence of PE.[ Consistent with these results, our study demonstrated that these MAPK signaling pathway members including epidermal growth factor receptor and ERK1 (MAPK3) may serve as hub genes in module brown. In addition, Chen et al[ found that hypoxia, accompanied by dysfunction of placental multipotent mesenchymal stromal cells, induced cAMP decrease and Rap1 inactivation in trophoblasts, which resulted in decreased trophoblast adhesion and migration during the PE process. Members of the fibroblast growth factor (FGF) family, which were detected in our hub gene analysis as FGF7, FGF9, and FGF10, were shown to participate in the MAPK/ERK/cAMP signaling pathway to modulate uterine–conceptus interactions.[ In a previous study, FGF10 expression was shown to be localized to first trimester extravillous trophoblasts, indicating its potential role in PE pathophysiology.[ In our study, older maternal age and shorter gestational age were simultaneously observed in subgroup II, suggesting that PE women with older age may manifest more serious disease leading to premature birth. Moreover, we also identified classical markers such as FLT1, FLT4, and FN1 as hub genes, and this suggests that the pathogenesis of PE in subgroup II patients may have a typical placental origin.

In subgroup III, genes in module blue and gray were upregulated and module blue showed a significant positive association with gestational age and infant z-score. Our gene enrichment analysis showed that these genes were enriched in chemokine signaling pathway and immune-related BPs such as nature killer (NK) cell-mediated cytotoxicity, neutrophil, T-cell activation, and lymphocyte proliferation. The hypothesis that the immune system is involved in the development of PE has already been confirmed by many reports.[ Autoimmune response may alter the maternal immune pathology, interfering with adaptions for placentation, and increase the incidence of PE as well as intrauterine growth restriction.[ A typical characteristic of the maternal immunological preparatory response to pregnancy is a T-helper type-2 lymphocyte (Th2) proportion increase simultaneously with Th1 upregulation.[ A suboptimal interaction between uterine NK cells and the trophoblast is believed to lead to maternal spiral artery modification alteration in patients with PE.[ Extensive infiltration of activated neutrophils in blood vessels of PE women may lead to systemic vascular inflammation and multiorgan involvement. Walsh et al[ discovered that cytokine interleukin-17A may increase neutrophil chemokines like CXCL5 and CXCL6 gene expression in vascular smooth muscles and result in neutrophil boost in PE pregnancy. The balance expression of chemokines at the maternal-fetal interface affects the immune cell profile and function in the decidua.[ In our analysis, most of detected hub genes are chemokines associated with multiple immunocytes, such as monocytes (CCR2), T cells (CCL5, CXCL11, and CXCL16), Th1 cells (CXCL9, CXCR4), and NK cells (CCL4). In addition, decreased expression of another hub gene, C3AR1, was also discovered in preeclamptic pregnancies with fetal growth restriction, indicating an inability of preeclamptic placentas to clear the excess of activated complement, promoting accumulation of immune complexes and inflammation.[ Interestingly, subgroup III showed mild PE clinical features while containing most subgroup-specific genes and differentially expressed genes compared with normal cases; this suggests the existence of a protective mechanism in these patients, which calls for further research. Accordingly, our results indicate that subgroup III may be interpreted as a maternal-fetal incompatibility coupled with a maternal immune response hyperfunction.

Subgroup I had the least number of differentially expressed genes with the normal controls, suggesting that subgroup I may exhibit normal-like expression patterns. Compared to subgroup III, no significant differences were observed with respect to gestational age or maternal age. However, the intrinsic biological characteristics in subgroups I and III showed great differences. Subgroup I showed a significant correlation with module turquoise and genes enrichment analysis showed that most genes were enriched at focal adhesion. Our hub gene analysis showed that all hub genes in module turquoise are collagen family member genes. Collagen is one of the most abundant components of extracellular matrix which could influence and regulate trophoblast invasion and participate in the remodeling of the decidua at the maternal–fetal interface.[ Studies have shown that excessive collagen deposition may affect smooth muscle growth and decrease the remodeling of spiral arteries, which may be an important factor in the pathogenesis of early-onset PE.[ A previous study identified aberrant expression of collagen fragments in PE; COL4A1 and COL4A2, which code for collagen type IV α chain 1 and chain 2, respectively, were identified as maternal PE susceptibility genes.[ Ohmaru-Nakanishi et al found significant upregulation of COL1A1 expression in fibroblasts from PE placentas compared with normal placentas; in addition, placental fibroblasts were also shown to increase the production of collagen type I and IV in hypoxic conditions in vitro.[ Thus, subgroup I may represent an earlier disease stage based on its normal-like expression patterns and early-onset collagen abnormal expression.

Finally, miRNA–mRNA network was constructed to detect possible miRNA involvement. miR-34a-5p and miR-106a-5p were found to be correlated with all 3 significant gene modules, indicating their potential involvement in the pathogenesis of PE. Xue et al[ reported upregulation of miR-34a-5p in placental tissues and serum of patients with severe PE; they found that miR-34a-5p inhibited the invasion and migration of trophoblast cells targeting Smad4. In previous studies, miR-34a elevation was shown to inhibit trophoblast proliferation, migration, and invasion and to simultaneously promote apoptosis.[ However, miR-106 was reported to promote proliferation, migration, and invasion of cytotrophoblast via blocking syncytiotrophoblast differentiation.[ Our results indicated that miR-19b-3b and miR-210-3p participated in the network of module blue and turquoise. Studies have shown downregulation of miR-19b-3p and upregulation of miR-210-3p in the peripheral circulation of women with PE compared with healthy pregnant women.[

Our findings suggest that PE patients belonging to different subgroups may have different genomic patterns. Further studies are required to validate these findings, especially for the miRNA–mRNA predictions. In conclusion, our study observed heterogeneity between PE subgroups and attempted to explain the heterogeneity according to different gene expression patterns. Our findings suggest that patients in unique subgroups may benefit from more personalized treatment. The predicted hub genes and miRNAs identified in this study may help develop new therapeutic targets in the future.

Author contributions

Conceptualization: Min Zhang, Xiaheng Deng, Ziyan Jiang.

Methodology and Administrative support: Zhiping Ge.

Data curation: Min Zhang, Xiaheng Deng.

Formal analysis and investigation: Min Zhang, Xiaheng Deng.

Data analysis and interpretation: Min Zhang, Xiaheng Deng, Ziyan Jiang.

Writing – review & editing: Min Zhang, Xiaheng Deng, Ziyan Jiang, Zhiping Ge.