Risk Scores Based on Six Survival-Related RNAs in a Competing Endogenous Network Composed of Differentially Expressed RNAs Between Clear Cell Renal Cell Carcinoma Patients Carrying Wild-Type or Mutant Von Hippel–Lindau Serve Well to Predict Malignancy and Prognosis

Clear cell renal cell carcinoma (ccRCC) carrying wild-type Von Hippel–Lindau (VHL) tumor suppressor are more invasive and of high morbidity. Concurrently, competing endogenous RNA (ceRNA) network has been suggested to play an important role in ccRCC malignancy. In order to understand why the patients carrying wild-type VHL gene have high degrees of invasion and morbidity, we applied bioinformatics approaches to identify 861 differentially expressed RNAs (DE-RNAs) between patients carrying wild-type and patients carrying mutant VHL from The Cancer Genome Atlas (TCGA) database, established a ceRNA network including 122 RNAs, and elected six survival-related DE-RNAs including Linc00942, Linc00858, RP13_392I16.1, hsa-miR-182-5p, hsa-miR-183-5p, and PAX3. Examining clinical samples from our hospital revealed that patients carrying wild-type VHL had significantly higher levels of all six RNAs than those carrying mutant VHL. Patients carrying wild-type VHL had significantly higher risk scores, which were calculated based on expression levels of all six RNAs, than those carrying mutant VHL. Patients with higher risk scores had significantly shorter survival times than those with lower risk scores. Therefore, the risk scores serve well to predict malignancy and prognosis.

Introduction

Kidney cancer is one of the top 10 most common tumors, the second-ranked malignant tumor in urologic system, and contributes 2% to 3% of the human malignant tumors worldwide (1, 2). There are 400,000 new cases of kidney malignant tumors of which 175,000 lead to death per year (3). Among patients with renal cancer at early stage and even some of renal cancer patients who do not have any specific symptom regardless of the stage of the disease, only about 30% of patients diagnosed with renal cancer was found to be at initial stage of malignancy, but the other 70% were found to be already at late stage of malignancy (4, 5). About 90% of the cases with kidney cancers are those with renal cell carcinomas, of which 70% are with clear cell renal cell carcinomas (ccRCCs), which are either hereditary or sporadic (6). Among many genetic factors that were found to be related to ccRCC, the inactivation of Von Hippel–Lindau (VHL) tumor suppressor gene occurs in 80% of patients and is the most common event (7–9), leading to the interruption of both vascular endothelial growth factor (VEGF) and mammalian target of rapamycin (mTOR) signaling pathways, which serve as the targets of most of the drugs currently used in clinical practice (10, 11). Although the inactivation of VHL gene alone is not sufficient to cause tumors (12, 13), a study of ccRCC patients from Germany suggested that low expression of VHL was identified as a risk factor for worse overall survival (OS) (14). However, it was reported that there is no difference in survival between patients carrying wild-type or mutant VHL or between patients with different levels of either VHL protein or messenger RNA (mRNA), and patients with nonsense mutations in exon 1 have a worse prognosis based on a study of ccRCC patients from Brazil (15), suggesting that loss of VHL function is directly related to the prognosis of ccRCC patients. A recent study indicated that about 40%–60% of patients with sporadic ccRCC carry wild-type VHL, and such tumors are more invasive and lead to a dramatic reduction of survival rates as compared with those carrying mutant VHL (16, 17). Therefore, it is more urgent to identify other factors that are closely associated with patients’ survival in addition to VHL mutation.

Long non-coding RNA (lncRNA) is a type of transcript with more than 200 nucleotides that does not encode a protein (18). It was initially considered to be “noise” generated by genome transcription without biological functions. More and more evidences show that lncRNA is involved not only in various physiological processes, such as cell proliferation, differentiation, and apoptosis, but also in pathological processes of various diseases (18). The long intergenic non-coding RNA (lincRNA) is the most common form of lncRNA. LncRNA is reported to be related to a variety of tumor diseases such as liver cancer (19), gastric cancer (20), oral squamous cell carcinoma (21), bladder cancer (22), prostate cancer (23), and kidney cancer (24). MicroRNAs (miRNAs) are small single-stranded non-coding molecules of about 22 nucleotides long and target the 3′UTR or CDS of mRNAs to form RNA-Induced Silencing Complexes (RISCs) to induce the degradation of the correspondent mRNA (25). Circular RNA (circRNA) is homologous to its target gene sequence and may act as a molecular sponge. Since lncRNA and mRNA share sequence homology, it is also possible to combine miRNA with lncRNA or circRNA to change expression levels of their downstream target genes. Therefore, miRNAs, lncRNAs, mRNAs, and other non-coding RNAs form a large-scale network of competing endogenous RNAs (ceRNAs), which are currently considered to be involved in the regulation of a variety of tumor diseases and play an important role in different stages of tumorigenesis (26).

Here, we constructed a ceRNA network after the sets of differentially expressed lncRNAs (DE-lncRNAs), miRNAs (DE-miRNAs), and mRNAs (DE-mRNAs) between patients carrying wild-type and those carrying mutant VHL, which were common in multiple database, were identified. Further identification of six survival-related RNAs including Linc00942, Linc00858, RP13_392I16.1, hsa-miR-182-5p, hsa-miR-183-5p, and PAX3, in the ceRNA network led to the building of a model in which risk scores were successfully applied to predict malignancy and prognosis of ccRCC.

Materials and Methods

Selection of Patient Data in The Cancer Genome Atlas Database

The selection of patient data in The Cancer Genome Atlas (TCGA) database was performed following the workflow as shown (Figure 1A). The expression data (read count) of ccRCC patients in TCGA were downloaded and integrated using R’s GDCRNATools package. The sample size for mRNA and lncRNA was 530, and miRNA was 516. The information about mutation of those samples was from cBioPortal (http://www.cbioportal.org/). The information about OS of those samples was from xenabrowser (https://xenabrowser.net/). Since RNA sequencing data were obtained directly from TCGA, approval by any ethics committee was not required.

FIGURE 1

Figure 1 A large number of DE-RNAs are identified between patients carrying either wild-type or mutant VHL. (A) A diagram showing the workflow of the study. (B–D) Volcanic maps showing the distribution of DE-lncRNAs (B), DE-miRNAs (C), or DE-mRNAs (D). Red, blue, and black dots represent those RNAs that are significantly upregulated, significantly downregulated, and not significantly changed, respectively. VHL, Von Hippel–Lindau; DE-RNAs, differentially expressed RNAs; DE-lncRNAs, differentially expressed long non-coding RNAs; DE-miRNAs, differentially expressed microRNAs; DE-mRNAs, differentially expressed messenger RNAs.

Identification of Differentially Expressed Long Non-Coding RNAs, Differentially Expressed MicroRNAs, and Differentially Expressed Messenger RNAs Between Clear Cell Renal Cell Carcinoma Patients Carrying Wild-Type or Mutant Von Hippel–Lindau

Among the 530 ccRCC samples with information of mRNA and lncRNA expression, samples of 446 patients were genotyped as either wild type (223 samples) or mutant (223 samples). Among the 516 samples with information of miRNA, samples of 220 patients were genotyped as wild-type, and 221 patients were genotyped as mutant VHL. The RNA read count of samples carrying wild-type and mutant VHL was standardized and analyzed with the DESeq2 package. DE-lncRNAs, DE-miRNAs, and DE-mRNAs were elected with |log2foldchange| > 1 and Adj.p.val < 0.05 as a threshold according to the method of Benjamini–Hochberg.

Construction of Competing Endogenous RNA Regulation Network

The interaction between lncRNA and miRNA came from miRcode database (http://www.mircode.org/). The interaction between miRNA and mRNA came from miRDB (http://mirdb.org/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), and TargetScan (http://www.targetscan.org/vert_72/) databases. In order to improve the reliability of the results, we selected the intersection of DE-miRNA–DE-mRNA interaction pairs found in the three databases as candidate genes to construct the ceRNA network. Finally, we used Cytoscape software to visualize the ceRNA network.

Enrichment Analyses of the Competing Endogenous RNAs

The Gene Ontology (GO) molecular function enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of DE-mRNA in the ceRNA network were analyzed with R’s ClusterProfiler package. p < 0.05 was considered statistically significant.

Survival Analysis of Differential RNAs in the Competing Endogenous RNA Regulation Network

The expression level of each DE-lncRNA, DE-miRNA, and DE-mRNA in the ceRNA network in each patient was classified into either high- or low-expression group with their median as a threshold. The survival curve for each group was plotted through the Kaplan–Meier (K–M) survival analysis. p < 0.05 was considered statistically significant.

Univariate and Multivariate Cox Regression Analyses

Univariate and multivariate Cox regression analyses were used to screen the survival-related differentially expressed RNAs (DE-RNAs) between patients carrying wild-type and those carrying mutant VHL to eliminate confounding factors, reduce the number of DE-RNAs, and calculate a hazard ratio (HR) and a 95% confidence interval (95% CI) for each variable. Genotypes with an HR greater than 1 or less than 1 and a Wald test p-value less than 0.05 were considered to be the ones that significantly affected patient survival.

External Verification of Survival-Related RNAs With Clinical Tissue Samples

Tissue samples were collected from 21 postoperative ccRCC patients who underwent surgery but did not receive any anti-tumor radiotherapy or chemotherapy before surgery at The Fifth Affiliated Hospital of Guangzhou Medical University from January 2018 to May 2020. All samples were frozen and stored in liquid nitrogen immediately after sampling and subjected to DNA sequencing to determine if their VHL gene was wild type or mutant. The expression levels of VHL protein in the samples were determined by immunohistochemistry analyses with a polyclonal antibody against VHL protein (Cat# 16538-1-AP, from ProteinTech^®). The expression levels of target RNAs in the samples were verified by qRT-RCR. Tumor samples were staged according to the 2010 American Joint Committee on Cancer (AJCC) TNM stage and clinical stage system, and their histopathological grades were determined with the Furman grade method. The study was approved by the Ethics Committee of The Fifth Affiliated Hospital of Guangzhou Medical University. Patients and their family members had been fully informed, and signed consent forms indicated that their samples would be used for scientific research.

Establishment of a Prognostic Risk Score System and Clinical Correlation Analysis for Differentially Expressed RNAs in the Competing Endogenous RNA Regulation Network

A risk score system was established using the DE-RNAs selected by multivariate Cox regression analysis. To predict the patient’s prognosis, risk scores were calculated with the formula “Risk score = ∑(βi * Xi)” in which “i” was the number of characteristic genes, “β” the correlation coefficient of mRNA in the multivariate Cox regression analysis, and X the levels of gene expression after log2 conversion. Pearson’s chi-squared test and Fisher’s exact test were applied to analyze the correlation between patient’s risk score and clinicopathological characteristics in TCGA data set. p < 0.05 was considered to be statistically significant. The t-test was used to compare the difference of risk scores between clinical features that were significantly related to OS. p < 0.05 was considered to be statistically significant.

Gene Set Enrichment Analysis (GSEA) was applied to analyze the signaling pathways associated with patients carrying wild-type VHL and high risk score. The gene expression data from patients carrying high risk score and low risk score in the c2.cp.kegg.v7.0.symbols.gmt database were subjected to estimate the normalized enrichment scores (NESs), p-values, and q-values after false discovery rates (FDRs) were adjusted for each signal pathway with the GSEA software.

Results

Multiple Differentially Expressed Long Non-Coding RNAs, Differentially Expressed MicroRNAs, and Differentially Expressed Messenger RNAs Are Identified Between Patients Carrying Either Wild-Type or Mutant Von Hippel–Lindau

We analyzed the levels of lncRNAs, miRNAs, and mRNAs between ccRCC patients carrying wild-type and those carrying mutant VHL in TCGA database; displayed the distribution of those DE-RNAs in form of volcano map; and defined those RNAs with |log2foldchange| > 1 and Adj.p.val < 0.05 as DE-RNAs. A total of 642 DE-mRNAs (635 upregulated and seven downregulated), 210 DE-lncRNAs (205 upregulated and five downregulated), and 19 DE-miRNAs (18 upregulated and one downregulated) were identified (Figures 1B–D).

A Competing Endogenous RNA Network Is Established

We searched the miRcode database, identified 210 miRNAs that were targeted by DE-lncRNAs, and found interactions between 10 DE-miRNAs and 98 DE-lncRNAs. We searched the miRDB, miRTarBase, and TargetScan databases to identify interactions between miRNA–mRNA interaction and found 10 DE-miRNA and 14 DE-mRNA interactions (Figure 2A). Eventually, a ceRNA regulation network including 290 edges and 122 nodes (98 lncRNAs, 10 miRNAs, and 14 mRNAs) was established between patients carrying wild-type or mutant VHL as shown in Figure 2B.

FIGURE 2

Figure 2 A ceRNA network is created. (A) A Venn diagram showing the miRNA–mRNA interaction as predicted based on three databases. (B) A diagram showing the established ceRNA network. Orange, red, and blue circles represent DE-lncRNAs, DE-miRNAs, and DE-mRNAs, respectively. Gray lines represent interactions among lncRNAs, miRNAs, and mRNAs. ceRNA, competing endogenous RNA; DE-lncRNAs, differentially expressed long non-coding RNAs; DE-miRNAs, differentially expressed microRNAs; DE-mRNAs, differentially expressed messenger RNAs.

The Signaling Pathways Affected by Genes in Competing Endogenous RNA Network Are Predicted

We performed GO enrichment and pathway enrichment analysis on 14 differential genes in the ceRNA network and displayed the top 10 candidate signaling pathways based on the rich factors (the number of differential genes in the GO term divided by the total number of genes in the GO term). The biological processes were mainly related to some kinase activity, growth factor receptor binding, receptor regulatory activity, and so on; cell components were mainly involved in the formation of cell adhesion, synapses and axon membrane structures, and mast cell granules; and molecular functions were mostly related to the kidney functions such as renal vesicle development, epithelial cell differentiation during kidney development, kidney capsule morphogenesis, and renal vesicle morphogenesis (Figure 3A). KEGG enrichment search identified not only some kidney-related pathways such as collecting duct acid secretion and vasopressin-regulated water absorption but also some cancer-related pathways such as PI3K-AKT signaling pathway and MAPK signal pathway as shown in the bubble chart (Figure 3B).

FIGURE 3

Figure 3 The signaling pathways affected by genes in the ceRNA network are identified. (A) A plot showing the signaling pathways affected by genes in the ceRNA network as detected via GO enrichment analyses. (B) A bubble chart showing the signaling pathways affected by genes in the ceRNA network as detected via KEGG pathways analysis. ceRNA, competing endogenous RNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. *p < 0.05, **p < 0.01.

Survival-Related RNAs in Competing Endogenous RNA Network Are Identified Through Univariate and Multivariate Cox Regression Analyses

We estimated the significance of 98 lncRNAs, 10 miRNAs, and 14 mRNAs in the ceRNA network with OS of patients and identified 18 lncRNAs, four miRNAs, and three mRNAs, which significantly impacted OS either positively or negatively (Table 1). Univariate regression analysis found that these RNAs were significantly related to the aggressiveness of disease. Multivariate Cox regression analysis found that LINC00942, LINC00858, RP13_392I16.1, hsa-miR-183-5p, hsa-miR-182-5p, and PAX3 were considered as independent factors that impacted the prognosis of patients (Table 1). Therefore, these six RNAs were selected to build the risk score system.

TABLE 1

Table 1 Univariate and multivariate cox regression analysis of survival-related lncRNAs, miRNAs, and mRNAs.

The Survival-Related RNAs Detected in Databases Are Verified With External Clinical Samples

A total of 21 ccRCC tissue samples were collected in our hospital and subjected to gene sequencing and immunohistochemistry analysis. Using three sets of primers covering the three exons of VHL gene (Figures 4A, B), we identified 18 samples carrying wild-type VHL and three samples carrying mutant VHL as shown in Figure 4C. The identification of high percentage of patients carrying wild-type VHL gene was a surprise to us and triggered us to investigate if the wild-type gene has been suppressed due to other epigenetic modifications. Patients carrying wild-type VHL gene expressed higher levels of VHL protein than those carrying mutant VHL (Figures 4D–G), suggesting that the expression of wild-type gene was not suppressed in general by other epigenetic events such as DNA methylation. The qRT-PCR analyses using primer sets listed (Figure 5A) demonstrated that levels of LINC00942, LINC00858, RP13_392I16.1, hsa-miR-183-5p, hsa-miR-182-5p, and PAX3 were significantly higher in patients carrying wild-type VHL than in patients carrying mutant VHL (Figures 5B–G). Therefore, the results obtained from databases were verified with external clinical samples.

FIGURE 4

Figure 4 Determination of VHL gene expression in the 21 patients. (A) A diagram showing the structure of VHL gene and positions of primers used for gene sequencing. (B) The list of primers used to sequence VHL gene. (C) The description of three detected VHL mutants. (D–F) Representative images showing the low expression levels of VHL protein in patients carrying mutant VHL (D) or wild-type VHL (E) and high expression levels of VHL protein in patients carrying wild-type VHL as detected by immunohistochemical staining in ccRCC patients enrolled in our hospital. (G) A plot showing the percentage of patients with low expression levels of VHL protein in patients carrying wild-type or mutant VHL gene. VHL, Von Hippel–Lindau; ccRCC, clear cell renal cell carcinoma. *p < 0.05.

FIGURE 5

Figure 5 The survival-related RNAs detected in databases are verified with external clinical samples. (A–G) A list of primers used to determine (A) and scatter diagrams showing (B–G) the expression levels of survival-related RNAs LINC00942 (B), LINC00858 (C), RP13.392I16.1 (D), hsa-miR-183-5p, hsa-miR-182-5p (F), and PAX3 (G) as detected by qRT-PCR in 21 clinical patients (wild-type n = 18 and mutant n = 3). *p < 0.05, *p < 0.01.

Risk Scores Are Significantly Different Between Patients with Different Clinical Characteristics

Based on a ratio of 3 to 1, we divided 217 ccRCC patients carrying wild-type VHL into a training set (163 cases) and a test set (54 cases), attributed each patient in the training set with a risk score, and separated patients into a high-risk group (81 patients) and a low-risk group (82 patients) with the median of risk scores 1.77 as the threshold. Patients in the high-risk group exhibited significantly worse prognosis than those in the low-risk group in both training set and test set as well as in all patients (Figure 6A–C). Receiver operating characteristic (ROC) curve analysis detected an area under the curve (AUC) value of 0.79, 0.77, or 0.78 and further confirmed the results from survival analyses (Figures 6D–F). Patients carrying wild-type VHL had significantly higher risk scores than those carrying mutant VHL (Figure 6G). There was no significant difference in risk scores between different age groups, but patients at higher and advanced ccRCC classification stages exhibited significantly higher risk scores (Table 2).

FIGURE 6

Figure 6 Risk scores predict patients’ survival. (A–C) Kaplan–Meier curves showing the relationship between risk scores and patients’ survival in the training set (A), test set (B), and all patients (C). (D–F) ROC curves showing AUC values of the training set (D), test set (E), and all patients (F). (G) A scatter diagram showing differences of risk scores between patients carrying wild type and those carrying mutant VHL. ROC, receiver operating characteristic; AUC, area under the curve; VHL, Von Hippel–Lindau.

TABLE 2

Table 2 The analyses of correlation between risk scores and clinicopathological characteristics of patients in TCGA data set.

Analysis of Signaling Pathway in Clear Cell Renal Cell Carcinoma Patients Carrying Wild-Type Von Hippel–Lindau and High Risk Scores

We conducted GSEA and identified 24 signaling pathways such as the “OLFACTORY_TRANSDUCTION” and “INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION,” which were enriched in the group of patients with high risk scores (Table 3).

TABLE 3

Table 3 GSEA analysis of signaling pathways associated with patients with high risk scores.

Discussion

A ceRNA network has been proposed to play important roles in ccRCC tumorigenesis and aggressiveness. Previous studies have revealed that the LINC01094 enhances the expression of miR-184 and inhibits the expression of SLC2A3, which suppresses the development of ccRCC (27). LINC01426 interacts with insulin-like growth factor 2 mRNA binding protein 1 (IGF2BP1) to enhance expression of CTBP1 in cytoplasm and promote the binding of CTBP1 to the promoter of miR-423-5p in the nucleus to recruit HDAC2 to synergistically inhibit the expression levels of miR-423-5p, leading to an elevation of expression levels of FOXM1 and a promotion of proliferation and migration of ccRCC cells (28). The ceRNA network constructed in this study yielded 98 DE-lncRNAs, 10 DE-miRNAs, and 14 DE-mRNAs. Analyses of the 14 DE-mRNAs revealed their involvement in functions related to kidney. Further analysis of the relationship of DE-RNAs in the ceRNA network with patients’ survival rates led to the identification of LINC00942, LINC00858, RP13_392I16.1, hsa-miR-182-5p, hsa-miR-183-5p, and PAX3, which were used to calculate risk scores for each individual patient as independent factors to successfully predict the malignancy and prognosis of patients.

Three lncRNAs was identified to be survival-related ceRNAs. Currently, there is no report related to any role of RP13_392I16.1 yet, so this lncRNA needs to be further studied. LINC00942 was reported to promote METT14-mediated M6A methylation and regulate the expression and stability of its target gene CXCR4 and CYP1B1 during the initiation and progression of breast cancer (29). In another report, LINC00942 was identified as one of seven lncRNAs in an lncRNA signature related to immune cell infiltration and immune checkpoint blockade of immunotherapy-related molecules and may serve as a prognostic biomarker of hepatocellular carcinoma (30). Similarly, LINC00942 was identified as one of four immune-related lncRNAs to predict prognosis and immunotherapy efficiency in bladder cancer patients (31). It also reported that LINC00942 was identified as one of the two components in a ceRNA network associated with gene number copy variation, which can be used to predict tumor response to drug treatment in patients with lung adenocarcinomas (32). Initially, LINC00858 was identified as a ceRNA impacting miR-422a to control the expression of kallikrein-related peptidase 4, and its high expression was found to be closely correlated to tumor progression of non-small cell lung carcinomas (33). It was reported to be one of the four lncRNAs delivered from database search that were further validated with clinical samples (34) and able to promote colorectal cancer by sponging miR-4766-5p (35) or miR-22-3P (36). Literatures published recently revealed that LINC00858 was found to play key roles in multiple types of cancers including Wilms’ tumor (37), ovarian cancer (38), colon cancer (39), gastric cancer (40), hepatocellular carcinomas (41), and osteosarcomas (42).

Both hsa-miR-182-5p and hsa-miR-183-5p belong to the miR-183/96/182 family because of their sequence homology and function similarity (43). They were found to be among the 10 miRNA signatures that significantly distinguish cancer tissues from adjacent normal tissues from patients with ovarian cancer (44), the five miRNAs associated with tumorigenesis in lung adenocarcinomas (45), the top three miRNAs that target the largest number of genes in atypical endometrial hyperplasia (46), and the five miRNAs that show optimal diagnostic biomarkers for hepatocellular carcinomas (47). The levels of hsa-miR-182-5p are elevated in tumor tissues from patients with ovarian cancer and may predict poor prognosis (48). hsa-miR-182-5p was reported to be associated with resistance of breast cancer cell lines to anti-tumor drug veliparib (49). The significance of hsa-miR-183-5p with additional 10 genes in survival of ccRCC patients has been reported, and high expression of hsa-miR-183-5p predicts a reduced OS rate and poor prognosis (50). It was reported that hsa-miR-182-5p is reduced in renal cancer tissues and cell lines and regulates the expression of DLL4 gene, which causes change of tumor microenvironment and tumor inhibition (51).

PAX3 gene (paired box gene 3) encodes a member of PAX family of transcription factors whose target genes impact proliferation, survival, differentiation, and motility (52). Rhabdomyosarcoma is one of the typical tumors in children and adolescence with a poor prognosis and an OS rate of only 20%–40% (53). PAX3-FOXO1 or PAX3-FKHR, a specific fusion gene that resulted from chromosomal translocation, is an important factor in the occurrence and development of rhabdomyosarcoma (54, 55). Such fusion gene was also reported to be associated with other types of cancers such as melanoma (56) and biphenotypic sinonasal sarcoma, a low-grade spindle cell sarcoma that affects middle-aged adults (57). PAX3 is highly expressed in prostate cancer tissues and cell lines and promotes the progression of prostate cancer by inhibiting the TGF-β/SMAD signal axis (58). Although there is no report about PAX3 in ccRCC, PAX3 gene is differentially hyper-methylated in chromophobe renal cell carcinomas compared with renal oncocytomas, a benign kidney neoplasm (59).

As discussed above, the relation between VHL mutation and patient survival is different in different reports. We identified LINC00942, LINC00858, RP13_392I16.1, hsa-miR-182-5p, hsa-miR-183-5p, and PAX3 as survival-related DE-RNAs between patients carrying wild-type or mutant VHL and were further used to calculate risk scores for each individual patient. The identified RNAs may not directly impact the VHL function. Because each of these DE-RNAs has already been shown to play significant roles in the tumorigenesis and aggressiveness of other types of cancers as discussed above, it is predicted that the risk scores may serve well as factors independent to VHL gene status to predict the malignancy and prognosis of ccRCC patients in the future.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

The studies involving human participants were reviewed and approved by The Ethics Committee of The Fifth Affiliated Hospital of Guangzhou Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

RZ, XL, ZC, LL, and GX have full access to all of the data in the study, take responsibility for the integrity of the data and the accuracy of the data analysis, and wrote the manuscript. SL, YY, YX, DL, HZ, WY, and JB contributed the reagents. All authors contributed to the article and approved the submitted version.

Funding

This study is partially supported by the National Natural Science Foundation of China 81974392 and Natural Science Foundation of Guangdong Province 2018A030313087 to GX, the National Natural Science Foundation of China 82103359 to ZC, the Health Technology Project of Guangzhou City 20201A011106 and 20211A011103 to JB, and the National Natural Science Foundation of China 81772931 to LL.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Clark DJ, Dhanasekaran SM, Petralia F, Pan J, Song X, Hu Y, et al. Integrated Proteogenomic Characterization of Clear Cell Renal Cell Carcinoma. Cell (2020) 180(1):207. doi: 10.1016/j.cell.2019.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ljungberg B, Campbell SC, Choi HY, Jacqmin D, Lee JE, Weikert S, et al. The Epidemiology of Renal Cell Carcinoma. Eur Urol (2011) 60(4):615–21. doi: 10.1016/j.eururo.2011.06.049

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Jiang Z, Chu PG, Woda BA, Liu Q, Balaji KC, Rock KL, et al. Combination of Quantitative IMP3 and Tumor Stage: A New System to Predict Metastasis for Patients With Localized Renal Cell Carcinomas. Clin Cancer Res (2008) 14(17):5579–84. doi: 10.1158/1078-0432.CCR-08-0504

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Hsieh JJ, Purdue MP, Signoretti S, Swanton C, Albiges L, Schmidinger M, et al. Renal Cell Carcinoma. Nat Rev Dis Primers (2017) 3:17009. doi: 10.1038/nrdp.2017.9

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Muglia VF, Prando A. Renal Cell Carcinoma: Histological Classification and Correlation With Imaging Findings. Radiol Bras (2015) 48(3):166–74. doi: 10.1590/0100-3984.2013.1927

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Cancer Genome Atlas Research, N. Comprehensive Molecular Characterization of Clear Cell Renal Cell Carcinoma. Nature (2013) 499(7456):43–9. doi: 10.1038/nature12222

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hsieh JJ, Chen D, Wang PI, Marker M, Redzematovic A, Chen YB, et al. Genomic Biomarkers of a Randomized Trial Comparing First-Line Everolimus and Sunitinib in Patients With Metastatic Renal Cell Carcinoma. Eur Urol (2017) 71(3):405–14. doi: 10.1016/j.eururo.2016.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Haase VH. The VHL Tumor Suppressor in Development and Disease: Functional Studies in Mice by Conditional Gene Targeting. Semin Cell Dev Biol (2005) 16(4-5):564–74. doi: 10.1016/j.semcdb.2005.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Leung JY, Kim WY. Bap1 and Pbrm1: Determinants of Tumor Grade and mTOR Activation in VHL-Deficient Mouse Models of Renal Cell Carcinoma. Cancer Discovery (2017) 7(8):802–4. doi: 10.1158/2159-8290.CD-17-0610

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kaelin WG Jr. The Von Hippel-Lindau Tumor Suppressor Protein and Clear Cell Renal Carcinoma. Clin Cancer Res (2007) 13(2 Pt 2):680s–4s. doi: 10.1158/1078-0432.CCR-06-1865

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Haase VH, Glickman JN, Socolovsky M, Jaenisch R. Vascular Tumors in Livers With Targeted Inactivation of the Von Hippel-Lindau Tumor Suppressor. Proc Natl Acad Sci U.S.A. (2001) 98(4):1583–8. doi: 10.1073/pnas.98.4.1583

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Mandriota SJ, Turner KJ, Davies DR, Murray PG, Morgan NV, Sowter HM, et al. HIF Activation Identifies Early Lesions in VHL Kidneys: Evidence for Site-Specific Tumor Suppressor Function in the Nephron. Cancer Cell (2002) 1(5):459–68. doi: 10.1016/s1535-6108(02)00071-5

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Hogner A, Krause H, Jandrig B, Kasim M, Fuller TF, Schostak M, et al. PBRM1 and VHL Expression Correlate in Human Clear Cell Renal Cell Carcinoma With Differential Association With Patient’s Overall Survival. Urol Oncol (2018) 36(3):94 e1–94 e14. doi: 10.1016/j.urolonc.2017.10.027

CrossRef Full Text | Google Scholar

15. Alves MR, Carneiro FC, Lavorato-Rocha AM, da Costa WH, da Cunha IW, de Cassio Zequi S, et al. Mutational Status of VHL Gene and its Clinical Importance in Renal Clear Cell Carcinoma. Virchows Arch (2014) 465(3):321–30. doi: 10.1007/s00428-014-1629-z

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Dagher J, Kammerer-Jacquet SF, Brunot A, Pladys A, Patard JJ, Bensalah K, et al. Wild-Type VHL Clear Cell Renal Cell Carcinomas Are a Distinct Clinical and Histologic Entity: A 10-Year Follow-Up. Eur Urol Focus (2016) 1(3):284–90. doi: 10.1016/j.euf.2015.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Herman JG, Latif F, Weng Y, Lerman MI, Zbar B, Liu S, et al. Silencing of the VHL Tumor-Suppressor Gene by DNA Methylation in Renal Carcinoma. Proc Natl Acad Sci U.S.A. (1994) 91(21):9700–4. doi: 10.1073/pnas.91.21.9700

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Goodall GJ, Wickramasinghe VO. RNA in Cancer. Nat Rev Cancer (2021) 21(1):22–36. doi: 10.1038/s41568-020-00306-0

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Yuan JH, Yang F, Wang F, Ma JZ, Guo YJ, Tao QF, et al. A Long Noncoding RNA Activated by TGF-Beta Promotes the Invasion-Metastasis Cascade in Hepatocellular Carcinoma. Cancer Cell (2014) 25(5):666–81. doi: 10.1016/j.ccr.2014.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Teng F, Zhang JX, Chen Y, Shen XD, Su C, Guo YJ, et al. LncRNA NKX2-1-AS1 Promotes Tumor Progression and Angiogenesis via Upregulation of SERPINE1 Expression and Activation of the VEGFR-2 Signaling Pathway in Gastric Cancer. Mol Oncol (2021) 15(4):1234–55. doi: 10.1002/1878-0261.12911

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wu C, Chen W, Yu F, Yuan Y, Chen Y, Hurst DR, et al. Long Noncoding RNA HITTERS Protects Oral Squamous Cell Carcinoma Cells From Endoplasmic Reticulum Stress-Induced Apoptosis via Promoting MRE11-RAD50-NBS1 Complex Formation. Adv Sci (Weinh) (2020) 7(22):2002747. doi: 10.1002/advs.202002747

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Miao L, Liu HY, Zhou C, He X. LINC00612 Enhances the Proliferation and Invasion Ability of Bladder Cancer Cells as ceRNA by Sponging miR-590 to Elevate Expression of PHF14. J Exp Clin Cancer Res (2019) 38(1):143. doi: 10.1186/s13046-019-1149-4

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wu X, Xiao Y, Zhou Y, Zhou Z, Yan W. LncRNA FOXP4-AS1 is Activated by PAX5 and Promotes the Growth of Prostate Cancer by Sequestering miR-3184-5p to Upregulate FOXP4. Cell Death Dis (2019) 10(7):472. doi: 10.1038/s41419-019-1699-6

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Yang W, Zhang K, Li L, Xu Y, Ma K, Xie H, et al. Downregulation of lncRNA ZNF582-AS1 Due to DNA Hypermethylation Promotes Clear Cell Renal Cell Carcinoma Growth and Metastasis by Regulating the N(6)-Methyladenosine Modification of MT-Rnr1. J Exp Clin Cancer Res (2021) 40(1):92. doi: 10.1186/s13046-021-01889-8

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gregory PA, Bert AG, Paterson EL, Barry SC, Tsykin A, Farshid G, et al. The miR-200 Family and miR-205 Regulate Epithelial to Mesenchymal Transition by Targeting ZEB1 and SIP1. Nat Cell Biol (2008) 10(5):593–601. doi: 10.1038/ncb1722

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA Hypothesis: The Rosetta Stone of a Hidden RNA Language? Cell (2011) 146(3):353–8. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jiang Y, Zhang H, Li W, Yan Y, Yao X, Gu W. FOXM1-Activated LINC01094 Promotes Clear Cell Renal Cell Carcinoma Development via MicroRNA 224-5p/CHSY1. Mol Cell Biol (2020) 40(3):e00357–19. doi: 10.1128/MCB.00357-19

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Jiang Y, Zhang H, Li W, Yan Y, Yao X, Gu W. LINC01426 Contributes to Clear Cell Renal Cell Carcinoma Progression by Modulating CTBP1/miR-423-5p/FOXM1 Axis via Interacting With IGF2BP1. J Cell Physiol (2021) 236(1):427–39. doi: 10.1002/jcp.29871

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Sun T, Wu Z, Wang X, Wang Y, Hu X, Qin W, et al. LNC942 Promoting METTL14-Mediated M(6)A Methylation in Breast Cancer Cell Proliferation and Progression. Oncogene (2020) 39(31):5358–72. doi: 10.1038/s41388-020-1338-9

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Xu Q, Wang Y, Huang W. Identification of Immune-Related lncRNA Signature for Predicting Immune Checkpoint Blockade and Prognosis in Hepatocellular Carcinoma. Int Immunopharmacol (2021) 92:107333. doi: 10.1016/j.intimp.2020.107333

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang L, Li L, Zhan Y, Wang J, Zhu Z, Zhang X. Identification of Immune-Related lncRNA Signature to Predict Prognosis and Immunotherapeutic Efficiency in Bladder Cancer. Front Oncol (2020) 10:542140. doi: 10.3389/fonc.2020.542140

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hu H, Xu H, Lu F, Zhang J, Xu L, Xu S, et al. Exploring the Effect of Differentially Expressed Long Non-Coding RNAs Driven by Copy Number Variation on Competing Endogenous RNA Network by Mining Lung Adenocarcinoma Data. Front Cell Dev Biol (2020) 8:627436. doi: 10.3389/fcell.2020.627436

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zhu SP, Wang JY, Wang XG, Zhao JP. Long Intergenic non-Protein Coding RNA 00858 Functions as a Competing Endogenous RNA for miR-422a to Facilitate the Cell Growth in non-Small Cell Lung Cancer. Aging (Albany NY) (2017) 9(2):475–86. doi: 10.18632/aging.101171

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Yamada A, Yu P, Lin W, Okugawa Y, Boland CR, Goel A. A RNA-Sequencing Approach for the Identification of Novel Long non-Coding RNA Biomarkers in Colorectal Cancer. Sci Rep (2018) 8(1):575. doi: 10.1038/s41598-017-18407-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhan W, Liao X, Chen Z, Li L, Tian T, Yu L, et al. LINC00858 Promotes Colorectal Cancer by Sponging miR-4766-5p to Regulate PAK2. Cell Biol Toxicol (2020) 36(4):333–47. doi: 10.1007/s10565-019-09506-3

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Sha QK, Chen L, Xi JZ, Song H. Long non-Coding RNA LINC00858 Promotes Cells Proliferation, Migration and Invasion by Acting as a ceRNA of miR-22-3p in Colorectal Cancer. Artif Cells Nanomed Biotechnol (2019) 47(1):1057–66. doi: 10.1080/21691401.2018.1544143

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhou D, Wang J, Xu S, Li Z, Kou D. LINC00858 Facilitates the Malignant Development of Wilms’ Tumor by Targeting miR-653-5p. Minerva Med (2020). doi: 10.23736/S0026-4806.20.06566-0

CrossRef Full Text | Google Scholar

38. Xue H, Wu Z, Rao D, Zhuo B, Chen Q. Long non-Coding RNA LINC00858 Aggravates the Oncogenic Phenotypes of Ovarian Cancer Cells Through miR-134-5p/RAD18 Signaling. Arch Gynecol Obstet (2020) 302(5):1243–54. doi: 10.1007/s00404-020-05722-z

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Xu T, Wu K, Zhang L, Zheng S, Wang X, Zuo H, et al. Long non-Coding RNA LINC00858 Exerts a Tumor-Promoting Role in Colon Cancer via HNF4alpha and WNK2 Regulation. Cell Oncol (Dordr) (2020) 43(2):297–310. doi: 10.1007/s13402-019-00490-8

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Qin L, Yao F, Yu H, Xu Y, Jin Y, Dong W. Up-Regulated lncRNA LINC00858 Facilitates Proliferation and Migration of Gastric Cancer Cells. Pol J Pathol (2020) 71(3):236–43. doi: 10.5114/pjp.2020.99790

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Qi WY, Mao XB, He YB, Xiao CH. Long non-Coding RNA LINC00858 Promotes Cells Proliferation and Invasion Through the miR-153-3p/Rabl3 Axis in Hepatocellular Carcinoma. Eur Rev Med Pharmacol Sci (2020) 24(18):9343–52. doi: 10.26355/eurrev_202009_23017

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Gu Z, Hou Z, Zheng L, Wang X, Wu L, Zhang C. Long Noncoding RNA LINC00858 Promotes Osteosarcoma Through Regulating miR-139-CDK14 Axis. Biochem Biophys Res Commun (2018) 503(2):1134–40. doi: 10.1016/j.bbrc.2018.06.131

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Dambal S, Shah M, Mihelich B, Nonn L. The microRNA-183 Cluster: The Family That Plays Together Stays Together. Nucleic Acids Res (2015) 43(15):7173–88. doi: 10.1093/nar/gkv703

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Wang L, Zhu MJ, Ren AM, Wu HF, Han WM, Tan RY, et al. A ten-microRNA Signature Identified From a Genome-Wide microRNA Expression Profiling in Human Epithelial Ovarian Cancer. PloS One (2014) 9(5):e96472. doi: 10.1371/journal.pone.0096472

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hsu YL, Hung JY, Lee YL, Chen FW, Chang KF, Chang WA, et al. Identification of Novel Gene Expression Signature in Lung Adenocarcinoma by Using Next-Generation Sequencing Data and Bioinformatics Analysis. Oncotarget (2017) 8(62):104831–54. doi: 10.18632/oncotarget.21022

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Tang S, Dai Y. RNA Sequencing Reveals Significant miRNAs in Atypical Endometrial Hyperplasia. Eur J Obstet Gynecol Reprod Biol (2018) 225:129–35. doi: 10.1016/j.ejogrb.2018.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zhao X, Dou J, Cao J, Wang Y, Gao Q, Zeng Q, et al. Uncovering the Potential Differentially Expressed miRNAs as Diagnostic Biomarkers for Hepatocellular Carcinoma Based on Machine Learning in The Cancer Genome Atlas Database. Oncol Rep (2020) 43(6):1771–84. doi: 10.3892/or.2020.7551

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Li Y, Li L. Prognostic Values and Prospective Pathway Signaling of MicroRNA-182 in Ovarian Cancer: A Study Based on Gene Expression Omnibus (GEO) and Bioinformatics Analysis. J Ovarian Res (2019) 12(1):106. doi: 10.1186/s13048-019-0580-7

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Uhr K, Prager-van der Smissen WJC, Heine AAJ, Ozturk B, van Jaarsveld MTM, Boersma AWM, et al. MicroRNAs as Possible Indicators of Drug Sensitivity in Breast Cancer Cell Lines. PloS One (2019) 14(5):e0216400. doi: 10.1371/journal.pone.0216400

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Qin S, Shi X, Wang C, Jin P, Ma F. Transcription Factor and miRNA Interplays Can Manifest the Survival of ccRCC Patients. Cancers (Basel) (2019) 11(11):1668. doi: 10.3390/cancers11111668

CrossRef Full Text | Google Scholar

51. Wang W, Hu W, Wang Y, An Y, Song L, Shang P, et al. Long non-Coding RNA UCA1 Promotes Malignant Phenotypes of Renal Cancer Cells by Modulating the miR-182-5p/DLL4 Axis as a ceRNA. Mol Cancer (2020) 19(1):18. doi: 10.1186/s12943-020-1132-x

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Boudjadi S, Chatterjee B, Sun W, Vemu P, Barr FG. The Expression and Function of PAX3 in Development and Disease. Gene (2018) 666:145–57. doi: 10.1016/j.gene.2018.04.087

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Sultan I, Qaddoumi I, Yaser S, Rodriguez-Galindo C, Ferrari A. Comparing Adult and Pediatric Rhabdomyosarcoma in the Surveillance, Epidemiology and End Results Program, 1973 to 2005: An Analysis of 2,600 Patients. J Clin Oncol (2009) 27(20):3391–7. doi: 10.1200/JCO.2008.19.7483

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Pappo AS, Shapiro DN, Crist WM, Maurer HM. Biology and Therapy of Pediatric Rhabdomyosarcoma. J Clin Oncol (1995) 13(8):2123–39. doi: 10.1200/JCO.1995.13.8.2123

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Crose LE, Galindo KA, Kephart JG, Chen C, Fitamant J, Bardeesy N, et al. Alveolar Rhabdomyosarcoma-Associated PAX3-FOXO1 Promotes Tumorigenesis via Hippo Pathway Suppression. J Clin Invest (2014) 124(1):285–96. doi: 10.1172/JCI67087

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Medic S, Ziman M. PAX3 Across the Spectrum: From Melanoblast to Melanoma. Crit Rev Biochem Mol Biol (2009) 44(2-3):85–97. doi: 10.1080/10409230902755056

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Wong WJ, Lauria A, Hornick JL, Xiao S, Fletcher JA, Marino-Enriquez A. Alternate PAX3-FOXO1 Oncogenic Fusion in Biphenotypic Sinonasal Sarcoma. Genes Chromosomes Cancer (2016) 55(1):25–9. doi: 10.1002/gcc.22295

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Zeng K, Xie W, Huang J, Yang J, Deng K, Luo X. PAX3 Silencing Inhibits Prostate Cancer Progression Through the Suppression of the TGF-Beta/Smad Signaling Axis. Cell Biol Int (2020) 44(10):2131–9. doi: 10.1002/cbin.11421

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Slater AA, Alokail M, Gentle D, Yao M, Kovacs G, Maher ER, et al. DNA Methylation Profiling Distinguishes Histological Subtypes of Renal Cell Carcinoma. Epigenetics (2013) 8(3):252–67. doi: 10.4161/epi.23817

PubMed Abstract | CrossRef Full Text | Google Scholar