Identification of novel prognostic biomarkers in renal cell carcinoma

Identification of DEGs, DEMs and prognostic biomarkers in RCC

We identified prognostic biomarkers by a standardized procedure (Figure 1A). The data in the GSE105288 dataset showed that there were a total of 480 DEGs in the RCC samples, including 164 upregulated DEGs and 316 downregulated DEGs (Figure 1B), while there were only 5 upregulated DEMs and 2 downregulated DEMs in the GSE116251 dataset (Figure 1C). Moreover, we obtained 1625 upregulated DEGs and 1321 downregulated DEGs in RCC tissues by using GEPIA2 database (Figure 1D and 1E).

Figure 1. DEGs/DEMs and prognosis-related genes/miRNAs identification in different datasets and databases. (A) A standardized procedure for the identification of prognostic biomarkers; (B) Volcano plot of DEGs from GSE105288; (C) Volcano plot of DEMs from GSE116251; (D) Intersection of the upregulated genes and (E) downregulated genes from the GSE105288 dataset and the GEPIA2 database; (F) Intersection of 49 prognostic miRNAs from the Oncolnc database and 7 DEMs from the GSE116251 dataset.

The top 500 genes (Supplementary Table 1) that were related to overall survival (OS) and disease-free survival (DFS) separately in GEPIA2 were then screened. Among them were 179 genes related to both OS and DFS (Figure 1D and 1E).

The overlapped genes in both GSE105288 and the above 179 genes in GEPIA2 included 12 ones as follows: WDR72, ALDH6A1, CDS1, SLC25A4, MTURN, ERBB2, NR3C2, OGDHL, BSPRY, HADH, DNASE1L3 and CLDN10 (Figure 1E).

Furthermore, as for the miRNAs OS profile a total of 448 miRNAs were downloaded from the Oncolnc platform. There were 49 miRNAs significantly related to OS (Figure 1F). Hsa-miR-21-5p was the only overlapped miRNA in both DEMs and among the above 49 miRNAs (Figure 1F).

Taken together, the 12 genes and hsa-miR-21-5p were selected as prognostic biomarkers in RCC.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment

To understand the underlying functions of the significant DEGs, Metascape was used for functional enrichment. The outcomes of GO analysis suggested that “extracellular structure organization”, “endoplasmic reticulum lumen”, and “response to oxygen levels” were the most important items for the upregulated significant DEGs (Figure 2A), while for the downregulated significant DEGs they were “apical plasma membrane”, “kidney development”, and “monovalent inorganic cation homeostasis” (Figure 2C). Moreover, in DEGs analyzed in the KEGG pathway analysis, the DEGs were involved in multiple pathways, including “Staphylococcus aureus infection”, “PPAR signaling pathway”, “Phagosome”, “Aldosterone-regulated sodium reabsorption”, and the “Collecting duct acid secretion” pathways (Figure 2B and 2D). These pathways may provide researchers with directions of further mechanistic investigations of these DEGs in RCC.

Figure 2. Results of DEGs for the GO and KEGG pathway analysis. (A) The enriched GO biological processes and (B) KEGG pathways of the significant upregulated DEGs; (C) The enriched GO biological processes and (D) the KEGG pathways of significant downregulated DEGs.

Hsa-miR-21-5p expression was negatively related to 11 of the coding genes in RCC

To further study the association between hsa-miR-21-5p and the 12 coding genes in RCC, we compared their levels between tumor and normal samples from UCSC Xena as well as the TCGA database. Consistent with the results of the GSE105288 and GSE116251 datasets, the levels of all 12 coding genes were lower in the RCC samples (Figure 3A), while the expression of hsa-miR-21-5p was significantly increased (Figure 3B). To figure out the expression relationship between hsa-miR-21-5p and the 12 coding genes, we next checked the coexpression of hsa-miR-21-5p and the 12 genes in StarBase. The results indicated that hsa-miR-21-5p was negatively related to these coding genes except for CLDN10 (Figure 3C), which was the only gene with a P value over 0.05.

Figure 3. Expressionss of the prognostic biomarkers and their coexpression in RCC. (A) The 12 coding genes were downregulated and (B) hsa-miR-21-5p was upregulated in tumor samples compared with normal samples; (C) Showing the coexpression of prognostic coding genes and hsa-miR-21-5p in RCC, among which 11 coding genes (WDR72, ALDH6A1, CDS1, SLC25A4, MTURN, ERBB2, NR3C2, OGDHL, BSPRY, HADH, DNASE1L3) were negatively related to the expression of hsa-miR-21-5p. CLDN10 was the only gene with a P value over 0.05.

Decreased level of the 12 coding genes and increased level of hsa-miR-21-5p showed poorer prognosis of RCC

To determine the exact prognostic significance of the identified genes, we evaluated the OS and DFS of the 12 coding genes as well as the OS of hsa-miR-21-5p in GEPIA2 and Oncolnc, respectively. The results showed that high expression of hsa-miR-21-5p was associated with poor OS and DFS (Supplementary Figure 1), and a decreased level of all 12 coding genes indicated a poor OS and DFS (Figure 4A and 4B), which suggested that these genes and hsa-miR-21-5p could indeed act as biomarkers in the prognostication of RCC.

Figure 4. Kaplan-Meier curves of OS and DFS of 12 prognostic coding genes. Lower expression of all 12 coding genes was relevant to both unfavorable OS (A) and worse DFS (B) in patients with RCC.

The prognostic biomarkers were related to the clinical stages of RCC

To investigate the association between the prognostic biomarkers and the clinical stages, the levels of the biomarkers at different stages were evaluated. As shown in Figure 5, the levels of the 12 coding genes were lower in higher stages, and some genes, for instance, WDR72, ALDH6A1, OGDHL, and DNASE1L3, were gradually decreased as stage increased. However, increased expression of hsa-miR-21-5p was positively related to increased stage.

Figure 5. The association of the prognostic biomarkers and clinical stages. Lower levels of the 12 coding genes were observed in higher clinical stages. However, the level of hsa-miR-21-5p was elevated in higher clinical stages.

Gene signatures predicted the OS and DFS of RCC

Multivariate Cox regression analysis was used to further explore the signature of the multiple genes regarding the survival of RCC. A group of seven genes (WDR72, ALDH6A1, CDS1, HADH, DNASE1L3, CLDN10, hsa-miR-21-5p) were finally exported into the model as a molecular signature to predict OS (Supplementary Table 2). Five genes (WDR72, ALDH6A1, OGDHL, HADH, and DNASE1L3) exported were found to be a molecular signature to predict DFS (Supplementary Table 3).

Based on the levels of the above seven biomarkers in the risk model, all patients were divided into two groups: low- and high-risk groups (Figures 6A, 7A). The patients in the low-risk group had a more favorable OS and DFS than those in the high-risk group (Figure 6B, 7B). The accuracy of this OS model (the area under the ROC curve) was 76.6%, 71.6% and 72.9% at 1, 3 and 5 years, respectively (P<0.001) (Figure 6C). Similarly, the accuracy of this DFS model (the area under the ROC curve) was 74.3%, 74.7% and 77.4% at 1, 3 and 5 years, respectively (P<0.001) (Figure 7C). These data suggested that the signature was quite valuable in predicting OS and DFS of RCC, and it was more convincing in predicting OS than DFS.

Figure 6. OS of the signature with multiple genes in RCC. (A) Patients were classified by risk score and their survival status; (B) OS of the signature of multiple genes in this model, which indicated that patients in the low-risk group had a more favorable OS than those in the high-risk group. (C) ROC curves suggested that the accuracy of this model was 76.6%, 71.6% and 72.9% at 1, 3 and 5 years, respectively.

Figure 7. DFS of the signature with multiple genes in RCC. (A) Patients were classified by risk score and their survival status; (B) DFS of the signature of multiple genes in this model, which indicated that patients in the low-risk group had a better DFS than those in the high-risk group. (C) ROC curves suggested that the accuracy of this model was 74.3%, 74.7% and 77.4% at 1, 3 and 5 years, respectively.

Increased WDR72 expression inhibited the survival and invasion of RCC cell lines in vitro

WDR72 was found to be the most significant prognostic gene among the identified coding genes in RCC. To further investigate its role in cell proliferation and invasion, WDR72 was successfully overexpressed in RCC 769P (clone #1 and #2) and 786-O (clone #1 and #2) cells by a lentivirus system (Figure 8A). CCK-8 assay displayed overexpression of WDR72 significantly inhibited cell proliferation (Figure 8B). Then, by cell invasion assay, we found that overexpression of WDR72 decreased the invasiveness of RCC cells (Figure 8C). These data suggested that increased expression of WDR72 inhibited the survival and invasion of RCC cell lines. Therefore, WDR72 might serve as a cancer suppressor in RCC.

Figure 8. Role of WDR72 overexpression in cell survival and invasion of RCC. (A) The results of western blot showed that WDR72 was successfully overexpressed in 769-P and 786-O cells; overexpression of WDR72 remarkably decreased the survival (B) and invasiveness (C) of 769-P and 786-O cells.

Differential gene expression analysis

The gene expression profile of GSE105288 [44] and the miRNA expression profile of GSE116251 [45] were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo). The GSE105288 dataset includes 35 RCC samples and 9 normal renal samples. GSE116251 is comprised of 18 RCC tissues and 18 normal renal tissues. The comparison was performed in GEO2R to identify DEGs and differentially expressed miRNAs (DEMs). Since the adj. P-value could amend the false positive rate, we selected the adj. P-value< 0.05 and |log₂FC|>1 as cut-off criteria. In addition, with q-value < 0.01 and |log₂FC|>1, the DEGs of RCC in the Gene Expression Profiling Interactive Analysis (GEPIA2) database were also downloaded using the ANOVA method. The DEGs from the GSE105288 dataset and the GEPIA2 database were considered as the significant DEGs, while the DEMs from GSE116251 were defined as significant DEMs. GEPIA2 is an open web version for visual analysis based on the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from the Cancer Genome Atlas (TCGA) and the GTEx projects [46]. GEO2R (https://https://www.ncbi.nlm.nih.gov/geo/geo2r), is an interactive web tool to identify DEGs through comparing grouped samples based on the GEO query and R-packages.

Functional and pathway enrichment analysis

Metascape is an efficient platform for gene function annotation analysis [47]. It is composed of functional enrichment, interactome investigation, genomic annotation and packages approximately 50 independent knowledge-bases into a single platform. Significant DEGs were input into the Metascape program to carry out GO and KEGG pathway analysis. Homo sapiens was selected and the P-value <0.05 was considered statistically significant.

Prognosis-related genes and miRNAs screening

Both the most significant OS-related genes and DFS-related genes in RCC were screened in the GEPIA2 database and the top 500 genes with P-value<0.05 were downloaded. Oncolnc (http://www.oncolnc.org) is an online program containing survival data for 8647 patients of 21 cancer types from the TCGA database. The clinical and expression data of the mRNAs, microRNAs, and long non-coding RNAs can be directly downloaded from Oncolnc. Thus, we downloaded the miRNAs OS profile of RCC and the BH-adjusted p-value <0.05 was thought to indicate prognosis-related miRNAs.

Prognostic biomarkers identification

The significant DEGs and significant prognosis-related overlapped genes were used to construct prognostic biomarkers in RCC. Similarly, the significant DEMs and prognosis-related miRNAs were also used to establish prognostic biomarkers in RCC. OS was defined as the time from treatment of RCC to death. Death could be either cancer-related or cancer-unrelated. DFS was defined as time from treatment of RCC to recurrence. A recurrence was considered as the development of RCC confirmed by pathological or radiological diagnosis at the operative site, in the regional lymph nodes, or a distant metastasis.

Verification of expression and prognostic significance of the prognostic biomarkers

UCSC Xena (https://xena.ucsc.edu/) is a visual exploration resource for both data downloading and online analysis based on Xena Browser. It also allows biologists to figure out the relationships between genomic and clinical data. StarBase (http://starbase.sysu.edu.cn/index.php) is a tool for exploring noncoding RNA functions and regulatory mechanisms. StarBase identifies millions of interactions among miRNA, ncRNA and mRNA using multiple dimensional data. The mRNA expression data of prognostic genes and prognostic miRNAs were downloaded from the UCSC Xena database and the TCGA database (https://www.cancer.gov/tcga), respectively. The prognostic significance of prognostic genes and miRNA were evaluated in the GEPIA2 and Oncolnc platforms separately. Furthermore, we analyzed the coexpression of the prognostic genes and prognostic miRNA in RCC with the StarBase version 3 database.

Association between prognostic biomarkers and clinical stage

We further analyzed the association between the expression of these prognostic biomarkers and the clinical stage of RCC patients. The expression levels of all the 13 biomarkers were compared in different stages (I, II, III, and IV). Differences between every two stages were analyzed by t-test. The ggplot2 package of R software (4.0.2) was used for result visualization. P value <0.05 was considered statistically significant.

Multivariate Cox regression analysis

To understand the combinatorial effect of the identified biomarkers, multivariate Cox regression analysis was performed and a risk model was established. Patients were divided into high-risk and low-risk groups by the cut-off value (1.0) of the risk score. Then, the survival of the two groups was evaluated by Kaplan-Meier survival curves with log-rank test. Time-dependent ROC curves were performed to assess the performance of the signature classifier at 1, 3 and 5 years.

Cell culture

The HEK293T, 786-O and 769P cell lines were purchased from the American Type Culture Collection (ATCC). The HEK293T, 786-O and 769-P cells were respectively cultured in RPM1-DMEM medium (with L-glutamine) and RPM1-1640 medium (with L-glutamine) (Fcmacs Biotech, Nanjing, China) in a 5% CO₂ incubator at 37° C. Both kinds of medium were complete medium supplemented with 10% fetal bovine serum (FBS) (Fcmacs Biotech) and 1% penicillin and streptomycin (Sangon, Shanghai, China).

Gene overexpression

Lenti-sgRNA^WDR72 was constructed by lenti sgRNA (MS2). The sgRNA^WDR72 plasmid was mixed with the REV, GAG and VSVG plasmids (Hanheng Biotech, Shanghai, China) to construct Lenti-sgRNA^WDR72. After transfection for 48 h, recombinant lentiviruses were obtained. Lenti-dCas9 and Lenti-MPH were constructed in the same way. The 786-O and 769-P cells were infected by the Lenti-dCas9 and Lenti-MPH to generate 769P^dCas9+MPH+ and 786-O^dCas9+MPH+ cells, which were finally infected with Lenti-sgRNA^WDR72. The 1640 medium containing puromycin was used to screen the cells for 3-4 days. The cells with stably overexpressed WDR72 were selected for further study. We also validated the expression of WDR72 by western blot.

Western blot

Radioimmunoprecipitation assay (RIPA) solution was utilized to extract the total protein of the cells and BCA reagent (Beyotime, Shanghai, China) was used to determine their concentration. Protein samples were separated by electrophoresis and then transferred to polyvinylidene difluoride membranes (Millipore, USA) for 2 h. Subsequently, the membranes were blocked with 3% bovine serum albumin (Sangon, Shanghai, China) and incubated with anti-WDR72 primary antibody (Cell Signaling, Danvers, MA, USA) at 4° C overnight. The next day, 1×TBST was used to wash the membranes 3 times. Then, a secondary antibody was used to incubate the membranes at room temperature for 1.5 h. Finally, the bands were scanned by a gel-imaging platform (UVP, CA, USA).

Cell proliferation assay

According to the instructions of the manufacturer, the Cell Counting Kit-8 (CCK-8) (Vazyme, Nanjing, China) assay was used to measure cell proliferation. Experimental cells were inoculated in the 96-well plate (each well with 1 × 10⁴ cells) and cultured in an incubator for 24 h, 48 h, and 72 h. Then, 10% CCK-8 reagent was added into each well and continuously cultured for 1 h. The optical density (OD) at 450 nm was measured. Three independent experiments were carried out.

Cell invasion assay

Invasion assay was conducted using 24-well Transwell chambers (Corning Incorporated, USA). 786-O^WDR72(+), 769-P^WDR72(+) cells as well as the control were seeded in the upper chamber (2×10⁴ cells) with 500 μL 1640 medium. Then, the cells were cultured in the incubator for 24 h. After wiping away the noninvading cells on the surface of the upper chamber, we stained the cells in the lower chamber with crystal violet (BBI life sciences, Shanghai, China) at room temperature for 30 min. The chamber was photographed before the cells were counted. Five independent areas of each well were acquired and the mean number of cells was used for graphing.

Statistical analysis

Every experiment was repeated no less than 3 times. The comparisons of different groups were performed with independent sample t-test. A two-sided P value<0.05 was defined as statistically significant. GraphPad prism 8.0 and R software (4.0.2) were used for graphing.