Integrative analysis of ceRNA network and DNA methylation associated with gene expression in malignant pheochromocytomas: a study based on The Cancer Genome Atlas.

BACKGROUND: Competitive endogenous RNAs (ceRNAs) have revealed a new mechanism of interaction between RNAs. Epigenetic regulation in the gene expression dynamics has become increasingly important in malignant pheochromocytomas (PCCs). We performed an integrative analysis of ceRNA networks and DNA methylation to identify key biomarkers and contribute to the understanding of the molecular biological mechanisms of malignant PCCs.
METHODS: Differentially expressed genes in malignant PCCs and controls were identified from The Cancer Genome Atlas database by using the Limma package in R (v3.4.4). An abnormal lncRNA-miRNA-mRNA ceRNA network was constructed for malignant PCCs, and function enrichment analysis was performed using the Database for Annotation, Visualization, and Integrated Discovery. For DNA methylation datasets, the methylation analysis package was used in identifying differential methylation genes, and potential prognostic genes were identified by Kaplan-Meier survival analysis.
RESULTS: A total of 447 lncRNAs, 26 miRNAs, and 1,607 mRNAs were found to be differentially expressed in malignant PCCs as compared with those in normal samples. We then constructed an abnormal lncRNA-miRNA-mRNA ceRNA network for malignant PCCs. The network consisted of 12 lncRNAs, 6 miRNAs, and 220 mRNAs. Functional enrichment analysis showed that differentially expressed mRNAs were particularly enriched in the biological process, cellular component, and molecular function. Furthermore, four differentially expressed mRNAs from ceRNAs were identified through the cross-analysis of gene expression and DNA methylation profiles. LncRNA C9orf147 and 6 out of 220 mRNAs were indicated as prognostic biomarkers for patients with malignant PCCs (P<0.05).
CONCLUSIONS: Our research increases the understanding of the pathogenesis of malignant PCCs and offers potential genes as underlying therapeutic targets or prognostic biomarkers. 2020 Translational Andrology and Urology. All rights reserved.

Introduction

Primarily arising from the chromaffin cells of the adrenal medulla (80–85%) or extra-adrenal paraganglia (15–20%) (1), pheochromocytomas (PCCs) are rare catecholamine-secreting tumors with a highly variable clinical presentation, including hypertension, headache, sweating, and palpitation (2,3). Approximately 10% of patients with PCC are malignant, which are defined by distant metastases, and the differential diagnosis between malignant and benign tumors remains challenging. Moreover, the risk of malignancy greatly exceeds the classical 10% in patients with extra-adrenal disease—the so called paragangliomas (1,4). Individuals with malignant PCCs have unfavorable 5 year survival rates that vary between 20% and 50% (5).

In the past decade, major advances in the understanding of the pathogenesis and methods for diagnosing malignant PCCs have been reported. In messenger RNA (mRNA) expression array data, contactin 4 (CNTN4) is more frequently expressed in malignant PCCs than in benign PCCs, and CNTN4 expression is consistently associated with malignant behavior in PCCs (6). In addition, differentially expressed microRNAs (miRNAs) in benign and malignant PCCs have been characterized and found to be markers of malignancy or disease recurrence (7). Recently, long noncoding RNAs (lncRNAs) hold great promise as novel biomarkers and therapeutic targets for various types of cancer because of their genome-wide expression patterns and tissue-specific expression characteristics (8). LncRNA-mediated gene expression involves various mechanisms, such as transcription regulation, translation, protein modification, and cell signaling pathways. However, few studies have directly explored the mechanisms of lncRNAs in malignant PCCs.

The abnormal expression of competitive endogenous RNAs (ceRNAs) is closely related to the occurrence, development, and prognosis of tumors (9). CeRNAs refer to transcripts, such as mRNA, tRNA, rRNA, lncRNA, pseudogene RNA, and circular RNA (10). CeRNAs act as molecular sponges for miRNA through miRNA response elements (MREs) and form complex regulatory networks (11). As ceRNAs, lncRNAs serve as miRNA sponges, playing important roles in the regulation of cancer-related genes. For example, lnc-MD1, a muscle-specific long noncoding RNA, governs the time of muscle differentiation by acting as a ceRNA in mouse and human myoblasts. Lnc-MD1 sponges miR-133 to regulate the expression of MAML1 and MEF2C, which are transcription factors that activate muscle-specific gene expression (12). Moreover, lnc-HULC acts as an endogenous sponge that downregulates a series of miRNA activities, including miR-372. MiR-372 inhibition reduces the translational repression of its target gene, PRKACB, which in turn induces CREB phosphorylation (13). In addition, the association of the RNA-binding protein HuR with lnc-p21 favors the recruitment of let-7/Ago2 to lnc-p21, and thus lnc-p21 stability is low (14). DNA methylation plays an important role in the regulation of gene expression through an epigenetic control in various biological processes and diseases (15). In recent years, the DNA methylation of specific CpG islands has been detected in many types of malignant tumors, and DNA methylation aberrancies are considered hallmarks of human cancers (16). Therefore, an integrative analysis of ceRNAs and DNA methylation associated with gene expression should be conducted to provide novel clues to the pathogenesis of malignant PCCs.

In the present study, we first investigated potential differentially expressed lncRNAs, miRNAs, and mRNAs in malignant PCCs on the basis of The Cancer Genome Atlas (TCGA) database. Subsequently, we constructed a differentially expressed lncRNA-miRNA-mRNA ceRNA network and performed methylation analysis and overall survival analysis on the components of the network to identify malignant PCC-related prognostic biomarkers. This study provided potential lncRNA, miRNA, and mRNA biomarkers, which are expected to contribute to the early diagnosis, treatment, and prognosis of malignant PCCs.

Methods

Analysis of differentially expressed genes on TCGA

The relevant data provided by TCGA are publicly available and open ended and do not require the approval of the local ethics committee (17). All data were derived from the TCGA database (https://tcga-data.nci.nih.gov/), including clinical sample information, mRNA-Seq, miRNA-seq, and lncRNA-seq data, from malignant PCCs. All data were acquired using the Illumina HiSeq 2000 RNA-Seq platform. Sequencing data were converted through Homo_sapiens. GRCh38.94.chr.gtf, which was downloaded from the Ensembl database (https://asia.ensembl.org/index.html). Basing on the Empirical Analysis of Digital Gene Expression Data package (edgeR) and the “limma” package in R, we normalized and analyzed the downloaded data to obtain differentially expressed mRNA, miRNA, and lncRNA molecules. The cut-off value of |log2 fold Change (FC)| was >2.0, and the adjusted P value was <0.01.

Prediction of lncRNA-miRNA and miRNA-mRNA interactions

In this study, we used miRcode (http://www.mircode.org/) to analyze the interactions between lncRNAs and miRNAs on the basis of the differential expression of malignant PCC-specific miRNAs. We used miRTarBase (http://mirtarbase.mbc.nctu.edu.tw), miRDB (http://www.mirdb.org/miRDB), and TargetScan (http://www.targetscan.org) online analysis tools to predict miRNA-targeted mRNAs. The mRNAs that appeared in at least 2 of 3 online prediction tools were selected to overlap with the differentially expressed mRNAs in malignant PCCs. The intersection between two groups were treated as the target mRNA.

Construction of the lncRNA/mRNA/miRNA ceRNA network

Basing on the differentially expressed lncRNAs, miRNAs, and mRNAs and the lncRNA-miRNA and miRNA-mRNA pairs, we used Cytoscape v3.6.1 to establish the lncRNA-miRNA-mRNA ceRNA network. Then, we employed the Database for Annotation, Visualization, and Integrated Discovery to obtain information for Gene Ontology (GO) biological processes. We used a Cytoscape plug-in (BNGO) tool, which involves the following three parts: biological process (BP), cellular component (CC), and molecular function (MF) to further understand the function of these differentially expressed genes in malignant PCCs. The selection criteria were set as an adjusted P<0.05 after FDR correction.

Methylation analysis of differentially expressed genes

DNA methylation data were downloaded from TCGA database and estimated by the beta value. As previously mentioned, the limma package in R was used to obtain differential methylation genes (DMGs) with a cut-off value of |log2FC| >1.0 and P value of <0.05 (18). An overlap was conducted between DMGs and differentially expressed genes (lncRNA, miRNA, mRNA) involved in ceRNA network. The Pearson correlation coefficient was calculated between the expression and methylation level of overlapping genes. Moreover, to further explore the relationship between gene expression and methylation, we searched the methylation sites that were negatively correlated with gene expression.

Relationship of differentially expressed genes with prognosis

The clinical data of patients, including prognosis and survival time, were obtained from TCGA. Survival packages in R and Kaplan-Meier curves were used to assess the relationship between the expression of differentially expressed lncRNAs, miRNAs, and mRNAs and overall survival in patients with malignant PCCs. A P value of <0.05 was considered statistically significant.

Results

Differential RNA expression analysis based on TCGA

We extracted and analyzed expression profiles of lncRNA, miRNA, and mRNA between 38 malignant PCCs and three control samples on TCGA. Heat maps on all linkage clustering of differentially expressed RNAs along with volcano plots were constructed by R (). In total, 447 differentially expressed lncRNAs and 26 differentially expressed miRNAs were screened out. Among the differentially expressed lncRNAs, 286 were downregulated, and 161 were upregulated. Differentially expressed miRNAs contained 21 downregulated miRNAs and five upregulated miRNAs. Moreover, 1,607 differentially expressed mRNAs were found, of which 660 were downregulated and 947 were upregulated.

Figure 1

Heatmaps and volcano maps of the differential expression of genes in malignant pheochromocytomas (PCCs) between 38 tumor tissues and 3 normal tissues. The color from green to red shows a trend from low expression to high expression. (A) lncRNAs; (B) miRNAs; (C) mRNAs.

LncRNA-miRNA and miRNA-mRNA interactions

Using the miRcode tool, we compared differentially expressed lncRNA-miRNA pairs and compared 12 differentially expressed lncRNAs and six differentially expressed miRNAs (). Subsequently, the target genes of the six miRNAs were predicted by miRTarBase, miRDB, and TargetScan online tools. Finally, 2,363 target mRNAs were identified. Furthermore, to obtain a differentially expressed miRNA-mRNA interaction pair in malignant PCCs, we matched the predicted target gene with differentially expressed mRNAs (). MiRNA-targeted mRNAs were excluded if they were not included in differentially expressed mRNAs. A total of 26 differentially expressed lncRNA-miRNA pairs and 366 differentially expressed miRNA-mRNA interaction pairs were obtained ().

Table 1

Representative interactions between lncRNAs and miRNAs for malignant PCCs

lncRNA	miRNA
IGF2-AS	hsa-mir-193b
AP002478.1	hsa-mir-195, hsa-mir-497, hsa-mir-508
FAM87B	hsa-mir-193b
TCL6	hsa-mir-507, hsa-mir-195, hsa-mir-497, hsa-mir-193b
BX255923.1	hsa-mir-506
FAM66C	hsa-mir-195, hsa-mir-497, hsa-mir-193b, hsa-mir-506
C9orf147	hsa-mir-507
NEXN-AS1	hsa-mir-193b
AC123595.1	hsa-mir-506
AL589765.1	hsa-mir-195, hsa-mir-497
ALDH1L1-AS2	hsa-mir-195, hsa-mir-497, hsa-mir-193b, hsa-mir-508
LINC00485	hsa-mir-195, hsa-mir-497, hsa-mir-193b

PCCs, pheochromocytomas.

Figure 2

The selection of differentially expressed miRNA-mRNA interaction pairs and construction of ceRNA network in malignant PCCs. (A) 2-set venn diagram presenting the intersection between the online predicted target genes and differentially expressed mRNAs. (B) The lncRNA-miRNA-mRNA ceRNA network. Diamonds, rectangles and circles represent lncRNA, miRNA and mRNA, respectively. All shapes in red and blue colors stand for upregulation and downregulation, respectively. PCCs, pheochromocytomas.

Table 2

Representative interactions between MiRNA and target mRNA in malignant PCCs

miRNA	mRNA
hsa-mir-193b	ALDH3A2, RAB27B, ADCYAP1R1, RPS6KA1, ZNF365, SH3BP5, DAAM2, PPARGC1B, NSG1, TGFBR3, REST, YBX3, GNAO1, LEF1, PLXNA2, CHDH, NOVA1, PER3, ATF7IP2, MYCN, PTK7, PLXNC1, STMN1
hsa-mir-195	PAG1, RUNDC3B, TRPC1, PCDHA1, CAMKV, MARCH4, SSBP2, DMTN, TMEM178B, PRKCD, FAM189A1, CLCN5, SEC14L4, INPP5J, GLP2R, FGFR4, UNC80, SIDT1, SEMA3D, CASKIN1, CACNB1, STX1A, ZC2HC1A, UNC13A, SCN3A, CDK5R1, KIF21A, ITGA2, XKR7, SEMA6D, ANK2, SPSB4, TSPAN5, CYB561A3, CLIP4, LRRTM2, ZNF275, REEP1, PLAG1, NSG1, RSPO3, PCDHAC1, SLC12A2, PFKFB4, SH3BGRL2, TUBB2A, RPS6KA6, TGFBR3, SCAMP5, PLSCR4, SYNJ1, PDE4D, THRB, DLGAP3, KIF5A, PCDHAC2, DZIP1, NSMF, PLXNA2, SLC4A8, DIRAS2, TNFSF9, PDCD4, CLCN4, PLD1, ALDH3B1, SH2D2A, RNF165, SIRT4, DNAJB4, ST8SIA3, RET, COBLL1, RAB15, C3orf70, ZNF704, HMGA1, RIMS3, MMD, MAP3K9, AQP11, MRAS, TMEM41A, MAMSTR, STXBP1, RAB9B, MGAT4A, NOTCH2, PTPN14, FAM110C, RNF43, SLC36A1, ABCG4, KL, FGFR1, TMEM138, RACGAP1, SCN2A, SYT4, CDS2, AATK, TMEM100, G0S2, NAV1, TUBA1A, ARHGEF9, RORA, DYNC1I1, MFAP5, SGK1, PHACTR2, SLIT2, POLE4, SGPL1, PRKG1, TPPP3, YWHAH, PITPNM3, PLXNC1, PTPRJ, MYT1L, PGM1, MASP1, ADAMTSL1, SYDE2, SCN8A, CLDN2, RASEF, RCAN3, TLCD2, STEAP3, AMOTL1, ANKRD13B, KATNAL1, RAB9A
hsa-mir-497	CYB561A3, PTPN14, PLXNC1, KIF21A, CLCN5, AQP11, SYNJ1, PAG1, SYT4, MGAT4A, TLCD2, MASP1, PLAG1, RPS6KA6, PGM1, MRAS, ANK2, CDK5R1, UNC80, RCAN3, KL, CLDN2, MMD, SCN3A, C3orf70, SYDE2, RAB15, PLD1, SLIT2, KATNAL1, RNF43, ST8SIA3, AMOTL1, RASEF, PLXNA2, KIF5A, RACGAP1, SLC12A2, TSPAN5, PCDHA1, SH3BGRL2, MFAP5, SH2D2A, ZC2HC1A, RAB9A, SSBP2, COBLL1, SCN8A, TMEM41A, MYT1L, RIMS3, CLCN4, REEP1, RUNDC3B, G0S2, RAB9B, STEAP3, DNAJB4, SPSB4, SLC4A8, DYNC1I1, NAV1, INPP5J, LRRTM2, CASKIN1, PRKG1, TMEM178B, NOTCH2, UNC13A, PTPRJ, MAMSTR, RET, FGFR1, ARHGEF9, MARCH4, SGK1, MAP3K9, SGPL1, SIDT1, SEC14L4, PFKFB4, TPPP3, XKR7, RNF165, RSPO3, YWHAH, FAM189A1, DZIP1, SIRT4, TRPC1, DMTN, ZNF275, DLGAP3, FAM110C, TGFBR3, GLP2R, SEMA3D, PDE4D, PCDHAC1, TUBA1A, PITPNM3, PCDHAC2, NSMF, TUBB2A, CLIP4, SCN2A, PLSCR4, ZNF704, SEMA6D, ADAMTSL1, NSG1, DIRAS2, RORA, AATK, ABCG4, CACNB1, ITGA2, FGFR4, TNFSF9, SLC36A1, TMEM100, STX1A, THRB, STXBP1, ALDH3B1
hsa-mir-506	KIAA1024, NIPA1, KCNK2, SNED1, PAK3, PHACTR2, PLIN3, FSTL3, SLC5A8, SERTAD4, NAP1L3, PCDHAC2, ITM2C, STEAP3, SLC40A1, ARHGEF37, TTL, SLITRK4, GATA6, TTC7A, CCR10, SLC18A2, TMEM41A, PAPSS2, PCDHAC1, SLC12A2, ZFR2, PCCB, IQGAP2, DHCR24, PPP1R13L, CPEB1, GRIA2, SLC31A1, KIAA1210, AMOTL1, GDAP1L1, KCNQ5, PCDH1, SEMA6D, KIF3A, VSTM4, PLXNA3, PCDHA1, SGPL1, ACADVL, RAB38, MYZAP, PGM1, ANGPTL1, ASB4, GRID1, TSPAN6, SDC1, RBM47, RYR2, P4HA2, PTPRJ, GNG2, RPS6KA1, KIF26A, TUB, TLN2, LDLRAP1, YAP1, ALDH1L2, H6PD, PTPN14, OAF, NBPF1, MGAT4A, LRRC58, SGK1, ILDR1, PRKG1, AKAP6, NYAP2, SPOCK3, RAB11FIP4, RAB3D, GPAM, SRPK3, SLC17A5

PCCs, pheochromocytomas.

Construction of ceRNA network in malignant PCCs

To further investigate the molecular mechanism of malignant PCCs, we constructed a ceRNA network on the basis of the interactions of 26 lncRNA-miRNA and 366 miRNA-mRNA pairs. Finally, 12 lncRNAs, 6 miRNAs, and 220 mRNAs were included in the ceRNA network (). Subsequently, we performed functional enrichment analysis of these differentially expressed mRNAs using the BNGO tool. Top six GO terms were identified in which differentially expressed mRNAs were greatly enriched in BP, CC, and MF (). Ion transport was the most significant GO term of BP, whereas plasma membrane was largely enriched in CC. For MF, purine nucleotide binding, protein domain specific binding, and nucleotide binding were the three most relevant terms.

Figure 3

The functions of DEmRNA in the ceRNA network were analyzed with DAVID, and prognostic significances of differentially expressed genes for overall survival were performed using Kaplan-Meier survival curves. (A) GO enrichment significance items of DEmRNA in different functional groups. (B) Distribution of DEmRNA in malignant PCCs for different GO-enriched functions. (C,D) Kaplan-Meier survival curves for lncRNAs (C9orf147) and mRNAs (KL, PCDHAC2, RAB11FIP4, RAB15, SCN8A, TTL) associated with overall survival in malignant PCCs (P<0.05). DEmRNA, differentially expressed mRNA; GO, gene ontology; PCCs, pheochromocytomas.

Survival analysis of vital differentially expressed RNAs in malignant PCCs

Kaplan-Meier curves analysis was used in analyzing the correlation of differentially expressed genes in the ceRNA network to overall survival. We found that seven differentially expressed genes were considered to play important roles in the survival outcome of malignant PCCs. LncRNA C9orf147 expression was significantly associated with poor overall survival (log-rank P=0.0343, ). Similar to the lncRNA, high expressions of six mRNAs (Klotho, PCDHAC2, RAB11FIP4, RAB15, SCN8A, TTL) significantly predict poor survival outcome ().

Cross-analysis of gene expression and DNA methylation profiles

A total of 192 differentially methylated genes were obtained from patients with malignant PCCs and compared with those of the controls (). Four of these differentially expressed mRNAs from ceRNAs overlapped with the methylation data (PLSCR4, GATA6, YWHAH, and RNF43). The Pearson's rank test for the genes showed significantly negative correlation between their expression and methylation levels (). Of the four genes, three (PLSCR4, GATA6, RNF43) were downregulated and hypermethylated and one (YWHAH) was upregulated and hypomethylated in tumor tissues. Moreover, the methylation sites of these four genes negatively correlated with expression level were sought ().

Figure 4

Heatmap of the differential methylation of genes in malignant PCCs between 38 tumor tissues and 3 normal tissues. The color from green to red shows a trend from low expression to high expression. PCCs, pheochromocytomas.

Figure 5

Cross-analysis of gene expression and DNA methylation profiles. (A) Four genes with different levels of expression and methylation between tumor and normal group; (B) four genes which expression levels negatively correlated with methylation levels.

Figure S1

Methylation sites which were negatively correlated with gene expression. (A) PLSCR4; (B) GATA6; (C) RNF43; (D) YWHAH.

Discussion

Predicting and managing the malignancy of PCCs remain as a considerable challenge in clinical practice (19). Thus, further understanding of the clinical and molecular characteristics of patients with malignant PCCs is urgently needed. In the present study, we performed a bioinformatics analysis on the differentially expressed lncRNAs, miRNAs, and mRNAs in malignant PCCs and compared them with nontumor tissues based on the TCGA database. We further analyzed the function of malignant PCCs-specific mRNAs by using the GO biological analysis. The GO analysis revealed that the functions of these genes were mainly enriched in the regulation of ion transport, plasma membrane, and nucleotide and protein domain-specific binding.

On the basis of the ceRNA network in malignant PCCs, we analyzed the potential mechanism of interaction among lncRNAs, miRNAs, and mRNAs. Of the 12 lncRNAs identified in the ceRNA network, high lnc-C9orf147 expression was significantly related to poor overall survival. According to our ceRNA network, lnc-C9orf147 is involved in the interaction with miR-507, which is a tumor suppressor in many cancer cells (20). However, the involvement of lnc-C9orf147 in other diseases or malignancies and its function have not been reported. Thus, we considered lnc-C9orf147 as a promising biomarker that needs further investigation and definition.

In addition, six differentially expressed mRNAs from the ceRNA network were significantly associated with the prognosis of malignant PCCs. One of the six mRNAs is Klotho (KL), which produces two types of proteins: α-Klotho and β-Klotho. These proteins are the essential components of endocrine fibroblast growth factor (FGF) receptor complexes. The FGF-KL endocrine system has a crucial role in the pathophysiology of aging-related disorders, including diabetes, cancer, arteriosclerosis, and chronic kidney disease (21). KL is dysregulated in several cancers, and its downregulation is related to the promoted proliferation and reduced apoptosis of tumor cells (22). Furthermore, the KL protein provides new insights into cancer target treatment. KL overexpression has synergistic effects with cisplatin to inhibit the proliferation of drug-resistant lung cancer cells in a dose-and time-dependent manner (23). Furthermore, KL overexpression inhibits the progression of ovarian cancer in nude mice partly via the inhibition of systemic inflammation (24). However, our result indicated that KL was upregulated in malignant PCCs and was significantly related to poor overall survival, which differed from other mentioned cancers. The upregulation of SCN8A-mediated invasiveness of cervical cancer cells involves MMP-2 activity, and SCN8A channels are considered therapeutic targets against cancer metastasis (25,26). Few studies have investigated the function and prognostic significance of the other four identified mRNAs.

Regarding the methylation profile analysis associated with gene expression in the ceRNA network, four mRNAs, namely, PLSCR4, GATA6, YWHAH, and RNF43, were identified. PLSCR4, GATA6, and RNF43 were downregulated and hypermethylated, and YWHAH was upregulated and hypomethylated in malignant PCCs. The expression of GATA6 was epigenetically silenced through promoter methylation in gastric cancer cell lines, and aberrant JAK/STAT signaling suppresses TFF1/2 partially through the epigenetic silencing of GATA6. Moreover, gastric cancer patients with high GATA6 methylation tend to have short overall survival (27). The promoter hypermethylation of GATA6 has high frequency in glioblastoma and is an independent predictor in glioblastoma patient outcome (28,29). GATA6 is also hypermethylated in both rhabdomyosarcoma cell lines and primary samples (30). RNF43 is a significantly mutated driver gene in gastric cancer. However, no relevant studies have focused on the gene methylation of PLSCR4 or YWHAH in human cancers and malignancies.

Although this study has shown outcomes with significant clinical significance, we still need to consider several limitations. First, the sample number of the control group from TCGA database was relatively small. Thus, more cases need to be enrolled for analysis. Second, the molecular mechanisms and prognostic roles of lnc-C9orf147 and vital methylated genes in malignant PCCs needs to be further studied in vivo and in vitro.

Conclusions

We constructed a ceRNA network with malignant PCC-specific lncRNAs, miRNAs, and mRNAs in this study. We clarified the mechanism of DNA methylation-associated genes in a ceRNA network in malignant PCCs. Our research increases the understanding of the pathogenesis of malignant PCCs and offers potential genes as underlying therapeutic targets or prognostic biomarkers.

The article’s supplementary files as