A Potential Three-Gene-Based Diagnostic Signature for Hypertension in Pregnancy

Introduction

Hypertensive disorders of pregnancy have become the second most common cause of direct maternal death, and affect approximately 5–10% of all pregnancies.^1,2 Chronic hypertension is defined by SBP ≥ 140mm Hg or DBP ≥ 90 mm Hg before pregnancy or before the 20th week of amenorrhea, and gestational hypertension is defined by SBP ≥ 140mm Hg or DBP ≥ 90mm Hg during or after the 20th week of amenorrhea.³ According to the National High Blood Pressure Education Program (NHBPEP) Working Group on High Blood Pressure in Pregnancy, methyldopa, slow release nifedipine, beta blockers (other than atenolol), labetalol, and a diuretic in pre-existing hypertension have been identified to be appropriate therapeutic drugs for hypertension in pregnancy.⁴ However, since the blood pressures of many patients are not controlled by these agents, or other intolerable side effects exist, hypertension in pregnancy is still a challenge for the health of pregnant women.⁵ Therefore, the identification of specific and sensitive diagnostic signatures is urgent, which can contribute to the early diagnosis of hypertension in pregnancy and even expand the window of treatment.

Previous studies have demonstrated that the blood pressure can act as a useful predictive marker for hypertension in pregnancy.⁶ Recently, a series of biological markers involved in hypertension in pregnancy are identified and may also be considered as potential diagnostic signatures. Sun et al identified four differentially expressed genes including KIR3DL2, AKR1C3, CHURC1, and SLC25A13 between the early onset preeclampsia patients and uncomplicated pregnancies controls, which might be beneficial for the diagnosis of early onset preeclampsia.⁷ Yuen et al found that TIMP3 was hypomethylated in early onset preeclampsia samples compared with control samples through DNA methylation profiling analysis and could be used as a promising signature for the diagnosis of early onset preeclampsia.⁸ Compared with normal controls, lncRNA MIR193BHG and lncRNA GATA3-AS1 showed obviously elevated expression levels in early onset preeclampsia samples and exhibited significant value in the diagnosis of this disease.⁹ Velicky et al demonstrated that pregnancy-associated diamine oxidase was significantly downregulated in early-onset preeclampsia and might be used for the early diagnosis.¹⁰ Chaiworapongsa et al analyzed the transcriptomic data of whole blood from patients with late-onset preeclampsia and normal controls, and identified significantly higher expression of CNTNAP3 and lower expression of CCR3 in late-onset preeclampsia samples than those in normal controls.¹¹ Pan et al investigated the relationship between polymorphisms of gene FOXO3 and late-onset preeclampsia, demonstrating that FOXP3 rs2232365 was related to the risk of late-onset preeclampsia, with a possible regulatory role for TG/HDL levels.¹² It was found that PTPRN2 had changes in 5-methylcytosine and 5-hydroxymethylcytosine and was associated with the pathogenesis of late-onset preeclampsia,¹³ providing some references for the diagnostics of this disease. Logistic regression model, often in the presence of other covariates, has been widely used to improve the probability of prediction compared with the predictive markers of single gene.^14,15 Although many potential diagnostic signatures have been identified, the research of classification model to diagnose hypertension in pregnancy is still lacking.

In this study, we identified 20 potential genes involved in hypertension in pregnancy though bioinformatics analysis. Then, a predictive model was constructed based on three optimal genes (LEP, PRL and IGFBP1), which could efficiently separate patients with hypertension in pregnancy from healthy subjects, and might be potentially applied for the early diagnosis of hypertension in pregnancy.

Materials and Methods

Data Collection

We downloaded two mRNA profiles including GSE75010 and GSE48424 from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/). GSE75010 included 80 placenta samples from patients with hypertension in pregnancy and 77 placenta samples from healthy subjects; GSE48424 included 18 blood samples from patients with hypertension in pregnancy and 18 blood samples from healthy subjects. The quantification of mRNA expression in GSE75010 was performed by using Affymetrix Human Gene 1.0 ST Array, and that in GSE48424 was performed by using Agilent-014850 Whole Human Genome Microarray 4x44K G4112F. The GSE75010 dataset was used to perform the differential expression analysis and the other analyses. The GSE48424 dataset was used as an independent validation set to test the reliability of our model.

Differential Expression Analysis

The mRNA profiles from the datasets were normalized by using Robust Multi-Array Average (RMA) method, and then the standardized data was obtained after log₂ transformation. In this study, the mRNA profile of GSE75010 was used to perform the differential expression analysis with the limma function package of R language version 3.5.2.¹⁶ The |log 2 (fold change [FC])|  >  1.5 and false discovery rate (FDR) ≤ 0.05 were considered as significant thresholds to screen differentially expressed genes (DEGs).

Functional Enrichment Analysis

Gene Ontology (GO) analysis (including biological process, molecular function and cellular component) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed based on the clusterProfiler function package of R language.¹⁷ And the DEGs were also subjected to fast gene set enrichment analysis (FGSEA), using fgsea package of R.¹⁸ P < 0.05 was considered as the threshold.

Proportion of Infiltrating Immune Cell

The relative proportions of 22 infiltrating immune cells in each sample of the GSE75010 dataset were evaluated using CIBERSORT software,¹⁹ which is a method to characterize the composition of immune cells with 547 preset barcode genes based on the deconvolution algorithm according to their gene expression profiles, and the sum of relative percentages of all estimated immune cell types in each sample is 1.

Protein–Protein Interaction (PPI) Networks

STRING database (https://string-db.org/,version 11.0) was used to analyze the functional connections and interactions among proteins,²⁰ and the PPI network was visualized by using Cytoscape (https://cytoscape.org/, version 3.7.2).²¹ Finally, the key genes of PPI network were screened by using the cytoHubba plug-in of Cytoscape software based on the Maximum neighborhood component (MCC) algorithm.

The Construction of Logistic Regression Model

In GSE75010, the multivariate logistic regression model was constructed by using the glm function package of R language,²² with the expression values of key genes as the continuous predictive variables and the sample type as the categorical responsive value (diseased or non-diseased samples). Then, the GSE48424 dataset was used as the validation set to determine the reliability of the logistic regression model via receiver operating characteristic (ROC) analysis.

Statistical Analysis

The difference in infiltrating immune cells among different groups was calculated by using Wilcoxon rank-sum tests, and P < 0.05 was considered as the significant threshold. Statistical analysis was performed with R software v3.5.2.

Results

The Identification of Differentially Expressed Genes in GSE75010

After standardization of the mRNA profiles in GSE75010, we found that there was no obvious deviation in overall mRNA level among all samples (Figure 1A), suggesting that the data could be used for subsequent analysis. In order to confirm the repeatability of the data within the group, principle component analysis (PCA) was performed and the result showed that the case group (samples with hypertension in pregnancy) and control group (healthy subjects) could be efficiently separated (Figure 1B), indicating a better reproducibility of the data. Then, the differential expression analysis in GSE75010 was performed and a total of 104 differentially expressed genes, including 76 upregulated genes and 28 downregulated genes, were identified between the samples with hypertension in pregnancy and healthy subjects (Figure 1C). Furthermore, the expressions of these 104 genes exhibited significant differences between case group and control group (Figure 1D). These results indicated that the 104 differentially expressed genes might be potentially associated with the development of hypertension in pregnancy.

Figure 1 The identification of differentially expressed genes in GSE75010 dataset. (A) The distribution of gene expression value in each sample after data standardization. The horizontal axis is the sample and the vertical axis is the gene expression value. (B) The PCA analysis. The points with different colors represent samples from different groups, and the closer the distance between two points, the more similar the expression of mRNA in the two samples. (C) The volcano plot of differentially expressed genes between case group and control group. The horizontal axis is Log2 FC, and the vertical axis is -log10 (FDR). The blue points represent upregulated genes, the red points represent downregulated genes. (D) The heatmap of differentially expressed genes between case group and control group. The horizontal axis represents samples, the vertical axis represents genes, the red color indicates high expression, and the blue color indicates low expression.

Functional Enrichment Analysis

To investigate the processes and metabolic pathways related to the progression of hypertension in pregnancy, the functional enrichment analysis was performed based on the 104 DEGs. The results indicated that there were 26 significantly enriched biological process (BP) terms including cell activation, regulation of p38 MAPK cascade and hormone secretion (P < 0.01), 13 significantly enriched molecular function (MF) terms including hormone activity, growth factor binding and monosaccharide binding (P < 0.05), as well as 7 significantly enriched KEGG pathways including cytokine–cytokine receptor interaction, HIF-1 signaling pathway and neuroactive ligand–receptor interaction (P < 0.05). The full list of significantly enriched GO terms and KEGG pathways is shown in Table S1. In addition, the top 20 most significantly enriched GO terms and related genes involved in the top 12 significantly enriched GO terms are shown in Figure 2A. The significantly enriched KEGG pathways are shown in Figure 2B. Additionally, based on the results of FGSEA, the DEGs in hypertension in pregnancy samples were significantly upregulated in Adipocytokine signaling pathway and other 21 pathways, significantly downregulated in Leukocyte transendothelial migration pathway and other 16 pathways (Figure 2C).

Figure 2 The functional enrichment and immune infiltration analyses in GSE75010 dataset. (A) The histogram of the top 20 most significantly enriched GO terms (the horizontal axis is the number of genes and the vertical axis is the name of GO terms), and the chord chart of genes enriched in the top 10 GO terms (the right semicircle represents the top 10 GO terms, and the left semicircle represents the genes enriched in the 10 GO terms). (B) The histogram of the top 7 most significantly enriched KEGG pathways. The horizontal axis is the number of genes and the vertical axis is the name of KEGG pathways. (C) The FGSEA results of DEGs in samples with hypertension in pregnancy. (D) The relative proportion of infiltrating immune cells in all samples with hypertension in pregnancy. (E) The violin plot of immune cells with significantly different infiltration proportions in samples with hypertension in pregnancy compared with healthy subjects. The horizontal axis is different groups, and the vertical axis is the relative infiltration proportion of immune cells. P value is determined by Wilcoxon test. (F) The PCA analysis of samples based on eight types of immune cells. The dots with different colors represent different types of samples.

It has been reported that p38 MAPK cascade is closely involved in inflammatory signaling and activated in response to various intracellular and extracellular stimuli, including oxidative stress, cytokines and growth factors.²³ Previous studies demonstrated that multiple pathways in the inflammatory microenvironment were used to adjust HIF-1alpha levels to affect macrophage biology, and HIF activation played a crucial role in mammalian immune cells through moderating inflammation.²⁴ In addition, the neuroactive ligand–receptor interaction was identified to participate in synaptic transmission and central nervous system (CNS) disorders, and the related genes could be induced by avian and human H5N1 infections.²⁵ Based on these reports, our results suggested that inflammatory response might be occurred during the progression of hypertension in pregnancy.

Immune Status of Patients with Hypertension in Pregnancy

Immune response occurs during the process of pregnancy. Compared with the non-pregnant women, the immune response of pregnant women is different, which also fluctuates during the process of pregnancy.²⁶ Increasing evidences indicate that immune disorders play important roles in the pathogenesis of hypertensive diseases during pregnancy.²⁷ Then we analyzed the differences of 22 infiltrating immune cell types between different samples from the GSE75010 dataset by combining LM22 eigenmatrix and CIBERSORT method.²⁸ The infiltration of immune cells in the 157 samples from GSE75010 is summarized in Figure 2D, suggesting that the changes in the proportion of infiltrating immune cells among different samples might represent the intrinsic characteristics of individual differences. In addition, eight types of immune cells including Dendritic cells activated, Macrophages M2, Mast cells resting, Neutrophils, NK cells resting, Plasma cells, T cells CD8 and T cells regulatory Tregs were significantly differentially infiltrated in samples with hypertension in pregnancy compared with healthy subjects (Figure 2E), which might account for the development of hypertension in pregnancy. Meanwhile, PCA analysis was performed and the result indicated that the samples could be efficiently divided into two categories (samples with hypertension in pregnancy and normal samples) according to the eight types of immune cells (Figure 2F), which further confirmed that the difference in the amount of immune cells might be the potential cause of hypertension in pregnancy.

The Construction of PPI Network

Next, STRING database was used to construct the PPI network based on the 104 differentially expressed genes, and the interaction pairs of proteins whose minimum required interaction score ≥0.4 were selected and used to reconstruct the PPI network. Then, the new PPI network was visualized by the Cytoscape software (Figure 3A). Then, the analysis of topology structure of the PPI network was performed by Cytoscape software, and each node in PPI network was scored based on the MCC algorithm. The top 20 genes selected according to the score are shown in Figure 3B, in which the darker the color, the more important the node in the network. The full list of these 20 genes is shown in Table S2. The results revealed that these 20 genes were significantly associated with the development of hypertension in pregnancy and might be potential biomarkers.

Figure 3 The construction of PPI network. (A) The PPI network was constructed based on the differentially expressed genes. Each dot represents a node. The more lines connected with the dot, the greater the degree of the node, the more important the gene in the network. (B) The top 20 genes with higher degree in PPI network were screened based on MCC algorithm. The darker the red, the higher the degree.

The Construction and Evaluation of Logistic Regression Model

Subsequently, the logistic regression model was constructed by using the samples from GSE75010 dataset based on the 20 key genes screened from PPI network, and we found that the P values of LEP, PRL and IGFBP1 were all less than 0.05 (Figure 4A), suggesting that the three genes contributed more to the model. Then, we re-constructed the logistic regression model with LEP, PRL and IGFBP1 and found that this final model conformed to the normality test (Figure S1A); meanwhile, there was a better linear relation between the independent variables (LEP, PRL and IGFBP1) and the dependent variables (sample type: diseased or non-diseased samples) (Figure 4B). In addition, there were no extreme points in the model that affected the accuracy of the model (Figure S1B). To confirm the accuracy of the model, the ROC curves of the model in two datasets were plotted and shown in Figure 4C. The area under curve (AUC) values in the training set (GSE75010) and validation set (GSE48424) were 0.8494 and 0.7346, respectively. These results suggested that the logistic regression model based on LEP, PRL and IGFBP1 could efficiently separate the patients with hypertension in pregnancy from the healthy subjects. Moreover, the model showed reliable value not only in the unattainable placenta samples but also in the accessible blood samples.

Figure 4 The construction of logistic regression model. (A) The volcano plot of 20 genes screened from PPI network in the model. P < 0.05 indicated that the gene contributed more to the model. (B) The component plus residual plot of LEP, PRL and IGFBP1 in the model. The obvious linear relationship between the horizontal axis and the vertical axis indicated that the independent variables were suitable to be included in the model. (C) The ROC curves of logistic regression model in GSE75010 and GSE48424 datasets. The horizontal axis indicates false positive rate (FPR), and the vertical axis indicates true positive rate (TPR). The AUC value is the area under the curve and can intuitively evaluate the quality of the model, the larger the AUC value, the better the model.

Discussion

According to the severity of disease, hypertension in pregnancy is divided into mild, moderate and severe hypertension in pregnancy, of which, severe pregnancy hypertension syndrome is a crucial feature of preeclampsia.²⁹ Placenta is the key organ involved in preeclampsia,³⁰ and the identification of differentially expressed genes between normotensive and preeclamptic placentas is crucial to understand the molecular mechanism of the disease, contributing to the development of efficient biomarkers or therapeutic targets. In this study, we analyzed the gene expression profiles of placenta samples between patients with hypertension in pregnancy and normal controls, and identified 104 differentially expressed genes, suggesting that these differentially expressed genes might be potentially associated with the occurrence and development of hypertension in pregnancy.

Increasing evidences have indicated that preeclampsia is a severe gestation disorder which is characterized by significant inflammation in the maternal body that affects the work of the monocyte-macrophage system.³¹ In the last decades, many risk genes regulated the process of inflammatory response have been found to affect the progression of hypertension in pregnancy. These findings may provide potential biomarkers for hypertension in pregnancy. For instance, the expression of nectin-4 was significantly elevated in preeclamptic placentas relative to uncomplicated pregnancies, and natural killer (NK) cells exhibited a potent cytotoxic effect against nectin-4 overexpressing trophoblast cells.³² Vinnars et al found that an inflammatory cytokine NKG2D was upregulated and combined with a dysregulated regulatory response in the placenta of preeclampsia.³³ Chen et al indicated that the serum level of cytokine IL-10 in preeclamptic patients was significantly lower than that in normal pregnant women, suggesting the contribution of regulatory T cells to immune tolerance in preeclampsia.³⁴ In addition, Chatterjee et al found that co-treatment with interleukin 4 and interleukin 10 could regulate immune cells and then prevent hypertension in pregnant mice.³⁵ These reports confirm that inflammatory response plays an important role in the progression of hypertension in pregnancy. Here, the GO term and KEGG pathway enrichment analyses were performed based on 104 differentially expressed genes, and the results showed several inflammation-related functions or pathways such as cytokine–cytokine receptor interaction, neuroactive ligand–receptor interaction and regulation of p38 MAPK cascade were enriched. Meanwhile, the relative proportions of eight types of infiltrating immune cells, including Dendritic cells activated, Macrophages M2, Mast cells resting, Neutrophils, NK cells resting, Plasma cells, T cells CD8 and T cells regulatory Tregs, exhibited significant differences in samples with hypertension in pregnancy compared with healthy subjects. All these data suggest that a strong inflammatory reaction may occur during the development of hypertension in pregnancy.

PPI network is a viable tool to understand cell functions, disease machinery, and drug design/repositioning.³⁶ In the present study, 20 key genes with higher scores in the PPI network were screened to construct the logistic regression model. To construct the model with stronger explanation, we optimized three genes LEP, PRL and IGFBP1 with P value less than 0.05. Leptin (LEP) has been demonstrated to induce hypertension through acting on transient receptor potential Melastatin 7 channel in the carotid body.³⁷ Hypomethylation of LEP gene was identified in early-onset preeclampsia, and the expression level of LEP was elevated.³⁸ The protein product leptin of LEP was obviously increased in the serum of patients with preeclampsia in comparison to those with normotensive pregnancy.³⁹ It was found that the hypomethylation of PRL could contribute to increased expression of prolactin in prolactinomas cases,⁴⁰ however, to our knowledge, PRL hypomethylation has not been reported in patients with hypertension in pregnancy. Prolactin, encoded by PRL, is a 23 kDa polypeptide hormone mainly secreted by the pituitary gland, and has multiple functions, such as regulating lactogenesis, balancing salt and water, modulation of cell proliferation, etc.⁴¹ Compared with normal controls, the prolactin level in preeclampsia patients was obviously elevated,^42–44 implying that prolactin was associated with the hypertension in pregnancy. As far as we know, the hypomethylated IGFBP1 has not been investigated in patients with hypertension in pregnancy. Fang et al demonstrated that the genetic variant of IGFBP1 was closely associated with an increased risk of hypertension, and lower IGF-binding protein 1 (encoded by IGFBP1) level in serum might predict an increased risk of hypertension.⁴⁵ IGF-binding protein 3 (encoded by IGFBP3) belongs to the family of IGF binding proteins and is elevated during the process of pregnancy, with suppressive effects on the growth of cells.⁴⁶ Compared with healthy controls, the serum level of IGF-binding protein 3 was obviously reduced in the patients with preeclampsia.⁴⁷ These studies indicate that the three genes are potentially related to the development of hypertension in pregnancy. Subsequently, the logistic regression model was constructed by bringing into the three genes and could efficiently separate patients with hypertension in pregnancy from healthy subjects for both placenta samples and blood samples.

Although our study suggests that the three genes LEP, PRL and IGFBP1 may be potential diagnostic signatures for hypertension in pregnancy, their specific functions need to be explored in the future.

Conclusion

In summary, based on the publicly accessible data in GEO database, we have identified 3 crucial genes, including LEP, PRL and IGFBP1, which probably play vital roles in hypertension in pregnancy. Moreover, a prognostic signature based on these genes is constructed.