Abstract

Background. Clear cell renal cell carcinoma (ccRCC), the most common subtype of renal cell carcinoma (RCC), usually is representative of metastatic heterogeneous neoplasm that links with poor prognosis, but the pathogenesis of ccRCC remains unclear. Currently, numerous evidences prove that long noncoding RNAs (lncRNAs) are considered as competing endogenous RNA (ceRNA) to participate in cellular processes of tumors. Therefore, to investigate the underlying mechanisms of ccRCC, the expression profiles of lncRNAs, miRNAs, and mRNAs were downloaded from the Cancer Genome Atlas (TCGA) database. A total of 1526 differentially expressed lncRNAs (DElncRNAs), 54 DEmiRNAs, and 2352 DEmRNAs were identified. To determine the connection of them, all DElncRNAs were input to the miRcode database. The results indicated that 85 DElncRNAs could connect with 9 DEmiRNAs in relation to our study. Then, databases of TargetScan and miRDB were used to search for targeted genes with reference to DEmiRNAs. The results showed that 203 out of 2352 targeted genes were identified in our TCGA set. Subsequently, ceRNA network was constructed according to Cytoscape and the targeted genes were functionally analyzed to elucidate the mechanisms of DEmRNAs. The results of survival analysis and regression analysis indicated that 6 DElncRNAs named COL18A1-AS1, WT1-AS, LINC00443, TCL6, AL356356.1, and SLC25A5-AS1 were significantly correlative with the clinical traits of ccRCC patients and could be served as predictors for ccRCC. Finally, these findings were validated by quantitative RT-PCR (qRT-PCR). Based on these discoveries, we believe that this identified ceRNA network will provide a novel perspective to elucidate ccRCC pathogenesis.

1. Introduction

Renal cell carcinoma (RCC), the 7th most common malignant tumor in male, was estimated at almost 63,000 newly cases and 14,000 deaths occurring in America in 2016[1]. Epidemiological researches reported that more than half of renal cell carcinomas have local invasive tumors and 17% of patients with distant metastasis at the time of diagnosis [2]. Clear cell renal cell carcinoma (ccRCC), the most frequent RCC subtype, accounts for approximately 90% of RCC [3]. Recently, a series of targeted therapies have improved ccRCCs median overall survival rate, while lots of patients ultimately become resistant to these advances due to the limited effects [4]. Although several evidences have documented that somatic genetic mutations were linked with the formation and progress of ccRCC [57], it is imperative to identify more effective biomarkers that could help us uncover more thorough pathogenesis regarding the biological characteristics of ccRCC.

With the development of next-generation sequencing technology to human whole genomes, accumulating analyses suggest that genomes encode proteins finally are less than 2% whereas more than 75% are transcribed into noncoding RNAs [8]. Long noncoding RNAs (lncRNAs), longer than 200 nucleotides in length accounting for most of these noncoding RNAs, have been verified as a new class of dramatically role involved in transcriptional and posttranscriptional regulation during various type of cancer progression [9, 10]. In 2011, the competing endogenous RNA (ceRNA) hypothesis proposed by Salmena et al. [11], reported that miRNA could suppress targeted gene transcript through complementation with the mRNA sequence at the posttranscriptional level, leading to a variety of RNAs silence. But lncRNA with partial similar miRNA sequence could function as a miRNA endogenous sponge to affect stability of binding with targeting mRNA via shared miRNA response elements, which might exert a significant effect in pathogenesis of cancer. Since then, this hypothesis has been emerged and validated by several experimental studies. Lv et al. identified an epithelial-to-mesenchymal transition lncRNA H19, which could sequester miR-29b-3p from regulating targeted gene DNMT3B by acting as a ceRNA [12]. Recently, a new LincRNA SNHG14 was upregulated and facilitated ccRCC cell invasion via sponging miR-203 [13]. However, to the best knowledge of the authors, little is known concerning the investigation of ceRNA related ccRCC based on high-throughput detection with large-scale samples.

The Cancer Genome Atlas (TCGA), a public available platform with more than 30 types of cancer from 11,000 patients at least and their clinicopathological information, has been widely applied by large numbers of researches to explore the genetic basis of tumor according to the high-throughput sequencing [14]. In our research, RNA expression profiles between 539 tumor samples and 72 matched normal samples of ccRCC were downloaded from TCGA project. As a result, 1526 aberrantly expressed lncRNAs were identified according soft R, and the miRcode, a transcriptome-wide microRNA target database was then used to screen the hub differential expressed lncRNAs (DElncRNAs) and DEmiRNAs as well as DElncRNAs-DEmiRNAs interactions. Moreover, the DEmiRNA targeted DEmRNAs and DEmiRNAs-DEmRNAs connections were screened performing with the two popular databases. CeRNA network was then constructed to elucidate the crosstalk of the hub DElncRNAs, DEmiRNAs, and DEmRNAs. Survival analysis and regression analysis were utilized to analyze the correlation between clinical traits of patients and related lncRNA. Finally, these bioinformatics discoveries were validated by quantitative RT-PCR (qRT-PCR).

2. Materials and Methods

2.1. Acquisition of Clinical Characteristics and RNAseq Data

A total of 539 ccRCC patient information were obtained from TCGA database. After excluding patients without entire follow-up information, 530 patients remained available for survival analysis. Besides, lncRNA, miRNA and mRNA expression files included tumor and corresponding normal samples were also downloaded from TCGA portal (https://cancergenome.nih.gov). Moreover, 30 frozen tissue specimens contained 15 tumor tissues and 15 matched adjacent nontumor tissues were obtained from ccRCC patients in Renmin Hospital of Wuhan University. All tissues were collected immediately after surgical resection, and snap-frozen in liquid nitrogen until RNA extraction. This study was approved by the ethics committee of Renmin Hospital of Wuhan University.

2.2. Data Procession of RNAseq

An EdgeR (empirical analysis of digital gene expression data in R) package [15], a package for examining differential expression of replicated count data based on Poisson mode, was used to identify the RNAseq and miRNAseq data of ccRCC. Only the differential expression genes (DEGs) with adj.P.Val < 0.05 and |log2 fold change (FC)| ≥ 2 were considered as significant.

2.3. Construction of lncRNA-miRNA-mRNA ceRNA Network

The ceRNA network was constructed according to the relationship between lncRNA, miRNA and mRNA. Firstly, aberrantly expressed lncRNA, miRNA, and mRNA were identified, and then lncRNA-miRNA connections were predicted by miRcode [16] (http://www.mircode.org/), a database with more than 10,000 lncRNAs which offers whole transcriptome human miRNA target predictions at the basis of comprehensive GENCODE (http://www.gencodegenes.org/) gene annotation. Subsequently, two famous online databases for miRNA targeted mRNA prediction and functional annotation, TargetScan [17] (http://www.targetscan.org/) and miRDB [18] (http://www.mirdb.org/), were used to predict the miRNA-mRNA interactions. Each targeted mRNA in the present research was then selected to intersect with DEmRNAs. Finally, Cytoscape software [19] (version 3.5.1) was constructed the ceRNA network.

2.4. Functional Enrichment Analysis of Targeted DEmRNAs

Database for annotation, visualization, and integrated discovery (DAVID), an online tool, has been used as functional enrichment analysis in many studies. In this paper, DAVID was used for gene ontology (GO) analysis to understand the functions of targeted genes in terms of biological process (BP), cellular component (CC), and molecular function (MF). In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) and an R Package, clusterProfiler [20], were distinguished pathway enrichment of each targeted gene. Cutoff value was set as P value < 0.05.

2.5. Quantitative Real-Time PCR (qRT-PCR)

Total RNA from specimens were extracted using TRIzol reagent (Invitrogen, Paisley, UK) and then reverse transcribed to cDNA by performing reverse transcription enzyme system (Thermo Fisher Scientific, Waltham, USA) according to the manufacturer instructions. RT-qPCR was utilized to detect selected lncRNA expression performing with SYBR Green qRT-PCR Master Mix (TaKaRa, Shiga, Japan). 2−ΔΔCt method was used to calculate the results of each lncRNA expression. All primer sequences are listed in Table 1, and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) is designed as an internal control to normalize the lncRNA expression levels.

2.6. Statistical Analysis

All results in this paper were expressed as the mean ± SD. The significance of the differences between tumor and normal group was performed with Student’s t-test, and the value of was considered significantly. The relationship between dysregulated lncRNAs expression level and patient’s survival was evaluated by univariate and multivariate Cox regression analysis. In the predictive model, a formula was calculated to predict survival in the training set as previously described [21, 22]. Each patient was assigned a risk score in the light of the Cox regression coefficients of significant survival lncRNAs, and divided into low-risk or high-risk groups based on the median risk score. The Kaplan-Meier (KM) method was used to estimate overall survival, and Receiver operating characteristic (ROC) analysis was performed to compare the sensitivity and specificity of survival prediction by R Soft (version 3.4.3) [23].

3. Results

3.1. Identification of DElncRNAs, DEmiRNAs, and DEmRNAs in ccRCC Patients

To determine the lncRNA, miRNA, and mRNA expression levels in ccRCC, 539 ccRCC patients and 72 normal samples were downloaded from TCGA project. Compared with normal patients, there were 1526 DElncRNAs, 54 DEmiRNAs, and 2352 DEmRNAs, respectively, among which 1083 lncRNAs, 33 miRNAs, and 1599 mRNAs were upexpressed as well as 443 lncRNAs, 21 miRNAs, and 753 mRNAs were downexpressed. To better distinguish the features of these RNAs, clustering analysis and volcano plots were constructed between tumor group and normal group (Figure 1). Each heatmap (Figures 1(a)1(c)) and volcano plot (Figures 1(d) and 1(e)) showed a distinct separation of distribution.

3.2. Prediction of DEmiRNA Targeting DElncRNAs and DEmRNA and Construction of ceRNA Network

To clarify the correlation between miRNA targeting lncRNAs in the ceRNA crosstalk, all DElncRNAs were input to the miRcode database. The results indicated that 85 DElncRNAs could interact with 86 miRNAs. However, only 9 out of 86 miRNAs were discerned in relation to our study (Table S1). We then performed with the databases of TargetScan and miRDB to search for targeted genes with reference to DEmiRNAs. The results showed that 203 out of 2352 targeted genes were identified in the TCGA set (Fig S1). Based on the data as shown above, Cytoscape software was applied to construct and plot lncRNA-miRNA-mRNA ceRNA crosstalk. A total of 85 lncRNAs, 9 miRNAs, and 203 mRNAs were involved in the network (Figure 2).

3.3. GO and KEGG Enrichment Analysis of Targeted mRNAs in ceRNA

To better comprehend the overall biological features and the signal pathways in relation to the ceRNA network, these mRNAs were investigated by DAVID database and a package, clusterProfiler in R software. As shown in Table 2 and Figure 3(a), these mRNAs mainly were enriched in sodium ion transport, cell-cell signaling, and ion transport about BP, and CC terms related to extracellular region part, plasma membrane, and plasma membrane part, as well as MF terms linked to sodium ion binding, voltage-gated sodium channel activity, and alkali metal ion binding. Pathway enrichment analysis suggested that these mRNAs were predominantly enriched in metabolic related pathway, cytokine-cytokine receptor interaction, and other important pathways (Figure 3(b)).

3.4. Evaluation of the Significant lncRNAs in ceRNA concerning ccRCC Patients

To determine whether those lncRNAs were correlative with the survival of ccRCC patients, we used KM method and univariate Cox proportional hazard regression analysis to test and screen the important lncRNAs linked with survival of the ccRCC patients in the R software. Among 85 aberrantly expressed lncRNAs, 34 lncRNAs with P value less than 0.05 were presented by the survival curve (Figure S2). In addition, multivariate Cox regression analysis was carried out to evaluate whether each of these candidates were acted as independent predictors of ccRCC patient’s survival. As indicated in Table 3, 6 out of 34 lncRNAs named COL18A1-AS1, WT1-AS, LINC00443, TCL6, AL356356.1, and SLC25A5-AS1 were served as covariates () in the training set. With the risk score formula as mentioned above, the 6-lncRNA expressed levels were then calculated for each patient in the training set, and patients with risk score were then divided into low-risk groups (n = 265) and high-risk (n = 265) groups via using the median risk score of the test as a cutoff point. The plot results indicated that patients with high-risk had significantly shorter overall survival than those in low-risk group (, Figure 4(a)), and ROC analysis revealed that area under receiver operating characteristic (AUC) curve of 6 significant lncRNAs was 0.716 (Figure 4(b)), suggesting that these lncRNAs might be functionally important in pathogenesis of ccRCC patients.

3.5. Validation of lncRNAs through qRT-PCR

In terms of the above analysis from microarray analysis, these 6 lncRNAs (COL18A1-AS1, WT1-AS, LINC00443, TCL6, AL356356.1, and SLC25A5-AS1) involving ceRNA network were regarded as potential biomarkers linked with ccRCC patients (Table 4). To verify the validity and reliability of the results, the 6 lncRNAs were selected to analyze their expression between 15 diagnosed ccRCC specimens and 15 adjacent nontumor samples through qRT-PCR. As indicated in Figure 5, the results of qRT-PCR were consistent with our findings in computer analysis, implying that these lncRNAs may be acted as biomarkers in the process of ccRCC deterioration.

4. Discussion

RCC, a heterogeneous kind of tumors, mainly originates in renal tubular epithelial cells, more common in males after the age of 60 [1]. It has been documented that ccRCC as the most common pathological subtype of RCC is still leading to a high number of mortality in urologic neoplasms [24]. Therefore, in order to improve this situation, it is important to elucidate the ccRCC related genes and underlying molecular mechanisms involving in the progression of ccRCC.

Nowadays, increasing studies have showed the crucial biological functions of lncRNAs involved in the development and process of human diseases by epigenetic regulation and transcription as well as posttranscription regulation. A growing number of evidences have revealed that lncRNAs harbored similar sequences to their targeted miRNAs to regulate mRNAs expression in many tumor [2527]. LncRNA ROR promotes radioresistance through acting as a ceRNA for microRNA145 to regulate RAD18 expression in hepatocellular carcinoma. LncRNA Gas5 serves as a ceRNA to activate PTEN/AKT pathway via sponging miR-222-3p in papillary thyroid carcinoma. Moreover, lncRNA PTENP1 could interact with miR-19b, upregulate the expression of p53, and downregulate p-AKT expression, inhibiting breast cancer cell growth and metastasis.

However, few researches between lncRNA and ccRCC have been achieved based on the large numbers of samples. In this regard, we comprehensively integrated 611 samples lncRNA, miRNA, and mRNA data from TCGA and investigated the lncRNAs related ceRNA network to uncover the underlying mechanism about ccRCC. In the present research, 34 DElncRNAs were identified in ccRCC samples compared with the normal samples. 6 out of them involved in the ceRNA network were significantly linked with overall survival might be considered as prognostic markers for ccRCC. Notably, two DElncRNAs, lncRNA WT1-AS and TCL6, interacted with most of the DEmiRNAs in ceRNA network. LncRNA WT1-AS, Wilms tumor 1 (WT1) gene antisense, locates in the upstream of WT1, which was originally found on the chromosome 11p and function as regulating in the development of paediatric kidney tumor originated from embryonal blastemal stem cells [28]. Kaneuchi et al. [29] found that the expression of WT1-AS was in parallel with WT1 gene in ovarian clear cell adenocarcinoma with poor prognosis. Moreover, WT1-AS could upregulate expression of WT1 to inhibit hepatocellular carcinoma cell proliferation and apoptosis via activating JAK2/STAT3 and MAPK signaling pathway [30]. Several studies have demonstrated that WT1-AS could be acted as a therapeutic target in the acute myeloid leukemia and gastric cancer as well as colon cancer [28, 31, 32]. Consistent with the researches above, lncRNA WT1-AS was also found aberrant expression in the 530 ccRCC samples compared with normal sample and the results of univariate and multivariate regression analysis showed that lncRNA WT1-AS could be acted as an independent predictor of ccRCC patients. Moreover, high expression of lncRNA WT1-AS was related to poor prognosis based on KM curve. Of interest is the fact that the results of RT-qPCR were in agreement with the bioinformatic analysis, confirming that this lncRNA could play an important role in the process of ccRCC. In addition, the result of ceRNA network showed that the targeted miRNAs of lncRNA WT1-AS were miR-141 miR-155 and miR-216b, suggesting that it may be acted as a sponge to compete with these three key DEmiRNAs participated in the pathogenesis mechanism of ccRCC carcinogenesis.

Another lncRNA TCL6, named T-cell leukemia/lymphoma 6 (TCL6), from a region on chromosome 14q locus was found in T-cell leukemia. The TCL6 gene could encode different open reading frames with no homology via through alternative splicing, leading to T-cell lymphoproliferative diseases [33]. Previous study has showed that overexpressed of TCL6 could suppress the proliferation and growth of ccRCC cells and the expression of TCL6 was negatively linked with TNM stage of tumor [34], which was similar to our results. The samples in our research were much more than their study with 71 tumor samples, Moreover, this variable in our univariate and multivariate regression analysis were both significantly lncRNAs with P value less than 0.05, and the results of ROC curve confirmed that lncRNA TCL6 could be considered as an independent predictor for ccRCC patients. However, they proved the hypothesis via some cellular functional experiments instead of further mechanism discussion. Our results showed that lncRNA TCL6 could be served as a ceRNA to in regulation the expression of targeted genes by sponging miR-210, miR-216b, and miR-122, further reinforcing their conclusion. Furthermore, other lncRNAs (COL18A1-AS1, LINC00443, AL356356.1, and SLC25A5-AS1) were also significantly aberrantly expressed in the ccRCC tumors compared with the normal samples, but there are scarcely reports about them until now. Therefore, further investigations are needed to identify the function and exact mechanism of these for ccRCC.

MiRNA still warrants to be considered as a crucial modulator in the regulation of tumors although lncRNA has obtained lots of attentions in these years. There is no doubt that targeted genes could be changed via interrupting miRNAs, resulting the aggressive tumorigenic state. In our research, the intersection results between two important lncRNA WT1-AS and TCL6 targeted miRNA involved ceRNA network showed that miR-216b played an important role for ccRCC. MiR-216b was first discovered to function as a tumor suppressor involving in regulating nasopharyngeal cell proliferation and invasion as well as tumor growth via inhibiting the KRAS/AKT and ERK pathways [35]. Subsequently, several researches further proved that miR-216b played significant in various tumors such as breast cancer, colorectal cancer, and hepatocellular carcinoma as well as pancreatic ductal adenocarcinoma [3639]. Consistent with previous researches, our results also showed that the expression of miR-216b in tumor samples was abnormal and it participated in the ceRNA network for ccRCC. In addition, the targeted genes based on the online databases analysis mainly involved in metabolic pathway and cytokine-cytokine receptor interaction, and the most important functional alterations were sodium ion transport and extracellular region part, as well as sodium ion binding. These findings suggesting that miR-216b might influence targeted genes via sodium ion transporting from extracellular position in ccRCC carcinogenesis.

Nevertheless, our study has several limitations. Due to all the samples and patient clinical information in our study were based on TCGA, another independent database includes a large-scale multicenter clinical samples should be further verified our findings. In addition, the mechanism of ceRNA network suggests that one lncRNA can connect with numbers of miRNAs and relevant targeted mRNAs [11]. Nevertheless, the miRNAs and targeted mRNAs involved in ceRNA network were based on TCGA database in our research. In this regard, our further research will focus on validating these miRNAs and mRNAs and uncover the mechanism of ceRNA in vivo and in vitro to prove our discovery, especially to confirm the crucial interactions we discussed above.

In conclusion, dysregulated 34 cancer-specific lncRNAs were screened between tumor and nontumor samples from TCGA database. Moreover, we successfully constructed lncRNA-associated ceRNA network by integrating aberrant lncRNA, miRNA, and mRNA expression profiles. In addition, six lncRNAs were identified to be remarkably correlative with overall survival of patients and could be considered as candidate predictors. Obviously, our findings will contribute to further understanding of ccRCC pathogenesis and bring to light a novel ceRNA regulatory mechanism in ccRCC.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Ethical Approval

The proposal and consent procedures of this study were approved by the Ethic Committee of Renmin Hospital of Wuhan University.

All patients provided their verbal informed consent.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Authors’ Contributions

Kang Yang and Xiao-fan Lu participated in the design of this study. Kang Yang and Xiao-fan Lu performed the statistical analysis. Kang Yang and Peng-cheng Luo carried out the study and collected important background information. Kang Yang wrote the article. Pengcheng-Luo and Jie Zhang conceived of this study and helped to draft the manuscript. Kang Yang, Xiao-fan Lu, Peng-cheng Luo, and Jie Zhang read and approved the final manuscript.

Acknowledgments

The authors would like to thank Dr. Yusha Xiao, Dr. Cong Qin, Dr. Yang Su, and Dr. Cheng-cheng Xiao for helpful advice. This work was funded by the National Natural Science Foundation of China (no. 81470923, no. 81770688, and no. 81770078).

Supplementary Materials

Supplementary 1. Figure S1: Venn diagrams of differential expressed mRNAs interacted with putative targeted mRNAs predicted by databases of TargetScan (v7.1), miRDB.

Supplementary 2. Figure S2: Kaplan-Meier (KM) analyses on 34 significantly lncRNAs with value less than 0.05. Survival curves for patients in two groups, solid red lines represented high expression and green lines denoted low expression of each lncRNA.