CeRNA network analysis and functional enrichment of salt sensitivity of blood pressure by weighted-gene co-expression analysis

Introduction

Salt sensitivity of blood pressure (SSBP), first proposed by Kawasaki in 1978 (Kawasaki et al., 1978), is defined as the blood pressure of some people exhibiting changes parallel to changes in salt intake (Elijovich et al., 2016). Hypertension that related with SSBP is called salt sensitive hypertension (SSH) (Mo & Ren, 2012). SSBP is one of the early damage markers of hypertension and an independent risk factor for cardiovascular morbidity and mortality (Franco & Oparil, 2006), especially for Asian populations (Kario et al., 2018). Understanding the key regulatory factors and pathogenic pathways of SSBP can help us screen the high-risk individuals and conduct personalized treatments.

Evidences have showed that microRNAs (miRNAs) (Moazed, 2009) are involving in the pathogenesis of SSH. For example, miR-429 (Zhu et al., 2017), miR-133a (Guo et al., 2014) and miR-29b (Liu et al., 2010) are differentially expressed in SSH rat model, and hsa-miR-361-5p is associated with SSH patients (Qi et al., 2017). However, there lacks evidences about the associations between long non-coding RNAs (lncRNAs) (Wang & Chang, 2011) and SSBP. A study revealed that the expression level of lncRNA sONE is significantly decreased in rat model of SSH (Zhang et al., 2015). It is necessary to explore the expression profile of lncRNAs in SSBP individuals. In addition, with the rising advent of the competitive endogenous RNA (ceRNA) theory, it has been discovered that lncRNAs and miRNAs function together to form biological modules (Xu et al., 2011). LncRNAs could competitively combine with the miRNA response element (MRE) and repress miRNA’s negative regulation of target mRNAs. This lncRNA-miRNA-mRNA triple network has been demonstrated in several diseases, but not in SSBP (Chen et al., 2018; Pan et al., 2017b; Xie et al., 2017; Yuan et al., 2017; Zhang et al., 2016). Thus, it is important to elucidate the ceRNA co-regulatory of SSBP with bioinformatic analysis and data mining, which could construct the systematic regulatory network and explore the interactive relationship of different biomarkers.

Weighted correlation network analysis (WGCNA) is a systematic biology-based approach that transforms gene expression profiles into co-expressed gene modules, which provides insight into biological pathways (Langfelder & Horvath, 2008). WGCNA has been widely applied in biological areas such as acute myeloid leukemia (Pan et al., 2017a), hepatocellular carcinoma (Pan et al., 2016), schizophrenia (Ren et al., 2015), as well as in coronary artery disease (Yan, 2018). Similarly, gene set enrichment analysis (GSEA) is a gene expression analysis method that creates a priori gene sets according to prior knowledge about biological states and can be used to gain biological insight into the transcript level (Subramanian et al., 2005).

Taken together, this study aims to screen for SSBP associated differentially expressed hub lncRNAs and mRNAs using WGCNA method, construct the co-expression network and ceRNA network, and provide insight into the mechanisms of SSBP on transcriptome level.

Material and Methods

Results

Differentially expressed lncRNAs and mRNAs

SBC ceRNA array v1.0 contained 68,423 lncRNAs and 18,853 mRNAs. Figure 1 describes the screening procedures of lncRNAs and mRNAs. According to the four inclusion criteria, we discovered 274 DE-lncRNAs (119 upregulated and 155 downregulated) and 36 DEGs (23 upregulated13 downregulated) (Supplementary Fig. 1).

LncRNA-mRNA co-expression network construction

WGCNA was applied to detect the potential interactions between lncRNAs and mRNAs. Sample dendrogram analysis showed no outliers, so all 20 samples were included for further analysis. We could know from Fig. 2A that SS subjects were related with higher TC, HDL and lower TG. We used β = 6 to get adjacency matrix and generated a TOM dendrogram after calculating scale independence and mean connectivity (Fig. 2B). The cluster dendrogram shows three colors of modules (blue, turquoise and grey). The blue module contained 46 lncRNAs and 9 mRNAs and the turquoise module contained 191 lncRNAs and 17 mRNAs. Grey module was unclassified genes. Thus only blue and turquoise modules were considered for further analyses. The top 10 weighted lncRNAs and mRNAs in blue and turquoise co-expression network modules were summarized in the Supplemental Information (Tables S2 and S3). The correlations between modules and clinical traits are depicted in Fig. 2C. We found that the blue module had a positive relationship with salt-sensitivity (R = 0.7, P < 0.01), HDL (R = 0.53, P = 0.02), and TC (R = 0.55, P = 0.01). The turquoise module was positively related with TG (R = 0.8, P < 0.01) and LDL (R = 0.54, P = 0.01). A cluster heatmap was performed to comprehensively reflect the adjacency and topological structures of genes (Fig. 2D). We found two clusters in the heatmap, which represent the two modules. Then, the co-expression network was visualized with Cytoscape software (Figs. 3A and 3B). Due to too many nodes, only top ten associated lncRNAs of each mRNA were performed in turquoise co-expression network (Fig. 3B).

CeRNA network construction

We ranked lncRNAs by module membership in blue and turquoise modules. LncRNAs with P < 10⁻⁷ were selected for target miRNAs and mRNA prediction using miRDB and miRmap databases. Finally, 17 of 46 (37.0%) lncRNAs from the blue module and 128 of 191 (67.0%) lncRNAs from the turquoise modules were chosen for further prediction analysis because their full sequences could be obtained from LNCipedia. After target genes prediction, 79 miRNAs and 220 mRNAs were included in the blue module, while 347 miRNAs and 607 mRNAs were included in the turquoise module. Then, we merged the mRNAs from prediction with co-expression and got 9 (MGAM2, ITGB4, KCNH7, CLEC1A, KIF6, LOC643802, DHRS13, IL1B, SCN9A) overlapping mRNAs in the blue module and 14 mRNAs (TARP, GPR68, CCR5, TUSC5, ZNF98, TMEM14B, FOSL1, PLG, FAM111B, DNAJB6, SEC14L2, PTGER3, GLIPR1L2, TKTL1) overlapping in the turquoise module. Due to PTGER3 and KCNH7 were associated with known SSBP mechanism pathways, such as sodium reabsorption and potassium ion-channel, which cause fluctuations of blood pressure (Kelly & He, 2012), these two DEGs were selected to perform ceRNA network of two modules (Figs. 4A and 4B). There were 10 miRNAs and 10 mRNAs in blue turquoise module, which forming 11 ceRNA loops. Similarly, there were 32 miRNAs and 51 lncRNAs in turquoise module, which forming 73 ceRNA loops.

Discussion

In this study, 274 DE-ncRNAs and 36 DEGs were identified through differential analysis. WGCNA analysis discovered the correlation between clinical traits and three modules (turquoise, blue and grey), and turquoise and blue modules were selected for enrichment analysis. Eighty-four ceRNA loops were identified and lnc-CD302-1:1 →hsa-miR-1283 →PTGER3 may be of great importance for participating in mechanism of SSBP.

GO and KEGG analyses revealed that the turquoise module was related with aldosterone synthesis and secretion, which was one of the pathogenesis pathways of SSBP (Dengel et al., 2001). Aldosterone is a kind of steroid hormone and mainly responsible for fluid homeostasis of the body. Excess aldosterone would result in the reabsorption of salt and water and thereby the elevation of blood pressure (Shibata & Fujita, 2011). The role of aldosterone is more obvious in SSH patients, who have already undergone damages of kidney and exhibit overproduction of aldosterone. The lncRNAs and mRNAs that we screened out may involve in renin-angiotensin-aldosterone (RAAS) pathway through regulating the expressions of corresponding genes.

Among the 14 mRNAs of turquoise ceRNA network, prostaglandin E receptor 3 (PTGER3) plays a prominent role in the cardiovascular system. They involved in biological activities such as the regulation of blood pressure, metabolism of sodium and water, inflammation, reangiostenosis, and so on (Ceddia et al., 2016; Grilo et al., 2011). Among the 32 miRNAs that had predictive relationships with PTGER3, hsa-miR-1283 was associated with essential hypertension. Yang et al. (2015) found that hsa-miR-1283 may regulate the expression of ATF1 through binding with the 3′UTR of the gene, which was abnormally expressed in essential hypertension. Due to SSBP is a clinical trait of hypertension, hsa-miR-1283 is probably associated with SSBP. Although there is no direct evidence showing that lnc-CD302-1:1 is directly associated with SSBP, the co-expression of lnc-CD302-1:1 and PTGER3 as well as the predictive relationship of lnc-CD302-1:1 to hsa-miR-1283 provide hypothesis for future investigation.

In blue module, PI3K-Akt signaling pathway was the most significant pathway for SSBP. Insulin could activate the PI3K-Akt pathway by stimulating nitric oxide (NO) synthesis and thereby promoting vasodilatation (Kobayashi et al., 2004). Insulin sensitivity is also related with salt sensitivity by hyperinsulinemia, over-activation of sympathetic nervous system, and reducing suppression of RAAS pathway (Yatabe et al., 2010). Thus, the active of PI3K-Akt pathway may play a protective role in prevention and treatment of SSBP, corresponding to the strong correlation of the blue module and SSBP.

GSEA showed that SSBP is mainly associated with basic biological processes including glycan biosynthesis, and inositol phosphate metabolism. However, the results of GO and KEGG that based on WGCNA screening showed more SSBP-relevant pathways. So, the process of screening the differentially expressed RNAs, WGCNA co-expression analysis, and target genes prediction could effectively discover the accurate pathway information for SSBP, and more importantly, explore the function of ceRNA theory in the pathogenesis of SSBP. As all the mRNAs are considered when conducting GSEA analysis, the results of GSEA could be more comprehensive and therefore helpful for discovering the unknown pathways. However, there are some limitations for the two methods. For WGCNA, the nodes with low connectivity degree would not be recognized by modules. The heterogeneities of samples would result in the difficulties in module recognition and sample size under 15 could not use WGCNA (Liu et al., 2017). For GSEA, although it is a group-set analysis method, there are a few false positives due to focusing on sample-level variations rather than pathway-levels (Yi & Stephens, 2008). Thus, new methods are necessary to be developed to overcome the limitations. For this study, the combination of WGCNA data mining with GSEA could screen for key regulators of diseases and provide ideas for exploring pathogenesis pathways of other diseases.

There are a few limitations in our study. The regulatory effects of lncRNA to miRNA and miRNA to mRNA need to be further validated by real-time polymerase chain reaction with a larger sample size and functional experiments. The representative of participants is limited because they are older age and living in Beijing. In addition, whether the microarray data mining process is applicable for other diseases is currently unknown and should be explored in future studies.

Conclusions

The microarray data mining process based on WGCNA co-expression analysis had identified 84 ceRNA loops for SSBP. The results of our study provide implications for understanding the pathogenesis of SSBP and facilitate the precise diagnosis and therapeutics.

Supplemental Information