A novel oxidative stress-related genes signature associated with clinical prognosis and immunotherapy responses in clear cell renal cell carcinoma

Background Oxidative stress plays a significant role in the tumorigenesis and progression of tumors. We aimed to develop a prognostic signature using oxidative stress-related genes (ORGs) to predict clinical outcome and provide light on the immunotherapy responses of clear cell renal cell carcinoma (ccRCC). Methods The information of ccRCC patients were collected from the TCGA and the E-MTAB-1980 datasets. Univariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO) were conducted to screen out overall survival (OS)-related genes. Then, an ORGs risk signature was built by multivariate Cox regression analyses. The performance of the risk signature was evaluated with Kaplan-Meier (K-M) survival. The ssGSEA and CIBERSORT algorithms were performed to evaluate immune infiltration status. Finally, immunotherapy responses was analyzed based on expression of several immune checkpoints. Results A prognostic 9-gene signature with ABCB1, AGER, E2F1, FOXM1, HADH, ISG15, KCNMA1, PLG, and TEK. The patients in the high risk group had apparently poor survival (TCGA: p < 0.001; E-MTAB-1980: p < 0.001). The AUC of the signature was 0.81 at 1 year, 0.76 at 3 years, and 0.78 at 5 years in the TCGA, respectively, and was 0.8 at 1 year, 0.82 at 3 years, and 0.83 at 5 years in the E-MTAB-1980, respectively. Independent prognostic analysis proved the stable clinical prognostic value of the signature (TCGA cohort: HR = 1.188, 95% CI =1.142-1.236, p < 0.001; E-MTAB-1980 cohort: HR =1.877, 95% CI= 1.377-2.588, p < 0.001). Clinical features correlation analysis proved that patients in the high risk group were more likely to have a larger range of clinical tumor progression. The ssGSEA and CIBERSORT analysis indicated that immune infiltration status were significantly different between two risk groups. Finally, we found that patients in the high risk group tended to respond more actively to immunotherapy. Conclusion We developed a robust prognostic signature based on ORGs, which may contribute to predict survival and guide personalize immunotherapy of individuals with ccRCC.


Introduction
Renal cell carcinoma (RCC) is a common malignant tumor of genitourinary system with an incidence only secondary to prostate cancer and bladder cancer, affecting nearly 431,000 new patients and 179,000 related deaths in 2020 worldwide (1). Clear cell renal cell carcinoma (ccRCC) is the most prevalent histological types of RCC, which accounting for about 70-80% of RCC (2). Clinically, approximately 30% of patients would experience tumor metastasis after curative surgical resection during the follow-up (3). The average amount of time from curative surgical resection to metastasis has not been reported in detail in the literature. There are no metastasis markers for ccRCC currently known. Imaging is used primarily to determine whether metastasis has occurred. In recent years, although targeted therapy and immunotherapy have greatly improved clinical outcome of patients with advanced ccRCC and become mainstays of treatment for advanced ccRCC (4, 5), however, a significant number of patients did not respond to these treatments (6). Due to heterogeneous disease, ccRCC patients with similar clinical condition may have distinctive prognosis. In clinical practice, it is a great challenge to early identify stratification of risk in patients with ccRCC and provide accurate clinical individualized therapeutic. With the development of sequencing, mounting evidences indicated that prognostic signature based on combining genes expression and clinical feature could be used as a new biomarker to optimize risk stratification, predict clinical prognosis and evaluate response to clinical treatment (7)(8)(9). These novel prognostic signatures could help doctors implemented individualized therapy to extend the survival of patients. Thus, establishing a reliable prognostic model is crucial to predict clinical prognosis and guide personal precision therapy for patients with ccRCC.
Oxidative stress is a common pathological phenomenon in the body, which is characterized by the imbalance between synthesis of oxidants and antioxidants, leading to the accumulation of large amounts of reactive oxygen species (ROS) (10). It has been reported that high levels of ROS contribute to tumorigenesis and tumor progression through a variety of pathways, such as tumor signaling pathways, tumor microenvironments, immune escape, metastasis, DNA mutations, and angiogenesis (11)(12)(13)(14). Moreover, high levels of ROS also could affect the tumor development through influencing chemotherapeutic resistance and inducing cell apoptosis (15)(16)(17). Recent studies have shown that ORGs signature could be used as a biomarker for predicting clinical outcome and treatment responses in many cancers (18)(19)(20)(21). In urologic cancer, ORGs signature was developed for predicting clinical outcome and immune status in patients with bladder cancer (22). Nevertheless, the clinical value of prognostic signature based on ORGs for ccRCC are needed investigated in depth.
In this study, a prognostic signature was established based on 9 ORGs to predict clinical outcome in individuals with ccRCC, and the effectiveness and reliability of the prognostic signature were further confirmed. In addition, immune cell infiltration and immunotherapy responses were comprehensively investigated.

Materials and methods
The complete procedures of this study was presented in Figure 1.

Raw data collection
Transcripts and the corresponding clinical materials of ccRCC were acquired from the TCGA (https://genomecancer.ucsc.edu) and ArrayExpress datasets (E-MTAB-1980 dataset, https:// www.ebi.ac.uk/arrayexpress/). Patients with survival less than one month or missing clinical information were excluded for our study.

Identification of ORGs
ORGs were systematically searched from the GeneCards database (https://www.genecards.org/), with the screening threshold as relevance score ≥ 7. Finally, 1399 ORGs were included in our study (Supplementary Table 1).

Screening for differentially expressed ORGs
Differentially expressed genes (DEGs) between tumor and normal tissues were analyzed by the Package "limma", with the screening threshold as |log2FC| > 1 and adjusted p < 0.05.

Establishment and validation of ORGs prognostic signature
Firstly, univariate Cox analysis was conducted to screen prognosis-related ORGs in the TCGA database. Next, the least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox regression analysis were performed to screen genes for developing the prognostic signature for ccRCC. Subsequently, risk score was calculated with following formula: ÂExpression mRNA3 + … + b mRNAn Â Expression mRNAn : According to the medium of risk score, individuals were classified into two groups (high risk group Vs low risk group). OS of the two groups was compared by K-M analysis. Finally, Receiver operating characteristic (ROC) curve was generated to assess the predictive accuracy and sensitivity of the prognostic signature.

Establishment of a predictive nomogram
A nomogram was established containing risk score and independent prognostic factors based on results of the univariate and multivariate Cox regression analyses. The calibration curve was drawn to evaluate the predictive capability of nomogram by "rms" R package.

Stratified analysis and comprehensive analysis of the ORGs signature
Stratified analysis was carried out to assess clinical value of the ORGs prognostic signature based on clinical features. In addition, to better evaluate the role of the ORGs signature in the ccRCC development, the differences in risk score were compared in different subgroups based on clinical features. The flowchart of the present study design.
Wu et al. 10.3389/fonc.2023.1184841 Evaluation of tumor immune microenvironment The ssGSEA algorithm was performed to calculate the infiltrating levels of 29 immune-related functional indicators, including immune cells and immune-related pathways. Next, the infiltration proportion of immune cell types was quantified with CIBERSORT algorithm.

Evaluation of the response to immunotherapy
The expression level of four major immune checkpoints, including PD-1, PD-L1, CTLA4 and LAG3, were analyzed in ccRCC tissues. Moreover, the association between risk score and the expression level of immune checkpoints was examined by the Spearman method.

Statistical analysis
Statistical analyses were implemented by R software (version R-4.1.2) and GraphPad Prism (version 8.0.2). The difference in the continuous data between two groups was analyzed using Student's t test. The correlation analysis was implemented by the Spearman method. A p value < 0.05 indicating a statistically significant (*p < 0.05, **p < 0.01, and ***p < 0.001).

Identification of oxidative stress-related prognostic genes
Among 1399 ORGs, 242 DEGs were existed between tumor tissues and normal tissues ( Figure 2A). Then, 99 ORGs associated with prognosis were identified based on the univariate Cox regression analysis ( Figure 2B). Interactions of these 99 genes were further visualized with protein-protein interaction (PPI) network ( Figure 2C).

Establishment of the ORGs prognostic signature
We identified 20 prognostic ORGs by Lasso Cox regression analysis ( Figures 3A, B). Subsequently, multivariate Cox regression was performed to further filter out the candidate genes that were significantly related to survival, and finally identified 9 genes (ABCB1, AGER, E2F1, FOXM1, HADH, ISG15, KCNMA1, PLG, and TEK) for construction of prognostic model. The following equation was adopted to calculate risk score: risk score = ( Patients were assigned into high risk group and low risk group based on the consideration of the median risk score ( Figure 3C). The survival scatter chart revealed that individuals in the high risk group displayed a higher risk of mortality ( Figure 3D). The expression level of E2F1, FOXM1, AGER, and ISG15 were higher in the high risk group, while the expression level of HADH, KCNMA1, ABCB1, PLG and TEK were higher in the low risk group ( Figure 3E). Compared to normal tissues, the expression level of E2F1, FOXM1, AGER, ISG15, and KCNMA1 were significantly higher in ccRCC tissues, while the expression level of HADH, ABCB1, PLG and TEK were significantly lower (Supplementary Figure 1). We performed survival analysis on the 9 genes and found that patients with high expression level of identified genes, such as E2F1, FOXM1, AGER and ISG15, displayed a significantly poor outcome, while patients with high expression level of identified genes, such as KCNMA1, HADH, ABCB1, PLG and TEK, displayed a significantly favorable outcome (Supplementary Figure 2). K-M analysis shown that patients with low risk score displayed a significantly favorable outcome compared with patients with high risk score ( Figure 3F). The area under the ROC curve (AUC) of 1 year (AUC=0.81), 3 years (AUC=0.76), and 5 years (AUC=0.78) were all larger than 0.70 ( Figure 3G), suggesting that the prognostic signature displayed a favorable accuracy in predicting survival of patients with ccRCC.

Validation of the ORGs signature in the E-MTAB-1980 cohort
To further test the reliability of the ORGs prognostic signature, the same analyses were implemented in E-MTAB-1980 cohort. Individuals with high risk score had an adverse survival status ( Figures 4A, B). Expression level pattern of these genes that make up the ORGs prognostic signature were consistent with those in the TCGA ( Figure 4C). In addition, significantly shorter OS of patients with high risk score was observed ( Figure 4D). The AUC for 1 year, 3 years, 5 years was 0.8, 0.82 and 0.83, respectively ( Figure 4E), demonstrating that the ORGs prognostic signature displayed a good accuracy in predicting OS for patients with ccRCC.

Stratified analysis
To further confirm the accurately and independently prognostic value of the ORGs signature in ccRCC, stratification analysis was performed in different subgroups based on clinical features. We found that significantly poor OS was observed in the high risk group in all subgroups, such as age ( Figures 5A, B (Figures 5K, L). These findings suggested that the ORGs signature had universal applicability in predicting OS for patients with ccRCC.

Correlation of the ORGs signature with clinical features
To further explore the correlation between the ORGs signature and clinical features, risk score was calculated in the different subgroups based on clinical features. From the TCGA cohort analysis results, it was shown that risk score was significantly higher in G3-4, T3-4, M1, and stageIII-IV subgroups than those in the corresponding early clinicopathological stage subgroups ( Figures 6A, B). In addition, risk score of patients in the E-MTAB-1980 cohort was significantly higher in T3-4, N1-2, M1, and stageIII-IV subgroups than those in the corresponding early clinicopathological stage subgroups (Figures 6C, D). These results indicated that ORGs signature might act as a malignancy biomarker for ccRCC.

Independence of the ORGs signature and nomogram construction
To further confirm whether risk score could be used as an independent prediction biomarker for ccRCC survival, univariate and multivariate Cox regression analyses were implemented. Results demonstrated that risk score of the ORGs signature could be used as an independent indicator for predicting survival of ccRCC ( Figures 7A-D). Additionally, nomogram was widely used in predicting OS of patients with cancer based on nomogram scores, therefore, a nomogram was constructed by integrating independent prognostic markers based on multivariate Cox regression analysis ( Figure 7E). The calibration curve illustrated that the nomogram had an excellent capability in predicting survival at 1, 3, and 5 years ( Figures 7F-H).

Relationship between the ORGs signature and immune infiltration in the TCGA cohort
The immune infiltration landscape was explored in patients by multiple algorithms. The sGSEA algorithm results shown that infiltration abundance of most immune cells, including TIL, Th2_cells, Th1_cell, Tfh, T_helper_cells, pDCs, Macrophages, CD8+_T_cells, and aDCs, were considerably higher in individuals in the high risk group, whereas infiltration abundance of Mast_cells and iDCs were obviously higher in individuals in the low risk group (Figures 8A, B). Additionally, significantly higher immune function scores, such as T_cell_costimulation, T_cell_co-inhibition, Parainflammation, Inflammation-promoting, Cytolytic_activity, check-point, and APC_co_stimulation, were observed in individuals in the high risk group ( Figure 8C). Meanwhile, the CIBERSORT algorithm was performed to analyze infiltration abundance of 22 types of immune cells in patients. Correlations of these immune cells were shown in Figure 8D. As shown in Figure 8E

Relationship between the ORGs signature and immunotherapy response
Numerous studies have shown that immunotherapy is one of the effective treatments for patients with advanced ccRCC, and high expression level of immune checkpoint genes related to a better response to immunotherapy (3)(4)(5). The expression level of some wellknown immune checkpoint genes, including PD-1, PD-L1, CTLA4 and LAG3 were analyzed in patients from TCGA cohort. As shown in Figures 9A, B, considerably higher expression level of CTLA4, LAG3 and PD-1 were observed in patients in the high risk group, whereas the expression level of PD-L1 was no substantial change between the two groups. As predicted, the expression level of CTLA4, LAG3 and PD-1 were positively correlated with the risk score ( Figures 9C-F). Taken together, these findings illustrated that patients with high risk score had a higher sensitivity for immunotherapy responses.

Discussion
Oxidative stress, caused by the accumulation of large amounts of ROS, plays a vital role in the several phases of development of tumors, such as tumor initiation, progression and metastasis (10). In addition, accumulating evidences indicated that oxidative stress not only associated with tumor progression but also affected both the tumor microenvironment and the immunotherapy responses (23-25). In recent years, many studies focused on constructing prognostic signature for predicting prognosis and sensitivity to treatment of cancers (18)(19)(20)(21). Hence, constructing a risk prognostic signature using ORGs for predicting clinical outcome of ccRCC may be promising. We built a prognostic signature using 9 ORGs, and finally confirmed its clinical value. We identified two risk groups based on the calculated risk score in patients with ccRCC, and observed that survival of patients in the high risk group was obviously shorter than that of patients in the low risk group. The predicting efficiency of the signature was verified using ROC and independent prognostic analysis. The results suggested that the ORGs signature had accurately and independently predictive ability. Further, stratification analysis was performed in different subgroups based on clinical features, which suggested that patients in the high risk group had poor outcome in all subgroups. Meanwhile, we also found that patients in the high risk group were correlated with worse clinical features in terms of the tumor grade, T stage, M stage, N stage, and pathological stage. In addition, we built a nomogram model according to the results of multivariate independent prognostic analysis; high predictive performance was observed from the calibration graph. Our signature consisted of 9 ORGs, including ABCB1, AGER, E2F1, FOXM1, HADH, ISG15, KCNMA1, PLG, and TEK. ABCB1 was designed to protect cells from damage caused by xenobiotic and toxic substances, including chemotherapy drugs (26). It was reported that ABCB1 could confer resistance to chemotherapy (27). A recent review reported that inhibiting ABCB1 could restore cancer cell susceptibility to chemotherapy drugs (28). AGER was a member of immunoglobulin superfamily of cell surface receptors. AGER expression was associated with immune inflammatory response and cancerogenesis (29). Many evidences shown that the expression level and mutation rate of AGER were increased in multiple cancers, such as esophageal cancer, breast cancer, gastric cancer and endometrial cancer (30)(31)(32)(33). E2F1, a member of E2F transcription factor family, was up-regulated and confirmed as an oncogene in multiple human cancers, including hepatocellular cancer, breast cancer and gastric cancer (34-36). Shen et al. found that expression of E2F1 was up-regulated in ccRCC, and E2F1 knockdown inhibited the proliferation and metastasis of ccRCC cells (37). The transcriptional factor Forkhead Box M1 (FOXM1) was reported to belongs to the Forkhead box (FOX) transcription factor family (38). Previous studies revealed that FOXM1 transcription could increase the expression of multiple genes important for cancers progression (39)(40)(41). Jiang et al. found that FOXM1 regulated LINC01094 expression to promote ccRCC progression (42). HADH was a very important enzyme in the b-oxidation of fatty acid (43). It was reported that HADH expression was down-regulated in many tumors and its low expression was associated with tumors progression (44,45). ISG15 (IFN-stimulated gene), a ubiquitinlike protein, has been shown to be a tumor-related gene involved in tumors pathogenesis (46, 47). In cervical cancer, ISG15 expression was up-regulated and knockdown ISG15 inhibited proliferation and invasion of cervical cancer cell (48). KCNMA1 was one of the high voltage-activated channel conductance for potassium ions (49). Previous studies reported that KCNMA1 may be involved in various human tumorigenesis processes, such as prostate cancer (50), breast cancer (51), cervical cancer (52), and colorectal cancer (53). Plasminogen (PLG) acted as an important role in inhibiting tumor progression due to its ability to inhibit angiogenesis (54). PLG expression was down-regulated in ccRCC, and its low expression was associated with poor clinical outcome (55). TEK, also known as TIE-2, was a tyrosine kinase receptor for endothelial cells, with the ability to regulate of angiogenesis and remodeling (56). Liao et al. (57) found that TEK was low expressed in ccRCC compared with normal tissues and downregulation of TEK correlated with a poor clinical outcome which was also confirmed in previous studies (58,59).
Immune infiltration was closely related to tumor progression and immunotherapy responses. Thus, we further evaluated the relationship between the signature and immune infiltration status of ccRCC. Results from ssGSEA algorithm indicated that compared with individuals in the low risk group, individuals in the high risk group displayed more accumulation of immune cell infiltration and higher activity of immunity-related pathways. The results of CIBERSORT algorithm shown that infiltration abundance of B cells memory, Plasma cells, T cells CD8, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, NK cells activated, and Macrophages M0 were obviously higher in individuals with high risk score, while infiltration abundance of B cells naive, T cells CD4 memory resting, NK cells resting, Monocytes, Macrophages M2, Dendritic cells resting, and Mast cells resting were more higher in individuals with low risk score. Both T cells regulatory and T cells CD8 infiltration were related to adverse clinical outcome in individuals with ccRCC (60-62). Zhang et al. found that ccRCC patients with favorable prognosis presented relatively The relationship between risk score and expression of immune checkpoint genes in the TCGA cohort. higher enrichment levels of B cells naive, T cells CD4 memory resting, NK cells resting, Monocytes and macrophages M2 (63). Immunotherapy had proven to be effective and could significantly improve the prognosis for patients with advanced ccRCC (3)(4)(5). In our study, patients in the high risk group presented significantly higher expression level of CTLA4, LAG3 and PD-1. Of note, the above mentioned immune checkpoint genes expression were all positively associated with risk score of the signature. These results implied that compared with patients with low risk score, patients with high risk score displayed higher sensitivity for immunotherapy, and thus, the ORGs signature might have a potential in guiding personalized immunotherapy for patients with advanced ccRCC.
A recent article by Ma et al. also reported that an oxidative stress signature predicted clinical outcome in ccRCC and 4 oxidative stress genes (UCN, PLG, FOXM1, HRH2) were selected to construct a prognostic signature (64). The AUC of the 4 oxidative stress genes prognostic signature was 0.77 at 1 year, 0.70 at 3 years, and 0.71 at 5 years in the TCGA. Zhang et al. constructed a signature based on oxidative-stress related lncRNA in ccRCC (65). The signature consisted of 7 lncRNAs, including SPART-AS1, AL162586.1, LINC00944, LINC01550, HOXB-AS4, LINC02027, and DOCK9-DT. Wu et al. constructed a signature using mitochondrial genes related to oxidative stress to predict clinical outcome for ccRCC in the TCGA (66). Further analysis identified 6 prognostic-related mitochondrial genes, including ACAD11, ACADSB, BID, PYCR1, SLC25A27, and STAR. The AUC of the 6 mitochondrial genes prognostic signature was 0.736 at 1 year, 0.707 at 3 years, and 0.758 at 5 years in the TCGA. In our study, we constructed a prognostic signature using 9 four oxidative stress genes (ABCB1, AGER, E2F1, FOXM1, HADH, ISG15, KCNMA1, PLG, TEK). The AUC of the 9 oxidative stress genes prognostic signature was 0.81 at 1 year, 0.76 at 3 years, and 0.78 at 5 years in the TCGA.
Nonetheless, some limitations were presented in our study. First, our conclusions were obtained based on bioinformatic analysis, and multicenter and large-cohort clinical trials validation are needed to verify the robustness of the prognostic signature. Second, more in-depth experiments about the detailed biological functions of these genes that make up the ORGs signature are necessary to be investigated.

Conclusions
A novel 9-gene signature was built on the basis of the ORGs for predicting the clinical prognosis of ccRCC. It was proven that the prognosis signature model had a good and independent prediction performance. In addition, the prognosis signature might have a potential in predicting the responses of immunotherapy for patients with ccRCC, which could help clinicians to make immunotherapy decisions in order to achieve personalized treatment.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.