Bioinformatics-based identification of miRNA-, lncRNA-, and mRNA-associated ceRNA networks and potential biomarkers for preeclampsia.

This study aimed to identify significantly altered long non-coding RNAs (lncRNAs), microRNAs (miRNAs), mRNAs, pathways in preeclampsia (PE), and to investigate their targeted relationships and biological functions.GSE96985 from Gene Expression Omnibus database was extracted, involving 3 PE and 4 normal tissues. After the differential expression analysis of miRNAs, lncRNAs, and mRNAs using the limma package, protein-protein interaction (PPI) network and module analyses were performed for differentially expressed mRNAs (dif-mRNAs). Combined with the miRanda and miRWalk tools, a regulatory relationship between dif-miRNAs and dif-mRNAs/lncRNAs (dif-mRNAs/dif-lncRNAs) was predicted. Finally, mRNA-miRNA-lncRNA regulatory network construction was performed using Cytoscape software.A total of 511 dif-mRNAs were screened in PE. The top 5 nodes in the PPI networks included up-regulated complement component 3 (C3), C-X-C motif chemokine ligand 8 (CXCL8), and fibronectin 1 (FN1). Three significant network modules were identified for dif-mRNAs. C3 and CXCL8 were identified in module A, and FN1 was identified in module C. A disintegrin and metalloproteinase with thrombospondin motifs 6 (ADAMTS6) was down-regulated by the miR-210-3p. Therefore, lnc-CTD-2383M3.1 functions as a competing endogenous RNA in ADAMTS6 expression regulation by competitively binding to miR-210-3p during the regulation process of PE.C3, CXCL8, FN1, and ADAMTS6 might be involved in the development of PE. The lnc-CTD-2383M3.1-miR-210-3p-ADAMTS6 axis might be a potential regulatory mechanism in PE.

Introduction

Preeclampsia (PE) is a common pregnancy-associated disorder and is one of the main causes of pregnancy-related maternal and neonatal morbidity with clinical characteristics of pregnancy-induced hypertension and proteinuria.[ Every year, approximately 76,000 women and 500,000 babies die from PE, which is one of the major causes of maternal mortality.[ PE is an unpredictable multifaceted syndrome that has been associated with a high risk of cardiovascular complications and diabetes mellitus. PE can rapidly develop into the life-threatening condition, eclampsia.[ However, few effective treatments are currently available for PE patients in clinic. The only available cure for PE is delivery of the fetus and placenta, which reduces proinflammatory substances in the maternal cardiovasculature.[ The cause of PE remains elusive and numerous research studies have been carried out to investigate it.[ With the emergence of clinical signs and symptoms of PE in the third trimester of pregnancy, identification of early diagnosis/prognosis bio-markers has become a relevant research focus.

In the past few years, numerous studies have been performed on the biomarkers involved in PE. Vascular endothelial growth factor has been considered a down-regulated gene that plays an essential role in the abnormality of placenta and vascular dysfunction of PE.[ Furthermore, the elements involved in the activity of placental unfolded protein response, including activating transcription factor-6, phosphorylated inositol-requiring transmembrane kinase/endoribonuclease 1α, transcription factor X-box-binding protein-1, and glucose-regulated proteins GRP78 and GRP94, have a higher expression in PE patients compared to normotensive controls.[ Previous studies have shown that microRNAs (miRNAs) had a crucial role in PE. MiR-141 is one of the most abundant miRNAs in placental plasma during pregnancy.[ Recently, miR-141 has been found to have a crucial function during intercellular communication between the fetal trophoblast and maternal immune cells. It can be used as a marker of potential diseases like PE during pregnancy.[ Non-coding RNAs account for 90% of the human genome and are the largest group of RNA transcripts. In particular, long non-coding RNAs (lncRNAs) can potentially sponge their target miRNAs and mediate post-transcriptional regulation of target mRNAs by impairing target miRNA ability.[

In the present study, a microarray dataset for PE was downloaded from a public database. Then, differentially expressed mRNAs (dif-mRNAs), differentially expressed miRNAs (dif-miRNAs), and differentially expressed lncRNAs (dif-lncRNAs) between PE and control tissues were investigated using enrichment analysis. To further explore the mechanism of key mRNAs, lncRNAs, and miRNAs in PE, protein-protein interaction (PPI) network and regulatory network analyses were performed. This study might provide novel therapeutic targets for PE.

Materials and methods

Data source

Microarray dataset GSE96985 (miRNA sequencing platform: GPL20712 Agilent-070156 Human miRNA [miRNA version], lncRNA and mRNA sequencing platform: GPL22120 Agilent-078298 human competing endogenous RNA (ceRNA) array V1.0 4X180K [Probe Name Version]) from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database was downloaded. The database included 3 PE patients and 4 normal tissues. GSE96985 is a public dataset. Therefore neither ethics committee approval nor patient informed consent was required for analyzing GSE96985 data.

Differential expression analysis

The original data were preprocessed using the limma package.[ Data preprocessing included background correction, normalization, and concentration prediction. Gene expression values were calculated by mapping probes to symbols utilizing the microarray dataset and annotation files from the chip platform. The average value was defined as the value of gene expression when multiple probes were matched to one symbol.

The dif-mRNAs, dif-miRNAs, and dif-lncRNAs between PE patients and normal tissues were analyzed using a typical Bayesian method and corrected by the Benjamini/Hochberg tests. The heat maps for dif-mRNAs, dif-miRNAs, and dif-lncRNAs were plotted using the R package in pheatmap software (version: 1.0.10, https://cran.r-project.org/web/packages/pheatmap/index.html).[ The dif-miRNAs were defined as miRNAs with |log fold change (FC)| > 1 and P value <.05. The dif-lncRNAs and dif-mRNAs were defined as lncRNAs or mRNAs with |log2 (FC)| > 2 and P values <.05.

PPI network and module analyses for dif-mRNA

Using the STRING database (version: 10.0, http://www.string-db.org/),[ interactions between dif-mRNAs were retrieved. A PPI score = 0.4 was defined as the median confidence. Cytoscape software[ (version: 3.2.0, http://www.cytoscape.org) was used to construct the PPI network based on obtained interactions. Significant network modules were analyzed using Molecular Complex Detection (MCODE) plugin (version: 1.4.2, http://apps.cytoscape.org/apps/MCODE)[ in Cytoscape software. Modules with scores ≥ 8 were considered to be significant.

Pathway enrichment analysis

Based on the clusterprofiler tool in R package[ (version: 2.4.3, https://bioconductor.org/packages/3.2/bioc/html/clusterProfilter.html), the Kyoto Encyclopedia of Genes and Genomes (KEGG)[ pathway enrichment analysis was performed for dif-mRNAs, mRNAs in significant modules, and those in the lncRNA-miRNA-mRNA network. Related genes involved in the ceRNA net were also evaluated using the KEGG pathway enrichment analyses. Significantly enriched pathways with P values <.05 and involving no less than two mRNAs were selected.

Regulatory network analysis of miRNA with mRNA/lncRNA

The dif-miRNA and dif-lncRNA pairs were predicted using a score ≥150 and energy < −20 based on the miRanda tool (version: 3.3a, https://omictools.com/miranda-tool).[

The miRWalk 2.0[ (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) was used to predict the targeted mRNAs for dif-miRNAs in the above relationship. The mRNA-dif-miRNA pairs were selected according to the condition that the pairs appeared at least 5 times among the following databases: miRWalk, miRanda, miRDB, miRMap, Pictar2, RNA22, and Targetscan. The mRNAs-lncRNA pairs were finally selected among the acquired regulation pairs of mRNAs-lncRNAs with the condition that mRNAs belong to the dif-mRNA group.

LncRNA-mRNA coexpression network analysis

The dif-mRNAs among dif-mRNA-dif-miRNA pairs and dif-lncRNAs among dif-lncRNA-dif-miRNA pairs were selected to construct relative pairs of dif-mRNA-dif-lncRNAs. The Pearson correlation coefficient for dif-lncRNA-dif-mRNAs was calculated. The threshold values for significantly correlated dif-lncRNA-dif-mRNA pairs were R-value for top 100 with a false discovery rate <0.05.

LncRNA-miRNA-mRNA integration analysis

The Cytoscape tool was used to construct a complex network of lncRNA-miRNA-mRNAs by integrating the obtained lncRNA-mRNA, miRNA-mRNA, and miRNA-lncRNA pairs.

Results

According to the screening criteria, there were 511 dif-mRNAs (364 up-regulated and 147 down-regulated), 251 dif-lncRNAs (93 up-regulated and 158 down-regulated), and 36 dif-miRNAs (14 up-regulated and 22 down-regulated) in PE tissues compared with normal tissues.

The bidirectional clustering heat maps for dif-mRNAs, dif-lncRNAs, and dif-miRNAs are shown in Figure 1. The difference analysis can significantly separate the 3 PE patients from the 4 normal tissues on the level of mRNA, lncRNA, and miRNA. Different analysis results were reliable.

Figure 1

Bidirectional clustering heat map for mRNAs, dif-lncRNAs, and dif-miRNAs. (A) Heat map of differentially expressed mRNAs from PE and normal tissues (|logfold change (FC)| > 1, adjusted P < .05). (B) Heat map of differentially expressed miRNAs from PE and normal tissues (|logfold change (FC)| > 2, adjusted P < .05). (C) Heat map of differentially expressed lncRNAs from PE and normal tissues (|logfold change (FC)| > 2, adjusted P < .05).

Enrichment analysis

Enrichment analysis for the up- and down-regulated mRNAs was conducted. The up-regulated mRNAs were implicated in 36 pathways. The top ten pathways are shown in Figure 2A. Down-regulated mRNAs were implicated in 6 pathways (Fig. 2B).

Figure 2

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for differentially expressed mRNAs (dif-mRNAs). (A) Top ten up-regulated dif-mRNAs. (B) Down-regulated dif-mRNAs.

PPI network and module analyses

The PPI network included 229 nodes and 605 interactions (Fig. 3). Three significant models were separately identified from the PPI network (Fig. 3). Model A had 17 nodes and 136 interactions, model B had 10 nodes and 45 interactions, and model C had 9 nodes and 36 interactions. Complement component 3 (C3), C-X-C motif chemokine ligand 8 (CXCL8), and fibronectin 1 (FN1) were identified in the modules and among the top 5 nodes (C3 and CXCL8 were identified in module A and FN1 was identified in module C). The KEGG enrichment analysis was conducted for the mRNAs in significant modules. The top 10 results are shown in Figure 4.

Figure 3

Protein-protein interaction (PPI) network for differentially expressed mRNAs and modules A, B, and C identified from PPI network for dif-mRNAs. Red circles represent up-regulated mRNAs, green diamonds represent down-regulated mRNAs. Node size and degree are proportional.

Figure 4

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for differentially expressed mRNAs involved in modules (top ten listed).

Regulatory predictive analysis

The miRanda tool was used to predict the regulatory relationship for all dif-miRNAs and dif-lncRNAs. A total of 71 lncRNA-miRNA interaction pairs were obtained, including 24 dif-miRNAs and 45 dif-lncRNAs. To predict the targeted mRNAs for the above dif-miRNA, the miRWalk2.0 database was utilized. A total of 439 miRNA-mRNA interaction pairs were identified, including 22 dif-miRNAs and 186 dif-mRNAs.

Furthermore, the 22 dif-miRNA-regulated mRNAs and lncRNAs were screened for co-expression. The top 100 positively correlated lncRNA-mRNA pairs were obtained based on the screening threshold.

Regulatory network analysis

LncRNA-miRNA-mRNA pairs were further sorted according to the information available for all the pairs. A total of 605 interaction pairs were obtained, including 100 lncRNA-mRNA co-expressed, 439 miRNA-mRNA interaction, and 65 lncRNA-miRNA interaction pairs. The lncRNA-miRNA-mRNA network contained 9 up-regulated miRNAs, 13 down-regulated miRNAs, 137 up-regulated mRNAs, 49 down-regulated mRNAs, 16 up-regulated lncRNAs, and 11 down-regulated lncRNAs (Fig. 5). The lnc-CTD-2383M3.1 might function as a ceRNA in regulating disintegrin and metalloproteinase with thrombospondin motifs 6 (ADAMTS6) expression in PE by competitively binding to miR-210-3p (Fig. 5). The KEGG pathway enrichment analysis was carried out for all dif-mRNAs in the network (Fig. 6).

Figure 5

The mRNA-miRNA-lncRNA regulatory network. Red circles and green diamonds represent up- and down-regulated mRNAs, respectively. Yellow hexagons and blue quadrilaterals represent up- and down-regulated lncRNAs, respectively. Red triangles and dark green arrows represent up- and down-regulated miRNAs, respectively.

Figure 6

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for differentially expressed mRNAs involved in the mRNA-miRNA-lncRNA regulatory network.

Discussion

In the present study, 511 dif-mRNAs (364 up-regulated and 147 down-regulated), 251 dif-lncRNAs (93 up-regulated and 158 down-regulated), and 36 dif-miRNAs (14 up-regulated and 22 down-regulated) in PE tissues were screened. In the dif-mRNA PPI network, 3 (module A, B, and C) significant network modules were identified for dif-mRNAs. In addition, C3, CXCL8, and FN1 were identified in the modules and among the top 5 nodes (C3 and CXCL8 were identified in module A, FN1 was identified in module C). In the regulatory network, miR-210-3p targeted ADAMTS6 and showed a negative correlation. The lnc-CTD-2383M3.1 was predicted to be involved in the ceRNA network of ADAMTS6 and miR-210-3p.

Pregnancy is an immunological process in which the immune system should allow the fetus to thrive and protect both the mother and fetus from pathogens. Even though PE is a vascular disease, immunological mechanisms have been implicated.[ Previous studies have suggested that inadequate control of the maternal complement system contributes to the etiology of PE.[ The complement system belongs to the immune system, which can be activated by lectin. Studies have shown that complement C3 is the central component of all activation pathways.[ Indeed, the abundance of C3 in human serum indicated that C3 is a key component of immunity in the process of discriminating between self and non-self.[ Studies have shown that in complement inhibitor-deficient mice, activation of the complement alternative pathway is the key mechanism for reproductive failure.[ Recently, it has been shown that an alternative complement pathway is activated in human pregnancies with severe PE.[ Eculizumab is a targeted inhibitor of complement protein C5 that was able to successfully treat a patient suffering from the HELLP syndrome, which is a life-threatening complication of PE. This demonstrated that the complement system could provide a promising target for drug development in severe PE.[ Due to regulation of protective haplotype, the extravillous trophoblasts (EVTs) invading the maternal tissue during placentation do not encounter a vigorous complement attack of C3 activation.[ Concurrently, with an increased levels of factor B in early pregnancy, the predisposing haplotype might result in an increased level of complement activation. This indicates that a complement attack of the EVT invasion would consequently damage the placental function, resulting in an increased occurrence of severe PE.[ In the present study, C3 was significantly up-regulated in PE. The KEGG pathway enrichment analysis showed that C3 was involved in the complement and coagulation cascade pathway. Thus, C3 might play a significant role in PE.

Studies have shown that PE pathogenesis is involved in insufficient trophoblast invasion, causing dysfunction in the endothelial system and an increased inflammatory response.[ In PE cases, insufficient invasion of the trophoblast leads to insufficient remodeling of the spiral artery, which might result in hypoxia-reperfusion injury and oxidative stress in the placenta.[ CXCL8, a C-X-C motif chemokine ligand, is present in the decidua and trophoblast, which also express CXCL8 receptors.[ CXCL8-stimulated migration and invasion of the trophoblast cells is aided by the increasing expression of matrix metalloproteinase (MMP)-2 and -9.[ FN might contribute to the maturation of the trophoblast cells, which is an important component of trophoblastic basement membranes.[ Studies have shown that an inadequate trophoblast invasion might induce the pathogenesis of PE. FN1 is primarily synthesized by endothelial cells and has been reported to be involved in regulatory functions of cell interactions. FN1 circulating values correspond to endothelial damage. In addition, PE is considered to be a pathology of the endothelial cells.[ FN tissue remodeling is an important process during embryonic development and wound healing. In the present study, CXCL8 and FN1 were significantly up-regulated in PE. The KEGG pathway enrichment analysis showed that CXCL8 was involved in the pathway of viral protein interaction with the cytokine and cytokine receptor FN1 was involved in the PI3K-Akt signaling pathway. Therefore, CXCL8 and FN1 might be involved in the development of PE.

The present study also found that ADAMTS6 expression was down-regulated by the miR-210-3p. Lnc-CTD-2383M3.1 might serve as a ceRNA in this process. ADAMTS genes are members of MMP, which play an important role in degradation and repair processes of extracellular matrix (ECM). ADAMTS roles in PE might include invasion regulation of the spiral artery and ECM arrangement of the placenta.[ An increasing number of studies have found that miR-210-3p is mainly expressed in the villous and tunica adventitia trophoblastic layers.[ MiR-210-3p is involved in the process of placenta. Through the mitogen-activated protein kinase-dependent mechanism, elevated levels of miR-210-3p decreased the invasion of EVT.[ Hence, it was suggested that lnc-CTD-2383M3.1 acted as a ceRNA in regulating ADAMTS6 expression by competitively binding to miR-210-3p in the progression of PE.

Conclusion

In conclusion, a series of bioinformatics analyses was conducted in PE patients. The results revealed that C3, CXCL8, and FN1 might be implicated in the pathogenesis of PE. In addition, lnc-CTD-2383M3.1 might function as a ceRNA in ADAMTS6 expression regulation in PE by competitively binding to miR-210-3p. These results must be validated and supported by further experimental studies.

Author contributions

Conceptualization: Ruo-an Jiang.

Data curation: Feifei Zhou.

Formal analysis: Ruo-an Jiang, Ting Wang, Jing He.

Investigation: Dong Xu.

Methodology: Ruo-an Jiang, Ting Wang, Yingsha Yao.

Resources: Feifei Zhou.

Supervision: Dong Xu.

Writing – original draft: Ruo-an Jiang.

Writing – review & editing: Ruo-an Jiang, Ting Wang, Feifei Zhou, Yingsha Yao, Jing He, Dong Xu.