Identification of RNA Transcript Makers Associated With Prognosis of Kidney Renal Clear Cell Carcinoma by a Competing Endogenous RNA Network Analysis

Objective This study aims to identify several RNA transcripts associated with the prognosis of kidney renal clear cell carcinoma (KIRC). Methods The differentially expressed mRNAs, lncRNAs, and miRNAs (DEmRNAs, DElncRNAs, and DEmiRNAs) between KIRC cases and controls were screened based on an RNA-seq dataset from The Cancer Genome Atlas (TCGA) database. Subsequently, miRcode, miRDB, and TargetScan database were used to predict interactions between lncRNAs, miRNAs and target mRNAs. Then, a ceRNA network was built using miRNAs-mRNAs and lncRNAs-miRNAs pairs. Functional analysis of mRNAs in ceRNA was performed. Finally, the survival analysis of RNA transcripts in ceRNA network and correlation analysis for key RNA regulators were carried out. Results There were 1527 DElncRNAs, 54 DEmiRNAs, and 2321 DEmRNAs. A ceRNA network was constructed among 81 lncRNAs, 9 miRNAs, and 197 mRNAs. Functional analysis showed that numerous mRNAs were significantly associated with regulation of cellular glucuronidation. In addition, 35 lncRNAs, 84 mRNAs and two miRNAs were significantly corelated to the survival of patients with KIRC (P < 0.05). Among them, miRNA-21 and miRNA-155 were negatively related to three lncRNAs (LINC00472, SLC25A5.AS1, and TCL6). Seven mRNA targets of miRNA-21 (FASLG, FGF1, TGFBI, ALX1, SLC30A10, ADCY2, and ABAT) and 12 mRNAs targets of miRNA-155 (STXBP5L, SCG2, SPI1, C12orf40, TYRP1, CTHRC1, TDO2, PTPRQ, TRPM8, ERMP1, CD36, and ST9SIA4) also acted as prognostic biomarkers for KIRC patients. Conclusion We screened numerous novel prognosis-related RNA markers for KIRC patients by a ceRNA network analysis, providing deeper understandings of prognostic values of RNA transcripts for KIRC.


INTRODUCTION
Kidney renal clear cell carcinoma (KIRC) is a major subtype of renal carcinoma and originated from renal epitheliums (Drucker, 2005). Numerous studies have estimated that the morbidity and mortality of KIRC have been gradually increasing in recent years, accounting for most cancer-related deaths (Hsieh et al., 2017;Wong et al., 2017). Clinically, surgical resection has been a powerful therapeutic option for localized KIRC. Unfortunately, approximately 30% patients with localized KIRC ultimately develop malignant metastases due to high recurrence and delayed diagnosis (Hutson and Figlin, 2007;Gupta et al., 2008). Moreover, metastatic KIRC is significantly resistant to radiotherapy and chemotherapy, resulting in a poor clinical prognosis (Zall et al., 2010;De Felice and Tombolini, 2018). Therefore, it is an urgent need to identify several novel prognostic makers, which will contribute to developing effective therapeutic strategies for improving the overall survival (OS) of patients undergoing KIRC.
Previous studies have reported that many non-coding RNA transcripts such as long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) were involved in molecular pathogenesis of several types of cancer, including KIRC (Khadirnaikar et al., 2019;Ye et al., 2019). More remarkably, an increasingly large number of research groups revealed that the interplays among messenger RNAs (mRNAs), lncRNA and miRNAs played essential roles in coordinating numerous complicated molecular processes, however, the dysregulation of biological pathways will lead to pathological conditions (Karreth and Pandolfi, 2013). Specifically, lncRNAs as molecular sponges can bind to miRNAs by miRNA response elements (MREs) and consequently inhibit the expression of mRNAs. Accordingly, a competitive endogenous (ceRNA) network theory is proposed and discussed the interactions among these three types of RNA transcripts (Salmena et al., 2011). Recently, the researchers have screened many RNA transcripts targets associated with the clinical outcomes of various cancers by ceRNA regulatory network analyses. For example, Wang et al. (2019) constructed a mRNA-miRNA-lncRNA network, in which several RNA transcripts were associated with the prognosis of patients with pancreatic cancer. Similarly, Lin et al. (2018) also established a ceRNA network by a genome-wide analysis of lncRNAs, miRNAs and mRNAs, which were all related to the survival of hepatocellular carcinoma patients. Moreover, they found that this network was linked with a poor prognosis by regulating the cell cycle . However, few studies considered to identify prognostic RNA transcripts for KIRC patients based on a lncRNA-miRNA-mRNA ceRNA network analysis.
In this study, we conducted a differential expression analysis of three types of RNA transcripts (lncRNAs, mRNAs and miRNAs) between KIRC patients and normal controls using an RNA-seq dataset from The Cancer Genome Atlas (TCGA) database. Then, a lncRNA-miRNA-mRNA ceRNA regulatory network was established to the characterize RNA transcripts. Finally, a survival analysis was carried out to evaluate the relationships of RNA transcripts and KIRC prognosis, which will provide new insights into prognostic values of RNA transcripts for KIRC patients.

MATERIALS AND METHODS
Acquisition of RNA-Sequencing (RNA-Seq) Data RNA-seq raw data from patients with KIRC was downloaded from TCGA database 1 (which consisted of three types of RNA (lncRNA, miRNA and mRNA). This dataset contained 611 samples (505 samples from KIRC patients and 106 sample from healthy individuals). Those KIRC patients with missing followup information were excluded from this study. Finally, RNA data and clinical data from 474 tumor samples and 106 normal samples were analyzed.

Data Pre-processing and Differential Expression Analysis
The EdgeR package in R software has been successfully used for differential expression analysis of RNA-seq based on a negative binomial distribution model 14. The raw read counts of all samples were merged in a single read count matrix. The matrix was normalized and standardized by TMM (Trimmed Mean of M-values) method before the differential expression analysis. The TMM normalization method is implemented in the edgeR package by means of the calcNormFactors function. Herein, we utilized the EdgeR package to identify the differentially expressed mRNAs (DEmRNAs), lncRNAs (DElncRNAs), and miRNAs (DEmiRNAs) between KIRC tissues and normal tissues. The corresponding P-value of differential gene expression was calculated by Student's t-test and adjusted by Benjamini and Hochberg method. The screening criteria for statistical significance were set as adjusted P < 0.05 and |log 2 fold change (FC)|>2. Besides, the Volcano plots were visualized to demonstrate the distribution of differential expression of DEmRNAs, DElncRNAs and DEmiRNAs using the ggplot2 packages in R.

Construction of lncRNA-miRNA-mRNA Network
Existing evidence has suggested that there were crosstalks among several types of RNA transcripts by MREs, including lncRNAs, miRNAs and mRNAs. Accordingly, numerous studied have demonstrated that lncRNAs act as ceRNAs sponge to bind to miRNAs through MREs, thereby regulating the expression level of mRNA (Tay et al., 2014). Therefore, a ceRNA regulatory network was established to investigate associations among DElncRNAs, DEmiRNAs, and DEmRNAs. LncRNA-miRNA-mRNA pairs should share the same miRNAs. Briefly, the miRcode 2 , a database including more than 10,000 lncRNAs and providing putative miRNA target sites on the basis of comprehensive GENCODE 3 gene annotation, was employed to predict DElncRNA-DEmiRNA pairs (Jeggari et al., 2012). Meanwhile, the interactions between DEmiRNAs and mRNAs were assessed using experimentally validated miRTarBase 4 and TargetScan 5 database (Tan et al., 2009;Hsu et al., 2014). The interactions that matched with DElncRNAs and DE-mRNAs were screened. Also, miRNA-lncRNA and miRNA-mRNA interactions that did not contain the same miRNA were eliminated. Subsequently, the intersecting mRNAs between predicted mRNA targets for DEmiRNAs and DEmRNAs were obtained. Finally, we used Cytoscape software 3.5.1 6 to construct and visualize the ceRNA network.

Functional Enrichment Analysis of DEmRNAs in ceRNA Network
To explore the underlying biological functions of DEmRNAs from the ceRNA regulatory network, the Gene Ontology (GO) functional annotation analysis was carried out by an web-based functional annotation tool, Database for Annotation, Visualization and Integrated Discovery (DAVID) 7  (Ashburner et al., 2000;Dennis et al., 2003). GO analysis involved in three categories of biological process (BP), cellular component (CC), and molecular function (MF). In addition, the clusterProfiler package in R was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEmRNAs in ceRNA network (Yu et al., 2012;Kanehisa et al., 2017). The significant GO terms and KEGG pathways were selected according to the cut of P < 0.05.

Survival Analyses
The correlations between DEmRNAs, DElncRNAs and DEmiRNAs of the ceRNA regulatory network and OS of patients suffering from KIRC were, respectively, assessed by Kaplan-Meier (KM) method and log-rank test using R Survival package 8 . The P < 0.05 was regarded as statistically significantly different. The DElncRNAs, DEmiRNAs and DEmRNAs that showed close association with 5-year survival were considered as key prognostic markers for KIRC patients.

Correlation Analysis
Increasing studies have reported that lncRNAs-miRNAs signatures were predominately correlated with the clinical prognosis of cancer patients. Moreover, lncRNAs participate   out by applying the BeyoRT II cDNA synthesis kits (Beyotime Biotechnology, shanghai, China). The cDNAs were confirmed to qPCR experiments with utilizing the SYBR Green qRCR kits (TaKaRa Bio Inc., Otsu Prince Hotel, Japan) according to the manufacturer's recommendations.

Cell Proliferation Assays
The cell growth was detected in accordance with CCK-8 kit (DOJINDO, kyushu, Japan). Briefly, 786O cells were placed into ninety-six well plates (2000 cells per well). At assigned timepoint, CCK8 reagents were added into the wells (10 µl per well). After incubated at 37 • C in 5 c /o CO 2 for 3.5 h, the absorbance was detected at 450 nm on a microplate spectrophotometer.  and 2321 DEmRNAs (1557 up-regulated mRNAs and 764 down-regulated mRNAs; Table 3) between KIRC group and control group were extracted according to pre-determined screening criteria. Volcano plots for DElncRNAs, DEmiRNAs, and DEmRNAs were visualized and showed their distributions (Figure 1). All of DElncRNAs, DEmiRNAs, and DEmRNAs were list in Supplementary Table S1.

Construction of ceRNA Network
The predictive analysis for interactions between DElncRNAs and DEmiRNAs revealed that there were 174 DElncRNAs-DEmiRNAs pairs, including 81 DElncRNAs and 9 DEmiRNAs. Furthermore, two databases (miRTarBase and TargetScan) were utilized to predict the mRNA targets of 9 DEmiRNAs and a total of 197 DEmRNAs were intersected with predicted mRNAs. Afterward, 212 pairs of DEmiRNAs and DEmRNAs were determined, including 197 DEmRNAs and 9 DEmiRNAs. Ultimately, the DElncRNAs-DEmiRNAs-DEmRNAs ceRNA network was constructed using Cytoscape 3.5.1, which contained 81 DElncRNAs, 9 DEmiRNAs and 197 DEmRNAs (Supplementary Table S2 and Figure 2).

Functional Analyses
The GO and KEGG functional analysis were carried out to understand potential roles of DEmRNAs on KIRC. We found that these mRNAs were involved in 64 GO-BP terms, 17 GO-CC terms and 19 GO-MF terms. Specifically, the top three GO-BP terms were xenobiotic glucuronidation, flavonoid biosynthetic process and negative regulation of cellular glucuronidation (Supplementary Table S3). Meanwhile, the overwhelming majority of mRNAs were predominately correlated with three GO-CC terms, including integral component of membrane, proteinaceous extracellular matrix and extracellular space (Supplementary Table S3). For GO-MF analysis, the mRNAs played essential roles in multiple significantly enriched terms, such as glucuronosyltransferase activity and retinoic acid binding (Supplementary Table S3). In addition, KEGG enrichment analysis suggested that the mRNAs were markedly enriched in 16 KEGG pathways, including pentose and glucuronate interconversions, ascorbate and aldarate metabolism and retinol metabolism pathway (Supplementary Table S3).

Survival Analysis
To further assess the potential prognostic values of lncRNAs, miRNAs, and mRNAs in ceRNA, the corresponding KM survival analyses were conducted. The results indicated that 35 lncRNAs, 84 mRNAs and two miRNAs (up-regulated miRNA-21 and miRNA-155) were dramatically related to 5-year survival of patients undergoing KIRC (P < 0.05). More notably, lncRNAs participate in molecular mechanisms of cancerogenesis through sponging miRNAs, in which lncRNAs are generally negatively associated with miRNAs. Herein, we performed a Pearson's correlation coefficient analysis between prognosis-related lncRNAs and miRNAs. The results revealed that there were significantly negative correlations between miRNA-21 and 14 lncRNAs. Moreover, miRNA-21-LINC00472 (P = −0.337), miRNA-21-SLC25A5.AS1 (P = −0.330) and miRNA-21-TCL6 (P = −0.284) were top three pairs with the strongest relationships (Supplementary Table S4 and Figure 3). For miRNA-155, we noted that it was inversely related to 10 lncRNAs (Supplementary Table S4).

Risk Evaluation of the Key lncRNA and miRNA With Overall Survival of KIRC Patients
Many previous studies suggested that LINC00472 and miRNA-21 served as significant predictors for the survival of KIRC patients (Zaman et al., 2012;Yang et al., 2018;Yu et al., 2018). We assessed the associations between key RNA transcripts (SLC25A5-AS1 and miRNA-21) and the clinical survival of patients with KIRC. The corresponding P-values with 95% CI were calculated to predict patients' survival risk. We found that SLC25A5-AS1 and miRNA-21 were strongly related to the risk of KIRC occurrence according to the univariate Cox analysis (SLC25A5-AS1, HR = 0.610, 95% CI: 0.438-0.848, P = 0.003; miRNA-21, HR = 1.908, 95% CI: 1.37-2.657, P < 0.001; Table 4). Furthermore, the multivariate Cox analysis indicated that SLC25A5-AS1 was the independent prognostic factor of KIRC (HR = 0.836, 95% CI: 0.587-1.190, P = 0.032; Table 4). Besides, the relationships between clinicopathological features and KIRC risk were also evaluated. Our findings showed that age, pathological stage and cancer status were dramatically associated with the risk of KIRC patients (Table 4). Afterward, the correlations between SLC25A5-AS1 and the clinical features of patients with KIRC were also evaluated. The results showed that SLC25A5-AS1 was related to clinical stage, pathological stage and tumor volume (Table 5). Accordingly, we observed that there was the highest expression of SLC25A5-AS1 in KIRC patients under pathological stage IV, suggesting that SLC25A5-AS1 may be a key target for predicting KIRC progression (Figure 6).

The Expression of SLC25A5-AS1, LINC00472, and TCL6 in KIRC Specimens
To verify the bioinformatics analysis, we performed qRT-PCR assays to evaluate the SLC25A5-AS1, LINC00472, and TCL6 levels in KIRC samples and adjacent normal tissues. The results illuminated that the SLC25A5-AS1 levels were downregulated in KIRC samples in comparison to the adjacent normal tissues ( Figure 7A). Nevertheless, LINC00472 and TCL6 were up-regulated in KIRC samples compared with adjacent normal tissues (Figures 7B,C).
Knocking Out SLC25A5-AS1, LINC00472, and TCL6 Will Inhibit or Promote the Proliferation of KIRC Cells To demonstrate the effect of SLC25A5-AS1, LINC00472, and TCL6 on cellular growth of KIRC cells, we carried out loss of function studies using SLC25A5-AS1, LINC00472, and TCL6 siRNAs. Real Time-PCR analysis showed that SLC25A5-AS1 siRNAs, LINC00472 siRNAs, and TCL6 siRNAs, which constructed si-lnc#1 and si-lnc#2 respectively, successfully knocked out the expression of SLC25A5-AS1, LINC00472, and TCL6 in KIRC cells ( Figure 8A). Subsequently, the results from CCK-8 experiments uncovered that the proliferation of KIRC cells were obviously suppressed by silencing SLC25A5-AS1 expression. Nonetheless, the growth of KIRC cells were notably promoted by depressing LINC00472 and TCL6 expression ( Figure 8B).

DISCUSSION
In the present study, we performed a comprehensive bioinformatics analysis using RNA-seq data from TCGA database to identify key RNA transcript signatures associated with the clinical survival of KIRC patients. Totally, 1527 DElncRNAs, 54 DEmiRNAs, and 2321 DEmRNAs were screened between KIRC group and control group. Subsequently, a lncRNA-mediated ceRNA network was built among 81 DElncRNAs, 9 DEmiRNAs and 197 DEmRNAs. Moreover, 35 lncRNAs, 84 mRNAs, and two miRNAs (up-regulated miRNA-21 and miRNA-155) were significantly associated with OS of patients with KIRC. Moreover, there were the strongest negative interactions between miRNA-21/miRNA-155 and three prognosis-related lncRNAs (LINC00472, SLC25A5.AS1, and TCL6). Accordingly, the Cox regression analysis also revealed that SLC25A5-AS1 and miRNA-21 were dramatically linked with KIRC risk. SLC25A5-AS1 is an anti-sense lncRNA and located in X chromosome q24. Our result suggested SLC25A5-AS1 was  Yang et al. (2018) previously also constructed a ceRNA network based on a TCGA dataset and evaluated prognostic values of significant lncRNAs. They emphasized that SLC25A5-AS1 was under-expressed in KIRC tissues in comparison to adjacent normal tissues and associated with OS of KIRC patients. Notably, their finding showed an increased expression SLC25A5-AS1 was related to a favorable prognosis. However, our survival analysis indicated that a lower expression of SLC25A5-AS1 exhibited a better clinical outcome for patients suffering from KIRC. Possible explanations were the differences from sample size and the duration of follow-up. Interestingly, we also noted that SLC25A5-AS1 expression was positively correlated with pathogenic stage. There was the highest level of SLC25A5-AS1 in KIRC patients under pathogenic IV. Additionally, Li et al. (2019) pointed out SLC25A5-AS1 sponging miR-19a-3p restrained gastric cancer cell growth but increased apoptosis through PTEN/PI3K/AKT signaling pathway. We noted that SLC25A5-AS1 negatively interacted with miRNA-21. Differential expression analysis suggested that there was a higher level of miRNA-21 in KIRC samples than normal controls. Moreover, elevated expression of miRNA-21 was correlated with good survival outcomes. A growing body of research has demonstrated that miRNA-21 acted as a key predictor for KIRC. For example, Szabó et al. (2016) revealed that miRNA-21 was over-expressed in KIRC and strongly linked with metastasis and survival of KIRC patients. Later on, Xie et al. (2018) argued that a four-miRNA (miRNA-21-5p, miRNA-9-5p, miR-149-5p, and miRNA-30b-5p) axis could predict the clinical outcomes of patients undergoing KIRC. These studies provided direct evidence for our conclusion that miRNA-21 was a promising prognostic factor for KIRC. We also found that miRNA-21 closely interacted with seven mRNAs (ADCY2, ABAT, FGF1, FASLG, TGFBI, ALX1, and SLC30A10) in ceRNA regulatory network. Accumulating evidence demonstrated that FASLG, a member of the tumor necrosis factor superfamily, participated in the pathogenetic mechanisms of KIRC (Elsässer-Beile et al., 2003;Donskov et al., 2004;Xu T.R. et al., 2019). Sejima et al. (2012) highlighted that there was a remarkable worse survival for renal cell carcinomas (RCC) patients with FASLG mRNA-positive expression after radical nephrectomy in comparison to those with FASLG mRNA-negative expression, which was similar to our finding that FASLG expression was increased in KIRC patients and decreased level of FASLG exhibited good clinical prognosis. ABAT is central to the catabolism of gamma-aminobutyric acid, an inhibitory neurotransmitter in the mammalian central nervous system (Owens and Kriegstein, 2002). Chen et al. (2018) firstly discovered ABAT was dysregulated in patients with RCC compared with healthy controls and involved in the propanoate metabolism pathway by analyzing two microarray datasets (GSE53757 and GSE40435). Herein, our analysis showed that ABAT was down-regulated and its under-expression of ABAT was related to favorable prognosis of KIRC patients. TGFBI is located on chromosome 5q31.1, which played pivotal roles in pathologic responses of various cancers including RCC (Matsuda et al., 2008;Boguslawska et al., 2016;Pan et al., 2018). We found that higher TGFBI expression showed a better clinical prognosis for patients with KIRC. Moreover, elevated levels of FGF1, ALX1, and SLC30A10 were correlated with oncogenic outcomes while decreased level of ADCY2 was linked with good outcomes for KIRC patients. However, few groups considered the potential associations between these gene makers and survival prognosis of patients suffering from KIRC patients.
Our results also revealed that other two key lncRNAs (upregulated LINC00472 and TCL6) may be underlying prognostic makers for the KIRC patient survival. The KIRC patients with lower expression of LINC00472 and TCL6 had better prognosis. Similarly, Wang et al. also argued that LINC00472 and TCL6 acted as independent prognostic factors for predicting clinical outcomes of KIRC patients (Wang et al., 2018). Interestingly, these two lncRNAs were strongly associated with miRNA-21 and miRNA-155. Previous studies suggested that miRNA-155 was FIGURE 6 | The violin plot of the relationships between SLC25A5.AS1 expression and pathological stages. The correlations between four pathological stages and the expression level of SLC25A5.AS1 were evaluated. The horizontal axis presents the different pathological stage while the vertical axis presents the expression level of SLC25A5.AS1.
implicated with pathogenesis of renal carcinogenesis (Klimczak et al., 2017;Shiomi et al., 2019). Chen et al. conducted an integrative bioinformatics analysis based on a Gene Expression Omnibus dataset and indicated that miRNA-15 was involved in initiation and development of KIRC (Chen et al., 2016). Our analysis showed that the over-expression of miRNA-15 was significantly linked with good prognosis. More remarkably, there was close relationships between miRNA-15 and 12 mRNAs (STXBP5L, SCG2, SPI1, C12orf40, TYRP1, CTHRC1, TDO2, PTPRQ, TRPM8, ERMP1, CD36, and ST9SIA4). Interestingly, increased expressions of STXBP5L, SCG2, SPI1, C12orf40, TYRP1, CTHRC1, TDO2, PTPRQ, and TRPM8 were correlated with good survival outcomes but elevated expressions of ERMP1, CD36, and ST9SIA4 were associated with poor prognosis. Xu et al. previously reported that CD36 mRNA expression was elevated and the high expression of CD36 represented poor prognosis in KIRC patients, which was consistent with our finding (Xu W.H. et al., 2019). However, the relevant investigations regarding the relationships between survival outcomes of KIRC patients and other gene prognostic makers are lacking. Therefore, additional analyses still need to be performed to verify our conclusion.
limitations. Firstly, a comprehensive bioinformatics analysis with a larger sample size is still required to carried out to validate our results. Secondly, the relevant experimental evidence was also provided for our findings. Thirdly, the more detailed clinical information needed to be integrated into a large-scale prognosis analysis in the following research.

CONCLUSION
In summary, we identified 1527 DElncRNAs, 54 DEmiRNAs, and 2321 DEmRNAs between KIRC group and control group. Furthermore, a lncRNA-mediated ceRNA network was constructed. Besides, two miRNAs (miRNA-21 and miRNA-155) and corresponding three negatively correlated lncRNAs (LINC00472, SLC25A5.AS1, and TCL6) were associated with OS of patients with KIRC. Finally, seven mRNA targets of miRNA-21 and 12 mRNAs targeted by miRNA-155 were also pivotal prognostic makers for KIRC patients. Collectively, this study revealed several novel prognosis-related RNA makers for KIRC based on a ceRNA regulatory network, which provided deeper insights into developing promising therapeutic strategies for KIRC treatment.

AUTHOR CONTRIBUTIONS
QY and WC designed the study strategy and prepared the draft of the manuscript. WY and YC acquired the data and interpreted the data. JC and JY analyzed the data and prepared the figures and tables. XP and CC managed the data and collected the references. SG and XC revised the draft of the manuscript. All authors have agreed to the final content.