Comprehensive analysis of competing endogenous RNA network and 3-mRNA signature predicting survival in papillary renal cell cancer.

Long non-coding RNAs (lncRNAs) can act as competing endogenous RNAs (ceRNAs) to exert significant roles in regulating the expression of mRNAs by sequestering and binding miRNAs. To elucidate the functional roles and regulatory mechanism of lncRNAs in papillary renal cell cancer (pRCC), we conducted a comprehensive analysis of ceRNA network and constructed a mRNA signature to predict prognosis of pRCC.We collected mRNAs and lncRNAs expression profiles of 289 pRCC samples and 32 normal renal tissues, and miRNA expression profiles of 292 pRCC samples and 34 normal samples from The Cancer Genome Atlas (TCGA) database. Differential expressions of RNAs were evaluated by the "edgeR" package in R. Functional enrichment analysis of DEmRNA was performed by DAVID 6.8 and KEGG, while PPI network of top 200 DEmRNAs was conducted using the STRING database. The univariate and multivariate Cox regression were conducted to figure out the candidate DEmRNAs with predictive values in prognosis. Receiver operator characteristic (ROC) curve estimation was performed to achieve the area under the curve (AUC) of the ROC curve to judge mRNA-associated prognosic model. A ceRNA network was established relying on the basis of combination of lncRNA-miRNA interactions and miRNA-mRNA interactions.A total of 1928 DEmRNAs, 981 DElncRNAs, and 52 DEmiRNAs were identified at significance level of |log2Fold Change |>2 and adjusted P-value < .01. A 3-mRNA signatures consisting of ERG, RRM2, and EGF was constructed to predict survival in pRCC. Moreover, a pRCC-associated ceRNA network was constructed, with 57 lncRNAs, 11 miRNAs, and 28 mRNAs.Our study illustrated the regulatory mechanism of ceRNA network in papillary renal cancer. The identified mRNA signatures could be used to predict survival of pRCC.

Introduction

Renal cell cancer (RCC) represents one of the most common malignancies in urological system, with 65,340 estimated new cases and 14,970 estimated deaths in 2018 in the United States nationally.[ Although RCC contains various histopathological types according to the 2016 World Health Organization (WHO) classification,[ clear cell renal cell cancer (ccRCC), papillary renal cell cancer (pRCC), and chromophobe renal cell cancer (chRCC) still remain the top 3 histological types. Compared with ccRCC, pRCC is the second most encountered subtype and characterized with a significantly higher rate of organ-confined tumor and a higher 5-year cancer specific survival rate.[

Radical surgeries exert the most effective role in therapeutic approaches for RCC although several options, including active surveillance, thermal ablation and surgery, exist for management of clinically localized RCC.[ Molecular targeted therapy against vascular endothelial growth factor (VEGF) and mTOR pathways has been developed for metastatic patients or locally advanced RCC cases, however therapeutic response is various and most patients eventually progress.[ Due to the limitations of current therapeutic options, further investigation into the underlying molecular mechanism of RCC is expected to provide promising strategy for cancer therapeutics.

Long chain noncoding RNAs (lncRNAs), as one of the ncRNAs, is characterized with the length longer than 200 nucleotides and deficiency in protein-coding capacity, which were ever considered as simply fake transcriptional noise.[ With the deeper investigation of underlying mechanisms of lncRNA involvement in oncogenic and tumor suppressive pathways, the new cancer diagnostic markers and novel therapeutic targets are emerging.[ Intriguingly, a competing endogenous RNA (ceRNA) hypothesis was proposed by Salmena et al in 2011, in which mRNA, lncRNA, pseudogenes, et al performed as natural miRNA sponges to occupy the miRNA response elements (MREs) to affect miRNA functions.[ lncRNAs can act as ceRNAs to exert significant roles in regulating the expression of mRNAs by sequestering and binding miRNAs. For example, lncRNA GAS5 exerted as a ceRNA to regulate mRNA levels of hZIP1 by sponging miR-223 in the progression of ccRCC. The GAS5/miR-223/hZIP1 axis may serve as a therapeutic strategy for patients with ccRCC.[

To elucidate the functional roles and regulatory mechanism of lncRNAs in pRCC, we conducted a comprehensive analysis of ceRNA network and constructed a mRNA signature to predict prognosis of pRCC.

Methods

Data collection

This study was not involved with new participants, so ethics committee approval is not required. Information on pRCC was retrieved from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/) data portal. RNA expression patterns and clinical data was obtained from TCGA in October 2018. A total of 321 samples including 289 pRCC samples and 32 normal samples were included to provide mRNAs and lncRNAs expression data. The raw miRNA expression data was collected from 292 pRCC samples and 34 normal samples.

Identification of differentially expressed RNAs

Empirical Analysis of Digital Gene Expression Data in R (edgeR) was used to identify differences of mRNA, lncRNA, and miRNA expression between pRCC and normal samples with thresholds of |log2Fold Change(FC)|>2 and adjusted P-value < .01.[

Functional enrichment analysis

We conducted functional enrichment analysis using DAVID 6.8 to reveal the functional roles of differentially expressed mRNAs (DEmRNAs).[ The biological functions, cellular components, and molecular functions of the top 200 DEmRNAs were analyzed, and the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/) were searched for pathways at the significance level set (P < .05 and enrichment score >1.5). Protein–protein interaction (PPI) network for top 200 DEmRNAs was constructed with the STRING database(Version: 10.5, https://string-db.org/ ,minimum required interaction score >0.4). Moreover, Cytoscape v3.7.0 software was applied to demonstrate the PPI networks visually.

Construction of mRNA-associated prognostic model

The univariate Cox regression was initially used to identify pRCC-associated mRNAs by assessing the associations between DEmRNAs and overall survival (OS) with P-value < .01. Then, the multivariate Cox regression was applied to remove the candidate DEmRNAs without independent predictive effects on prognosis and establish a prognostic signature. Receiver operator characteristic (ROC) curve estimation was performed to achieve the area under the curve (AUC) of the ROC curve to judge mRNA-associated prognostic model.

CeRNA network construction

The identified pRCC-associated DEmRNAs, differentially expressed lncRNAs (DElncRNAs) and differentially expressed miRNAs (DEmiRNAs) were included to construct a ceRNA network. Firstly, the interactions between DElncRNAs and DEmiRNAs were evaluated based on the miRcode database (http://www.mircode.org/). Then DElncRNA-associated DEmiRNAs were used to predict targeted mRNAs by the combination of miRDB, miRTarBase, and TargetScan. The interactions between DElncRNA-associated DEmiRNAs and DEmRNAs were evaluated to construct the ceRNA network. A ceRNA network was established relying on the co-expression network of DElncRNAs-DEmiRNAs-DEmRNAs interactions. Cytoscape v3.7 software was used to demonstrate this network visually.

Survival analysis of key members in the ceRNA network

The prognostic values of DEmRNAs, DEmiRNAs, and DElncRNAs in the ceRNA network were evaluated by “survival” package in R. The survival analysis was assessed by Kaplan–Meier method with significance level P-value < .05.

Results

DEmRNAs, DElncRNAs, and DEmiRNAs in pRCC

We identified 1928 DEmRNAs between tumor tissues and normal tissues from TCGA, of which 877 mRNAs were downregulated and 1051 mRNAs were upregulated, as demonstrated in the volcano plot (Fig. 1A). A total of 981 DElncRNAs were obtained with 552 up-regulated lncRNAs and 429 down-regulated lncRNAs (Fig. 1B). With regard to the miRNA expression data, 92 significantly differentially expressed miRNAs (DEmiRNAs) were identified between tumor and normal tissues, among which, 52 miRNAs were upregulated while 40 miRNAs were downregulated (Fig. 1C). The heat maps of DEmRNAs (Fig. 1D), DElncRNAs (Fig. 1E), and DEmiRNAs (Fig. 1F) in pRCC were also demonstrated in Figure 1.

Figure 1

Volcano plots of DEmRNAs (A), DE lncRNAs(B) and DEmiRNAs (C), and red points represent up-regulated RNAs while green points represent down-regulated RNAs with with thresholds of |log2FC|>2 and adjusted P-value < .05. Heatmaps of DE mRNAs (D), DElncRNAs (E), DEmiRNAs (F).

Functional enrichment analysis of DEmRNAs

To understand the functional implications of DEmRNAs, the top 200 DEmRNA were selected to evaluate biological functions, cellular components, and molecular functions by DAVID 6.8. The top 10 DEmRNAs-associated biological functions included chloride transmembrane transport, excretion, transport, sodium ion homeostasis, ion transport, negative regulation of cysteine-type endopeptidase activity, sodium ion transmembrane transport, regulation of membrane potential, ion transmembrane transport, and renal water homeostasis, as demonstrated in Figure 2A. Figure 2B showed the top 10 cellular components including extracellular region, extracellular space, apical plasma membrane, integral component of plasma membrane, basolateral plasma membrane, plasma membrane, extracellular exosome, platelet alpha granule lumen, blood microparticle, and neuronal cell body. The top 10 molecular functions were shown in Figure 2C, including cysteine type endopeptidase inhibitor activity, alkane 1 monooxygenase activity, transporter activity, anion: anion antiporter activity, chloride channel activity, oxidoreductase activity, protease binding, extracellular ligand-gated ion channel activity, ligand-gated ion channel activity, and arachidonic acid epoxygenase activity. The potential pathways involved was discovered by KEGG analysis. These DEmRNAs were enriched in some carcinogenesis related pathways including Wnt signaling pathway and arachidonic acid metabolism (Fig. 2D). The top 200 DEmRNAs were also used to construct a protein–protein interaction (PPI) network with the String database and visualized by Cytoscape v3.7.0 (Fig. 3A).

Figure 2

Bubble charts of Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. (A) The top 10 DEmRNAs-associated biological functions. (B) The top 10 cellular components. (C) The top 10 molecular functions. (D) The potential pathways by KEGG analysis.

Figure 3

(A) Protein–protein interaction (PPI) network of the top 200 DEmRNAs in pRCC. (B) Competing endogenous RNA (ceRNA) network in pRCC. The red indicates upregulated RNAs, and the blue indicates the downregulated RNAs in pRCC. The diamond represents lncRNAs, round rectangle represents miRNAs, and ellipse represents mRNAs.

Identification of mRNA-associated prognostic signature

Univariate Cox regression analysis was conducted to identify prognostic DEmRNAs with regard to overall survival of pRCC patients from TCGA cohort. Then multivariate Cox regression was performed, and 3 mRNAs including ERG, RRM2, EGF were finally identified to establish a prognostic signature. The survival risk score was calculated as follows: survival risk score = (0.34745) × ERG + (0.76291) × RRM2 + (0.33446) × EGF. The prognostic power of the 3-mRNA signature was assessed based on the area under the receiver operating characteristic curve value. The sensitivity and specificity of the 3-mRNA signature in assessing prognosis were both good with AUC 0.815 while using a median risk score of 0.835 as the cutoff (Fig. 4A). The pRCC patients were finally divided into high risk group (risk score ≥0.825, n = 143) and low risk group (risk score <0.825, n = 144). The 2 groups demonstrated significantly different prognosis in the survival curve (Fig. 4B). The 5-year overall survival rates were 65.5% (95%CI: 54.9%–78.1%) in the high-risk group and 84.5% (95%CI: 74.8%–95.4%) in the low-risk group respectively.

Figure 4

Prognostic evaluation of the 3-mRNA signature in pRCC. (A) ROC curve analysis of the 3-mRNA signature. (B) Kaplan–Meier analysis of overall survival in high-risk group and low-risk group assigned by prognostic index with the 3-mRNA signature.

Construction of lncRNA-miRNA-mRNA ceRNA network

To explore the underlying interactions among DEmRNAs, DElncRNAs, and DEmiRNAs in pRCC, a ceRNA network was constructed using Cytoscape. Firstly, the DElncRNA-DEmiRNA interactions were predicted among the 981 DElncRNAs and 92 DEmiRNAs. Consequently, 59 DElncRNAs were identified to interact with 15 DEmiRNAs based on miRcode database (Table 1). Then, the identified 15 DElncRNA-associated DEmiRNAs were used to identify targeted mRNAs through miRDB, miRTarBase, and TargetScan database. After exclusion of DEmiRNA-associated mRNA not containing in DEmRNAs, only 11 DEmiRNAs were predicted to exert interaction with 28 DEmRNAs (Table 2). Moreover, a lncRNA-miRNA-mRNA network consisting of 11 DEmiRNAs, 28 DEmRNAs, and 57 DElncRNAs was constructed based on the combination of the lncRNA-miRNA and miRNA-mRNA interactions as demonstrated in Fig 3B.

Table 1

Interactions between DEmiRNAs and DElncRNAs.

Table 2

Correlation between DEmiRNAs and DEmRNAs.

Survival analysis with key members in ceRNA network

To evaluate the associations between the key members of ceRNA network and prognosis of pRCC patients, Kaplan–Meier method was used to evaluate the relationship between expression profiles of key members and overall survival. 12 of 57 DElncRNAs were associated with the prognosis of pRCC, of which 6 lncRNAs including AP000525.1, DNM3OS, GDNF-AS1, GLIS3-AS1, LINC00310, and LINC00462 were positively connected with overall survival of pRCC (Fig. 5), while the other 6 lncRNAs including COL18A1-AS1, CRNDE, GAS6-AS1, GPC5-AS1, LINC00327, and SACS-AS1 were negatively correlated with OS (Fig. 6). 5 of 11 DEmiRNAs were associated with OS, including hsa-mir-145, hsa-mir-211, hsa-mir-214, hsa-mir-216a, hsa-mir-217 (Fig. 7).14 of 28 DEmRNAs were connected with OS, containing E2F1, E2F2, ERG, GPC5, HAS2, HOXC13, IL11, LDLR, POLQ, RAPGEF4, RRM2, SELE, SFRP1 and SLC7A11 (Fig. 8). Notably, ERG and RRM2 were also included in the 3-mRNA prognostic signature.

Figure 5

Kaplan–Meier survival curves for 6 lncRNAs positively associated with overall survival of pRCC.

Figure 6

Kaplan–Meier survival curves for 6 lncRNAs negatively associated with overall survival of pRCC.

Figure 7

Kaplan–Meier survival curves for 5 miRNAs associated with overall survival of pRCC.

Figure 8

Kaplan–Meier survival curves for 14 mRNAs positively associated with overall survival of pRCC.

Discussion

Papillary renal cell cancer, which accounts for 15% to 20% of RCC, is highly heterogeneous, characterized with various morphologic features and biologic behaviors.[ Carcinogenesis and progression of pRCC are complex processes, involving numerous genetic and epigenetic alterations, and create a diversity of biological effects, which contribute to the molecular basis of metastasis and drug resistance.[

Genomic technology development revealed that protein-coding RNAs only take the small percentage 2% of human genome while the majority was transcribed into non-coding RNAs.[ Increasing evidence indicates the important role of lncRNAs in the regulation of proliferation, apoptosis, metastasis, and metabolism.[ The lncRNAs-miRNAs-mRNAs regulation network is associated with cancer occurrence and progression. Yin et al conducted an integrated analysis of ceRNA in ccRCC and constructed a 9-lncRNA signature to predict prognosis.[ With respect to the limited knowledge of ceRNA on pRCC from the perspective of TCGA data, this study demonstrated the interactions and potential crosstalk between RNAs and built pRCC-associated ceRNA network as well as a mRNA signature to predict survival in pRCC.

In this study, 1928 DEmRNAs were identified between tumor tissues and normal tissues from TCGA. The top 200 DEmRNAs were chosen to conduct GO and pathway analysis, in which Wnt signaling pathway and arachidonic acid metabolism were associated with carcinogenesis of pRCC. Univariate and multiunivariate Cox regression analysis were conducted to figure out the correlations between DEmRNA expression profile and overall survival of patients. A 3-mRNA signature was established, including ERG, RRM2, EGF. All coefficients of the 3 mRNAs were positive, indicating that higher expression of them associated with worse prognosis of pRCC.

ERG is a widely investigated oncogene, involving in hematopoiesis, chondrocyte maturation, and bone development and in apoptosis and cell migration.[ Due to the gene fusion with the promoter region of the androgen-induced TMPRRSS2 gene, ERG has become highly correlated with prostate cancer development.[ RRM2, as the small subunit of ribonucleotide reductase, is an essential enzyme involved in DNA replication and repair.[ Overexpression of RRM2 has been showed to be associated with aggressiveness and prognosis of bladder cancer,[ head and neck cancer,[ adrenocortical cancer,[ breast cancer,[ and pancreas adenocarcinoma.[ EGF encodes a member of the epidermal growth factor superfamily and exerts an important role in the growth, proliferation and differentiation of numerous cell types.[ Dysregulation of EGF has been recognized as an essential molecular event in carcinogenesis.[ All of the 3 indicators were widely investigated and proved to be involved in the carcinogenesis and cancer progression, and the signature consisting of the 3 mRNAs help to predict prognosis in pRCC.

With respect to the key members of lncRNAs associated with survival of pRCC, we retrieve literature investigating the molecular mechanism of previous lncRNAs.DNM3OS was found to confer significant radio-resistance by regulating DNA damage response in esophageal squamous cell carcinoma.[ DNM3OS overexpression in ovarian cancer was showed to significantly correlate to worse overall patient survival, due to its contribution to epithelial-to-mesenchymal transition (EMT) associated with metastasis.[ Circulating long non-coding RNAs (lncRNAs) may act as diagnostic markers of cancer and predictors of their histological classification and survival. GLIS3-AS1 was proposed as one of 8-lncRNA signature, which helps in distinguishing malignant Intraductal Papillary Mucinous Neoplasms of the Pancreas (IPMNs) from benign ones with greater accuracy than standard clinical and radiologic features.[ LINC00310 was elucidated to promote cell proliferation by regulating c-Myc expression and serum LINC00310 expression exerted diagnostic significance in breast cancer.[ LINC00462 promotes hepatocellular carcinoma progression through the PI3K/AKT signaling pathway.[ In addition, LINC00462 also act as a target of miR-665 and activate the SMAD2/3 pathway by enhancing the expression levels of TGFBR1 and TGFBR2, and promote tumorigenesis in pancreatic cancer.[ COL18A1-AS1was found to exert as a ceRNA in Cholangiocarcinoma[ and in clear cell renal cell cancer.[ Analysis based on 3 human exon arrays available from the public Gene Expression Omnibus, suggested that expression of lncRNA CRNDE was increased in both ccRCC and metastatic ccRCC samples, which indicated the role of it in cancer progression.[ CRNDE has been investigated in various kinds of cancer and serve as a new predictive factor for prognosis and diagnosis in different types of cancers.[ Low expression of lncRNA GAS6-AS1 could act as a poor prognostic indicator in patients with non-small cell lung cancer.[ To the best of our knowledge, investigations in AP000525.1, GDNF-AS1, GPC5-AS1, LINC00327 are still lacking. However, there are still a deficiency of reports on correlations between AP000525.1, GDNF-AS1, GPC5-AS1, LINC00327, and cancer. Not only do the lncRNAs act as ceRNA to regulate mRNA expressions, however all of them were associated with overall survival in pRCC, so they have more values to investigate in.

There are still some advantages and limitations in this study. The TCGA provide enough sample size and comprehensive measurement data, providing a relatively reliable basis for bioinformatics analysis. The application of various bioinformatic analysis methods can provide mRNA-dependent predictive model and explore the potential ceRNA network, giving a direction for further mechanism research. However, several lncRNAs with clinical significance in pRCC are still required to be validated for the underlying molecular mechanism.

Conclusion

In conclusion, our study illustrated the regulatory mechanism of ceRNA network in pRCC and constructed a 3-mRNA predictive signature. The identified mRNA signature could be used to predict survival of pRCC.

Author contributions

Conceptualization: Xin Zhu, Mi Zhou.

Formal analysis: Xin Zhu, Mi Zhou.

Methodology: Xin Zhu, Jianyu Tan.

Software: Xin Zhu.

Writing - original draft: Xin Zhu, Mi Zhou.

Data curation: Jianyu Tan, Zongjian Liang.

Investigation: Mi Zhou.

Project administration: Mi Zhou.

Writing – review & editing: Mi Zhou.