A novel gene signature related to oxidative stress predicts the prognosis in clear cell renal cell carcinoma

Data acquisition

All data about ccRCC were acquired from TCGA and ICGC. The TCGA KIRC cohort including 530 ccRCC samples and 72 normal samples was selected as a training set and was downloaded from the University of California, Santa Cruz (UCSC) Cancer Genomics Browser (https://genome-cancer.ucsc.edu/). Principal component analysis (PCA) was utilized to compare normal and cancer transcriptomic profiles in the TCGA KIRC cohort (Fig. S1). Additionally, gene expression and clinical data of 91 ccRCC samples in the validation set were acquired from the ICGC database (https://dcc.icgc.org/projects/RECA-EU).

Clustering of ccRCC samples

The consistent cluster were conducted on the ccRCC samples from the TCGA database via the R package ConsensusClusterPlus. After that, all samples were divided into seven clusters form k = 2 to k = 8 through Pearson correlation coefficient and Hierarchical clustering algorithm. Results were obtained from 100 repeated sampling on 80% of the sample. Finally, we determined the optimal cluster number (k = 3) according to the consistent cumulative distribution function (CDF) and the delta area diagram, which met the criteria of high consistency and regional stability under CDF curve.

Statistical analysis

The differentially expressed oxidative stress genes was analyzed by the R package DEseq2. The Kaplan–Meier method and log-rank test were used to generate survival curves. Cox regression models were constructed to calculate the hazard ratio (HR). Subsequently the following formula were utilized to calculate the risk score of each sample for further prediction of overall survival (OS): risk score = (exp mRNA 1 × coef mRNA1) + (exp mRNA 2 × coef mRNA2) + ... + (exp mRNAN × coef mRNAN). Moreover, the receiver operating characteristic (ROC) analysis was applied to access the efficiency of the prognostic signature and was validated as well as in the dataset ICGC. Prognostic signature and some clinical features were combined to construct a nomogram to further evaluate the survival probability of ccRCC patients.

All statistical tests were two-side, and P <0.05 was regarded as a statistically significant difference.

Tumor mutational burden (TMB)

Mutation data of ccRCC in the TCGA database were downloaded from UCSC Xena. The maf files and count the number of variants in ccRCC patients were identified through the R package “Maftools”. We then separated the patients into mutation and wild groups according to the mutation of gene, and compared the difference in prognosis.

Evaluation of tumor microenvironment

Here, one of the most commonly used tools for analyzing immune infiltration CIBERSORT was applied to evaluated the ccRCC microenvironment. Briefly, we calculated the respective proportions of list of 22 immune cells combined with the gene expression matrix in the different risk group patients.

Sensitivity analysis of chemotherapy

In this study, the NCI-60 cell line was utilized as a cancer cell sample group for anti-cancer drug testing. RNA-seq of 60 different cancer cell lines and drug activity data were obtained from the CellMiner (Reinhold et al., 2012). This was followed by the Pearson correlation analysis on the selected gene expression and sensitivity of 792 drugs.

Cell lines

The 786-O, OS-RC-2 and HEK-293T cell lines were obtained from the Shanghai Cell Bank Type Culture Collection Committee. Cells were cultured as previously described (Sun et al., 2021).

Plasmid and stable transfected cells construction

Plasmid, lentivirus and stable transfected cells were constructed as previously described (Sun et al., 2021). Specially, target DNA sequences were inserted into pLKO.1 plasmid.

RNA extraction, reverse transcribed PCR,quantitative real- time PCR (qRT-PCR)

Total RNA was extracted with TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). mRNA was reverse transcribed into cDNA via the Prime-Script™ RT Reagent Kit (TAKARA, Beijing, China). Quantitative PCR was performed with 2 × ChamQ Universal SYBR qPCR Master Mix* (Vazyme, Nanjing, China). The target gene UCN expression was normalized to the ratio of GAPDH. The specific primers (5′–3′) used for qPCR were as follows. (GAPDH-F: ACAACTTTGGTATCGTGGAAGG, GAPDH-R: GCCATCACGCCACAGTTTC; UCN-F: CAACCCTTCTCTGTCCATTGAC, UCN-R: CGAGTCGAATATGATGCGGTTC)

Cell proliferation assays

Cell proliferation assay was carried out via the CCK-8 kit according to the manufacturer’s instruction. Cells were seeded in a 96-well flat-bottomed plates, and each well contained 1,500 of 786-O cells or 2,000 of OS-RC-2 cells in 100 µL of cell suspension. After 24, 48, 72, and 96 h in culture at 37 °C, cell viability was measured through CCK-8 assays (Dojindo, Kumamoto, Japan). Each experiment consisted of five replicates and was repeated at least three times.

Additionally, cell proliferative ability was as well as measured by EdU incorporation assay (RiboBio, Guangzhou, China). The proliferative nuclei were stained with red fluorescence, while all nucleus were blue fluorescent light.

Colony-formation assays

786-O (1,000 cells/well) cells and OS-RC-2 (1,500 cells/well) cells were seeded into 6-well plate, After the cultivation for two weeks, the colonies were washed with cold 10% PBS twice, fixed with methano, and stained with crystal violet for 1 h at room temperature, followed by washing with water. Finally, the colony number in each well was counted and analyzed.

Transwell migration assay

For transwell migration assays, 6.0 × 10⁴ 786-O cells or 7.0 ×10⁴ OS-RC-2 cells were used for each well. The specific method was as previously described (Sun et al., 2021).

Intracellular ROS

The dichlorodihydrofluorescein diacetate (DCFH-DA) (Beyotime Institute of Biotechnology, China) was utilized to test intracellular ROS. The final concentration of DCFH-DA was adjusted to 10 mol/L via serum-free medium with 1:1,000. ccRCC cells transfected with sh-UCN and shControl plasmids were inoculated into 96-well plates (10,000 cells per well). Cells were washed with serum-free medium for three times after incubation at 37 °C for 30 min. Finally, the excitation fluorescence values were measured at 488 nm, while the emission fluorescence values were measured at 525 nm through a microplate reader (Thermo Fisher Scientific, Waltham, MA, USA).

Intracellular superoxide levels

The dihydroethidium (DHE; Beyotime Institute of Biotechnology, Jiangsu, China) was utilized to detect Intracellular superoxide levels. ccRCC cells transfected with sh-UCN and shControl plasmids were inoculated into 96-well plates (10,000 cells per well). Cells were washed with serum-free medium for three times after incubation at 37 °C for 30 min. Finally, the excitation fluorescence values were measured at 300 nm, while the emission fluorescence values were measured at 610 nm.

Gene mutations are identified in oxidative stress related genes in ccRCC

A total of 1,399 oxidative stress related genes were collected from the GeneCards and matched with the mRNA matrix of 530 ccRCC samples from the TCGA database. Then the oxidative stress related genes were analyzed by TMB. Some vital gene mutations with the ability to alter the corresponding protein sequence were considered, including missense-mutation, frame-shift-variant, splice-site–variant, nonsense-mutation, inframe-variant, translation-start-site-variant and nonstop-mutation. The top 50 oxidative stress genes with most mutation frequency was shown in Fig. 1A. Furthermore, we divided ccRCC patients into mutated and wild groups and compared the effects of gene mutations on the prognosis of patients. The results revealed that the gene mutation in oxidative stress gene set can reduced patients’ OS time. OS curves of two representative genes (MALAT1, RYR3) whose P values were the most significant are presented in Fig. 1B. It suggested that oxidative stress related genes were closely associated with genetic mutations in ccRCC.

Consistent clustering of ccRCC by oxidative stress genes

In order to research the possible effect of oxidative stress gene set on prognosis and clinical characteristics of ccRCC patients, we conducted an unsupervised clustering analysis and determined the optimal number of clusters by CDF. Considering both CDF and delta area, we noticed that CDF descending slope tended to be stable and had the optimal consistency and clustering confidence when K = 3 (Figs. S2A, S2B). In addition, the item-consensus graph also manifested that the sample classification region was stable enough when the cluster number was 3 (Fig. S2C). The matrix heatmap visually exhibited the distribution of the cluster samples (Fig. S2D). Ultimately, the heatmap of 1,399 oxidative stress related genes in 3 clusters was presented in Fig. 2A.

Subsequently, we utilized Kaplan–Meier curve analysis to study the effects of three clusters on survival time and the result indicated that there were significant differences in OS among the three clusters and the cluster 3 had the worst prognosis while the cluster 2 had the best prognosis compared to other clusters (Fig. 2B). Furthermore, we compared the clinical features among the three clusters including stage, grade, pathologic_TNM stage (Figs. 2C, 2D, Figs. S3A–S3C). In a brief, the clinical features of patients in different clusters varied visibly.

Oxidative stress genes are differently expressed in ccRCC

A total of 531 tumor tissues and 72 adjacent normal tissues in TCGA KIRC cohort were selected to study the oxidative expression profile of stress genes and explore the aberrant oxidative stress genes in ccRCC. 139 differentially expressed oxidative stress genes were extracted after the first round of screening, in which 82 oxidative stress genes were upregulated while 57 oxidative stress genes were downregulated in tumor tissues compared with adjacent normal tissues, and the screening criteria was P value <0.05 and fold change >1.0 (Table S2).

Secondly, 62 prognostic genes were identified through Univariate Cox analysis (P-value <0.05) among 139 differentially expressed oxidative stress in ccRCC, and the specific expression of these prognostic genes were shown in Figs. 3A, 3B. Moreover, the prognostic genes were found enriched in many pathways related to oxidative stress such as response to oxidative stress, regulation of cellular response to oxidative stress through GO and KEGG analysis (Figs. 3C, 3D).

Construction of the oxidative stress gene signature for predicting survival in ccRCC

As is shown in Fig. 4A, least absolute shrinkage and selection operator (LASSO) regression (Figs. S4A, S4B) and multivariate Cox analysis (Table 1) were applied to determine statistically significant oxidative stress genes which were related to survival time. Finally, four genes (UCN, PLG, FOXM1, HRH2) with prognostic significance were filtered to construct a prognostic signature after four rounds of screening. Subsequently, a prognostic model was established for the following analysis: risk score = (0.30010645 × UCN expression) - (0.16955462 ×PLG expression) + (0.23343259 ×FOXM1 expression) - (0.12178059 × HRH2 expression). To assess the prognostic value of four selected genes, the ccRCC patients in TCGA KIRC cohort were divided into high-risk group and low-risk group according to the median risk score. In addition, the Kaplan–Meier curve demonstrated that the high-risk group patients tended to have worse OS than the low-risk group patients (Fig. 4B). Additionally, area under the receiver operating characteristic (AUROC) curves for the four oxidative stress gene prognostic models were plotted, and the area under the time-dependent ROC curve for 1-, 3- and 5-year OS were 0.77, 0.70, and 0.71 severally, indicating that this model has great predictive value (Fig. 4D). The heatmap of four oxidative stress genes between high and low risk groups further visually evidenced that a higher risk score predicted a worse prognosis (Fig. 4F). Meanwhile, in order to more exactly assess the survival probability, we constructed the prognostic nomogram for OS at 1, 3, and 5 years (Fig. 4H).

Furthermore, we verified the predictive feasibility of the prognostic signature by 91 ccRCC samples from the ICGC database. All these samples were divided into high-risk and low-risk group by the above formula. As expected, the prognosis was significantly worse in the high-risk group than in the low-risk group (Fig. 4C). The AUCs of the four-gene signature for the 1-, 3- and 5-year OS were 0.677, 0.673 and 0.681 (Fig. 4E). Similarly, we plotted the heatmap of four genes in the ICGC and constructed the prognostic nomogram for OS at 1, 3, and 5 years (Figs. 4G, 4I).

We explored the relationship between oxidative stress gene signature and some clinical characteristics of ccRCC patients. The results revealed that men tended to have higher risk scores compared with women (Fig. S5A). Meanwhile, higher risk scores were closely associated with poorer disease-specific survival (DSS) outcomes and lymph node metastasis (Figs. S5B, S5C).

Immune infiltration in high risk and low risk group

Oxidative stress had been reported to be related to immune infiltration (McGarry et al., 2018; Mendiola et al., 2020). Immune infiltrates in the tumor microenvironment had been verified to exert important influence in tumor development and had the ability to affect the clinical outcomes of variety types of cancer patients (Zhang & Zhang, 2020). Here, we evaluated the immune infiltration status of 22 different immune cell types in the high risk and low risk ccRCC samples via the method of CIBERSORT. It was noteworthy that the high-risk group was more likely linked with higher percentage of CD8+T cells, Tregs, macrophages, T cells follicular helper and plasma cells. In contrast, high-risk group presented significantly higher proportions of Mast cells resting, B cells naïve and Dendritic cells resting (Figs. 5A, 5B). Generally, oxidative stress genes were closely related to immune infiltration. The difference in prognosis among ccRCC patients might be attributed to the different distribution of immune cells.

UCN independently predict poor prognosis and may be involved in oxidative stress

Subsequently, we chose the most significant P-value and highest hazard ratio gene UCN among the panel of four oxidative stress genes for next study. In TCGA KIRC cohort, patients with high UCN expression had poorer OS, DSS and progression-free survival (PFS) than low expressed patients (Figs. 6A–6C). We found that UCN was overexpressed in tumor samples compared to normal samples in both TCGA and ICGC cohorts (Figs. 6D, 6E). Then, the nomograms of UCN for OS were constructed to evaluate the survival probability of ccRCC patient (Fig. 6F). Furthermore, we performed GSEA analysis to look for the cancer hallmark pathways of UCN and the top enrichment pathways were shown in Fig. 6G. The result suggested that UCN exert non-negligible effect on oxidative phosphorylation. Generally, the level of oxidative stress was positively correlated with the content of ROS and superoxide in vivo. In order to investigate whether UCN impacted oxidative stress in tumor cells actually, we constructed a stable knocking down UCN ccRCC cell lines (Fig. 7A) and employed a series of specific kits to determine the ROS through DCFH-DA and the superoxide content by dihydroethidium. The result of the ROS and superoxide content further demonstrated that the UCN was of great importance in oxidative stress (Figs. 6H, 6I).

UCN promote proliferation and migration of ccRCC

Next, we implemented some functional experiments to determine the role of UCN in ccRCC cells. Comparing with the control cells, we observed that colony formation was markedly reduced in silencing UCN cells (Fig. 7B). Similarly, we also observed that the proliferation of silencing UCN cells was significantly weakened (Figs. 7C, 7D). Furthermore, to explore whether silencing UCN could as well as affect the migration ability of tumor cells, we performed the transwell migration assay and proved that knocking down UCN could vastly attenuated tumor migration (Fig. 7E). In brief, UCN exerted important influence in the development and migration of cancer cells.

Finally, the data of selected genes in the NCI-60 cell line was acquired from the CellMiner database, and the relationship between the UCN expression and drug sensitivity was analyzed. It was worth noting that, the sensitivity of many drugs was associated with the UCN (P < 0.001) (Fig. 8A). For instance, the high expression of UCN indicated the high IC50 levels of Gemcitabine, Elliptinium Acetate, Clofarabine, etc. (Fig. 8B). Hopefully, these drugs were promising as approach for tumor-targeted therapy in the future.