Gene Expression Network Analysis Identifies Potential Targets for Prevention of Preeclampsia

Introduction

Preeclampsia (PE) is a common hypertensive disorder complicating pregnancy, with an incidence rate of approximately 2–5% worldwide.^1,2 While 14% of all maternal deaths are caused by PE, it becomes the main reason for the increasing maternal and perinatal mortality.³ PE leads to maternal heart impairment. Compared with pregnant women with normal blood pressure, those with a history of gestational hypertension have a significantly increased risk of kidney and cardiovascular diseases and even impaired cognitive function in the later stage.^4,5 In addition, PE may aggravate fetal intrauterine hypoxia, resulting in premature delivery,⁶ fetal growth retardation, or defective fetal heart development leading to fetal cardiac dysfunction.⁷ The etiology of PE is multifactorial and polygenic, which remains elusive despite intensive research. According to the clinical features, several subtypes of PE can be sorted. In general, the development of PE can be divided into two stages: asymptomatic and symptomatic stages.⁸ Abnormal placental development in early pregnancy causes placental dysfunction and later leads to the release of factors into maternal blood, resulting in hypertension and organ impairment. The occurrence of PE is closely related to placental dysfunction. The researchers conducted in-depth research into the placenta as the root cause of PE, and as the research progressed, a more refined classification of PE stages was made.⁹

However, the molecular mechanisms underlying poor placentation remain largely unknown, although increasing research studies have applied RNA sequencing and other analytical approaches to explore placenta development. Some studies are sorting out the relationship between PE subtypes and specific genes.^10,11 By analyzing the RNA profiles of 302 human placenta samples, Gong et al¹² demonstrated that elevated serum levels of follistatin-like 3 (FSTL3) in pregnant women were predictive of subsequent PE and fetal growth restriction (FGR). Moreover, reported by recent studies, plasma cell-free RNA (cfRNA) could exhibit specific patterns to indicate normal pregnancy progression and serve as the biomarker to detect the risk of PE months before clinical manifestations. While the blood sample collection only requires venipuncture for once, cfRNA signatures can track pregnancy progression at the placental level, on the maternal or fetal sides, and can effectively predict the occurrence of PE, with a sensitivity of 75% and a positive predictive value of 32.3%.¹³ However, a comprehensive analysis on the clinical data and genomic profiling in PE is still lacking. To protect the maternal and fetal health against PE, early detection of PE and timely treatment for pregnant women with PE are crucial, thereby avoiding adverse pregnancy outcomes. Early detection and timely treatment of pregnant women with PE are crucial for avoiding adverse pregnancy outcomes. Accordingly, it is necessary to identify early diagnostic markers and develop effective treatment methods. In this study, based on comparing the transcriptional profile of placental samples and the clinical information, the purpose was to evaluate the relationship between specific genotypes and clinical features of PE, which might contribute to the identification of candidate biomarkers for the prediction and prevention of PE.

Data and Methods

Data Download and Processing

Clinical information and transcriptional data were downloaded from GEO database (https://www.ncbi.nlm.nih.gov/gds). Inclusion criteria were as follows: 1) Pregnancy-induced hypertension, including gestational hypertension, preeclampsia, eclampsia, chronic hypertension with preeclampsia, and pregnancy complicated with chronic hypertension; 2) Samples with complete clinical information and human transcriptional data. Perl (V 5.18.4) was used to convert the probe name into gene name in downloaded data. R (R 4.1.1) language in combination with “limma” and “SVA” packages were employed to merge the downloaded data, take the intersection of the sorted genes, compute the mean value of the same gene, take log2 for the data with large gene expression, and carry out batch correction to eliminate the impact of batch effect. R language “ggplot2” package for PCA analysis was run to verify the effect of batch correction.

Classification, Correlation Analysis of Continuous Variables, and Classification Difference Analysis

The integrated data were classified by using “consumusclusterplus“ package. The clinical data were sorted, and the correlation between age and gestational week of delivery was analyzed. According to the typing results, taking the mean difference >0.2 and P < 0.05 as the screening criteria, R language ”limma” package was used to take the genes of different typing samples and were compared with those of normal samples, and potential biomarkers in different groups were screened.

GSEA and WGCNA

The differential genes were obtained, enriched and analyzed by gene set enrichment analysis (GSEA) (V4.1.0).¹⁴ Thereafter, the differences in specific differential genes between each typing and normal samples were compared. In the meantime, R “WGCNA” package was used to calculate the threshold (optimal weight) of the co-expressive gene correlation network. Then, gene function modules were identified by dynamic segmentation algorithm (cutree dynamic), and genes with a similar expression pattern were assigned into one module. Finally, the dissimilarity between modules was calculated based on the module feature vector. The similar modules were combined, and the correlation between WGCNA gene function module and clinical features of samples was analyzed.

GO and KEGG Pathway Enrichment Analyses

GO and KEGG pathway enrichment analyses were performed by using R language software, Bioconductor plug-in, and clusterprofiler. Data were processed, and genes without ID were deleted. P < 0.05 was set as the standard for screening. The enrichment of different modules and genes was visually displayed, and R was run to obtain the bubble diagram of GO and KEGG enrichment analysis.

Results

Data Inclusion and Processing

According to the inclusion criteria, two PE chip datasets GSE60438 and GSE75010 were obtained from the GEO database and included in this study. The basic information of the two gene matrices were listed in Table 1, and the clinical data were listed in Tables 2 and Tables 3.

Table 1 Gene Expression Dataset of PE in GEO Database

Table 2 Clinical Data of the GSE60438 Dataset

Table 3 Clinical Data of the GSE75010 Dataset

PCA Verification and Sample Typing

We conducted PCA analysis on the corrected data. As depicted in Figure 1, samples from GSE60438 and GSE75010 prior to correction were distributed on both sides of the figure (A), while those after correction were randomly distributed (B). This observation indicated that the influence of batch effect has been eliminated. We, therefore, typed the data by using R language. According to the consistent scoring results, the genes were divided into two groups. The analysis identified a weak correlation between different genotypes (C).

Figure 1 PCA validation and sample typing. (A) PCA diagram before correction; (B) PCA diagram after correction; (C) The cluster of PE subtypes.

Correlation Analysis of Continuous Variables

We next employed R language to analyze the correlation between the age of pregnant women, gestational week of delivery and the classification. As shown in Figure 2, there were differences in the age of pregnant women and gestational week of delivery among different groups. According to the grouping results, we organize and analyze the genes of different groups. In group 1, the following genes had smaller p-values: SASH1, PIK3CB, FLT-1, FLNB, FAM120A, ENG, MICAL3, SH3PXD2A, MYO7B. While in group 2 is SMARCA1, TSC22D1, NTRK2, HSD17B1, GSN, EPHB6, APLN, SULF2, CETN2.

Figure 2 Correlation Analysis between maternal age, delivery gestational week, and classification. (A) Correlation analysis between maternal age and classification; (B) Correlation between gestational age and classification. Note: *P < 0.05.

GSEA and WGCNA

We performed GSEA to examine specific differential genes in each group. As illustrated in Figure 3, the differential genes were different between genotyping samples and normal samples (A, B). We then constructed a gene co-expression module network by WGCNA, based on the differential genes. A total of 1356 differential genes in 205 samples were included in this study. TOM cluster analysis revealed an association of the gene module with the known phenotype, and the correlation coefficient between the feature vector of the module and the phenotype was calculated. Notably, we found that ME turquoise module was positively correlated with maternal age, maximum systolic blood pressure, maximum diastolic blood pressure, and the number of urinary proteins (C). Moreover, a heatmap of differential genes in the module showed that while the group 1 genes were up-regulated in turquoise module, the group 2 genes were up-regulated in blue and brown modules (D). These data indicated that genes in turquoise module and group 1 were closely related to PE. Among those genes, according to the low p-values, three candidate genes were selected for subsequent analysis: SASH1, PIK3CB and FLT-1.

Figure 3 GSEA and WGCNA. (A) GSEA analysis on specific differential genes between group 1 and normal samples; (B) GSEA analysis on specific differential genes between group 2 and normal samples; (C) Module–trait relationships. (D) A heatmap of module genes.

GO and KEGG Pathway Enrichment Analyses

We further performed GO and KEGG pathway enrichment analyses by using R language (Figure 4). GO enrichment analysis revealed that differential genes in the blue module and brown module were mainly involved in biological processes such as T cell activation and those including extracellular matrix organization, respectively, while the genes in the turquoise module are mainly associated with biological processes such as small molecular catalytic process (A). Meanwhile, KEGG pathway enrichment analysis showed that while differential genes in the blue module and brown module were highly enriched in pathways such as those involving cell adhesion molecules and pathways including those involving apoptosis, respectively, pathways of differential genes in turquoise module mainly included MAPK signaling pathway, Rap1 signaling pathway, and pathways involving valine, leucine and isoleucine degradation (B).

Figure 4 Bubble Diagram of GO and KEGG enrichment analysis. (A) GO enrichment analysis bubble chart; (B) KEGG pathway enrichment analysis bubble chart.

Discussion

It has been shown that under the condition of PE, great changes such as oxidative stress and inflammatory response occur in the intrauterine environment of pregnant women, leading to alterations in various signaling pathways. PE has a great impact on growth and development of the fetus, impairs the structure and function of fetal cardiovascular system or kidney, and causes hypertension.¹⁵ Current studies suggest that the pathogenesis of PE may be related to decidual macrophage dysfunction,¹⁶ placental angiogenesis disorder,¹⁷ placental aging,^18,19 protein aggregation, and defective placental aggregation phagocytosis.²⁰ At the molecular level, the microRNA spectrum of blood exosomes has been used to evaluate the pathophysiology of preeclampsia and predict the disease,²¹ and there are more and more studies at gene level.^22–25

Zhang and Horvath²⁶ conducted the first WGCNA analysis in 2005 and then published WGCNA R software package in R based on theory of the research method in 2008. As an efficient and accurate bioinformatics approach, WGCNA has been widely used for the gene network analysis.²⁷

Herein, we intensively analyzed the published sequencing data in selected two datasets, including 205 samples (100 non-PE and 105 PE cases) with complete clinical phenotype information and transcriptome sequencing data. In the study, we constructed the co-expression network modules and performed clinical phenotype correlation analysis on the mRNAs obtained by placental tissue sequencing. The analysis revealed that MEturquoise module was positively correlated with maternal age, maximum systolic blood pressure, maximum diastolic blood pressure, and the number of urinary proteins. Further studies showed that the patients with group 1 are consistent with this law. The present study provided evidence that the up-regulated differential genes SASH1, PIK3CB and FLT-1 may affect the progression of PE through MAPK signaling pathway and Rap1 signaling pathway.

Serebrov et al²⁸ reported that rs34845949 of the SASH1 gene was significantly correlated with PE. The expression levels of SASH1 in placental tissues were increased in patients with PE compared with healthy pregnancy,^10,28,29 suggesting that SASH1 gene might play a role in placental dysplasia in patients with preeclampsia. Furthermore, reported by published studies, the overexpression of SASH1 could inhibit the proliferation of HTR-8/SVneo cells, possibly due to the inhibition on CCND1 and CCND3, which arrested the cells in G1 phase. And the expression levels of CCND1 and CCND3 were associated with MAPK signaling.^29,30 Ras/Raf/MEK/ERK pathway in MAPK can be activated by a variety of stimuli, and the activated Ras/Raf/MEK/ERK pathway regulate numerous physiological processes such as cell growth, development, division, and apoptosis. It has been demonstrated that SASH1 interacts with MAP2K2 and MAP4K4. While MAP4K4 indirectly activates MAP2K2 via MEKK1, it may affect receptor signaling through Rac1-MEKK1-MAP4K4-JNK pathway and play an important role in cell migration and invasion.³¹ Human placental pluripotent mesenchymal stromal cells regulate trophoblast migration via Rap1 activation.³² Previous studies have shown that inhibition of MAPK signaling pathway suppresses trophoblast proliferation and migration, probably promoting the progression of PE.^33,34 Based on the above observations, we speculate that SASH1 may be involved in the progression of PE by affecting the proliferation and migration of trophoblasts via MAPK signaling pathway and Rap1 signaling pathway. At present, the molecular mechanism enabling SASH1 to participated in the development of PE is still unclear. It is necessary to further study the specific role of SASH1 in PE, especially by analyzing the differential protein expression levels. Based on the relationship between the differential expression patterns of SASH1 and the specific clinical features in PE, the clinical significances of SASH1 can be evaluated, including the potential as a biomarker for early diagnosis and as the target for novel PE therapy.

While PIK3 protein family plays a role in the regulation of various cell functions including cell growth, proliferation and survival, PIK3 is involved in controlling the migration and penetration of human cancer cells and trophoblasts.³⁵ Among the protein family, PIK3CBBCL-2 is considered the core regulatory gene of decidua in eclampsia. Aberrant expression of PIK3CBBCL-2 causes abnormal decidual apoptosis related to the pathogenesis of eclampsia.³⁶ Studies have confirmed that the expression of PIK3CA gene is up-regulated in PE patients as compared to the normal group (P < 0.05).³⁷

Among all potential target genes identified in the present study, FLT-1 gene has been widely used in the animal modeling of PE and clinical trials. In this case, sFLT-1 overexpression is commonly used for generating the animal model of PE, whereas inhibiting FLT-1 can prevent PE.^38,39 It has been shown that sFLT-1/PlGF ratio is valuable in predicting the occurrence of PE.^40,41 Particularly, determination of sFLT-1/PlGF ratio becomes a valuable auxiliary means for PE diagnosis in the medical insurance system of France; a low ratio is regarded as a highly negative predictive value of PE, and unnecessary hospitalization and premature delivery can be avoided.⁴² Numerous studies analyzed the different critical values of the sFLT-1/PlGF ratio, thereby suggesting the ratio as an effective predictor for PE and PE-related complications.^43–47 Moreover, some preceding studies suggested that the pathogenesis of PE might be related to inhibition of MAPK pathway, thereby affecting the invasion of trophoblast cells and promoting the secretion of sFLT-1/sEng.⁴⁸ FLT-1 was found to be closely correlated with the occurrence and development of PE, and FLT-1 was significantly elevated in PE patients.⁴⁹ At present, the effective treatment of preeclampsia mainly relied on the termination of pregnancy. In recent years, some researchers explored the mechanism to reduce the level of sFLT-1 in pregnant women with preeclampsia, so as to restore the balance of angiogenesis. According to previous studies, blood sFLT-1 apheresis, a proton pump inhibitor, might interfere with placental renin-AGT uptake/transport, thereby reducing angiotensin-induced sFLT-1 synthesis;⁵⁰ aspirin nourishes HTR-8/SVneo, by inactivating NF-κB signaling in lamellar cells, attenuated the role of sFLT-1 in oxidative stress and endothelial dysfunction, and reduced trophoblast apoptosis,⁵¹ suggesting that reducing sFLT-1 in PE patients could be applicable and effective. Concentration-related therapeutic approaches is expected to be important applications in the future clinical treatment for the prevention of preeclampsia.

In conclusion, WGCNA and the correlation analysis on published placental mRNA data and clinical information of PE patients revealed that genes in MEturquoise module of group 1 patients were positively correlated with maternal age, maximum systolic blood pressure, maximum diastolic blood pressure, and the number of urinary proteins. Further studies showed that genes SASH1, PIK3CB and FLT-1 were closely related to PE. Moreover, we found that MAPK signaling pathway and Rap1 signaling pathway were highly enriched in KEGG. FLT-1 was conventionally used to predict PE risk, and sFLT-1 could also be used as indicator to evaluate PE treatment effect. As a candidate biomarker for predicting PE, SASH1 could participate in the proliferation, migration, invasion and epithelial mesenchymal transformation of human trophoblast cells, by regulating MAPK pathway and Rap1 signaling pathway, thus affecting the progression of PE. The mechanism allowing PIK3CB to regulate PE development was not clear, although the gene could be another candidate biomarker for PE prediction. This is an exploratory study, and our findings still required verification by further studies.

Furthermore, this study was subject to limitations in the design. Firstly, due to the limitations of the GEO database, the current WGCNA analysis showed that group 1 had a higher correlation with clinical data in genetic data. Therefore, this study only discussed the objects related to group 1, while the genes of group 2 could be further studied in the future; secondly, this study was based on the previous research, which also required further investigation at the molecular level, so the mechanism underlying the regulatory roles of SASH1 and PIK3CB in PE could be further proved in the future.