A Novel Ferroptosis-Related Pathway for Regulating Immune Checkpoints in Clear Cell Renal Cell Carcinoma

Ferroptosis is a novel form of cell death and plays a role in various diseases, especially tumors. It has been reported that ferroptosis is involved in the growth and progression of clear cell renal cell carcinoma (ccRCC); however, the specific molecular mechanisms are still unclear. In this study, we constructed a four-gene signature (FeSig) of ferroptosis-related genes via Cox regression analysis. ROC and survival analyses indicated that FeSig had good diagnostic and prognostic value. Further analysis revealed that ferroptosis was associated with tumor immunity in ccRCC. Next, weighted gene co-expression network analysis was performed to identify the potential regulatory mechanisms. Combined with correlation and survival analyses, the TAZ/WNT10B axis was identified as a tumor immune-related regulatory pathway. In conclusion, these findings suggest that ferroptosis is correlated with tumor immunity. The TAZ/WNT10B axis may be a novel biomarker and therapeutic target for immunotherapy in ccRCC.


INTRODUCTION
Renal cell carcinoma (RCC) is a common neoplasm that originates from renal tubular epithelial cells and accounts for 2-3% of all malignant tumors in adults (1). In the USA, it is estimated that there will be approximately 70,000 new cases and 14,000 deaths due to this cancer in 2020. The 5-year relative survival is 75.2% (2). Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype and accounts for over 70% of RCCs (1). Due to high resistance to conventional chemoradiotherapy, only limited therapies are currently available (3,4). At present, immunotherapy, as a novel treatment, is gradually being applied in ccRCC therapy (5). Immune checkpoint inhibitors (ICIs), such as PD-1 and CTLA-4 inhibitors (6,7), have improved clinical responses and quality of life for some patients. However, some patients are insensitive to ICIs (8). Therefore, it is necessary to further study the regulatory mechanisms of immune checkpoints in ccRCC.
Ferroptosis is a novel form of cell death caused by irondependent oxidative damage (9). Due to the failure of glutathione peroxidase (GPX4), a large number of reactive oxygen radicals (ROS) accumulate on membrane lipids (10). In recent years, many studies have verified that ferroptosis plays a vital role in degenerative diseases (i.e., Alzheimer's and Parkinson's diseases) (11,12), ischemia-reperfusion injury (13), and tumors (14). Zhang et al. showed that BAP1 inhibits tumor progression by promoting cellular ferroptosis (15). Furthermore, it has been reported that ferroptosis is closely associated with the therapeutic effect of immunotherapy. Wang et al. verified that CD8+ T cells enhanced the effect of immunotherapy by promoting ferroptosis of tumor cells (16). Interestingly, Lang et al. indicated that immunotherapy sensitizes tumors to radiotherapy by inhibiting SLC7A11 and reducing cystine uptake (17). However, the interaction between ferroptosis and immunotherapy in ccRCC is still unclear. In this study, we constructed a ferroptosis-related gene signature (FeSig) and analyzed the correlation between FeSig and immune checkpoints to identify novel therapeutic targets in ccRCC.

METHODS AND MATERIALS Data Download and Study Design
The gene expression data and clinical data were obtained from The Cancer Genome Atlas (TCGA-KIRC, https://www.cancer.gov/ tcga) and International Cancer Genome Consortium (ICGC-RECA-EU, https://dcc.icgc.org/) databases. An immunotherapy (immune checkpoint inhibitor, ICI) dataset of clear cell renal cell carcinoma (ccRCC) with clinical data and gene expression data was obtained from the study of Miao et al. (18). In addition, the study process is shown in a flow chart ( Figure 1).

Protein-Protein Interaction Network Analysis
The DE-FRGs were inputted into the STRING online tool (https://string-db.org) to construct the PPI network. Then, Cytoscape software (21) was used to analyze the PPI network.

Cox Regression Analysis and Construction of a Proportional Hazards Model
The "survival" package (https://CRAN.R-project.org/package= survival) was utilized to perform univariate and multivariate Cox regression analyses. The "glmnet" (22) and "survival" packages were used to perform least absolute shrinkage and selection operator (LASSO) regression analysis. Through integrated Cox analysis, four key FRGs were screened to construct the risk model (FeSig). The risk score (RS) formula was as follows: Where Coef (i) represents the coefficient, and X(i) represents the expression of selected genes.

Principal Component Analysis and T-Distributed Stochastic Neighbor Embedding Analysis
The "stats" and "limma" packages (23) were used to perform PCA. The "Rtsne" package (https://CRAN.R-project.org/ package=Rtsne) was utilized to perform t-SNE analysis. PCA and t-SNE analysis were both used to explore the distribution of different groups.

Survival Analysis and Receiver Operator Characteristic Curve Analysis
According to the median gene expression/risk score, ccRCC patients were divided into a high group and a low group. Then, survival curves of overall survival (OS) and disease-free survival (DFS) were drawn by the "survival" package in R and GraphPad software (version 7.0). A p-value <0.05 was considered statistically significant. Moreover, the "survivalROC" (https://CRAN.Rproject.org/package=survivalROC) package was used to generate a time-dependent ROC curve to evaluate the predictive value of the risk model.

Gene Set Enrichment Analysis and Gene Set Variation Analysis
All patients were divided into two groups (high group and low group) according to the median gene expression/risk score. GSEA (24) was performed to discover potential mechanisms and downstream signaling pathways. Moreover, the "GSVA" package (25) was utilized to find differential signaling pathways between the high group and the low group.

Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis
The "clusterProfiler" package (26) was used to conduct KEGG enrichment analysis of DEGs. The results were visualized by the "ggplot2" package (https://CRAN.R-project.org/package=ggplot2) in the R programme. A p-value <0.05 was selected as the cut-off point.

Weighted Gene Co-Expression Network Analysis
The "WGCNA" package (27) was applied to construct the weighted gene co-expression network of DEGs. First, the TCGA-KIRC cohort was evaluated via sample clustering to detect outliers with a height cut-off point = 60. Then, Pearson's correlation was calculated between each of the gene pairs. Second, a matrix of adjacencies was constructed according to Pearson's Ferroptosis-Related Pathway Regulates Immune-Checkpoints correlation. Then, the adjacencies achieved scale-free topology based on the soft threshold power b ( Figure 7A). Third, the adjacencies were transformed into a topological overlap matrix (TOM). Genes with a high absolute correlation were clustered into the same module. Finally, combined with clinical traits, we calculated the correlation between the module eigengenes (MEs) and clinical traits to screen the clinically significant modules. After that, the correlation between genes and clinical traits (cor.geneTraitSignificance) and the correlation between genes and MEs (cor.geneModuleMembership) were conducted to identify hub genes. In this study, cor.geneTraitSignificance >0.2 and cor.geneModuleMembership >0.6 were selected as the cutoff criteria.

Human Clinical Specimens
A total of 19 pairs of ccRCC samples and adjacent normal samples (4 cm away from the margin of the tumor tissues) were obtained from Wuhan Union Hospital between 2018 and 2020. The study was approved by the Human Research Ethics Committee of Huazhong University of Science and Technology (HUST), and all patients signed the informed consent.

Statistical Analysis
In this study, all data are presented as the mean ± SD. SPSS (version 22.0) and GraphPad Prism (version 7.0) were used to analyze the data. Student's t-test, Mann-Whitney test, and Pearson's c 2 test were used to conduct statistical analyses. A p-value <0.05 was considered statistically significant.

Identification of DE-FRGs and Construction of a Proportional Hazards Model
The gene expression data were obtained from the TCGA-KIRC cohort. According to the cut-off criteria, 46 DE-FRGs were identified ( Figure 2A). The protein-protein interactions of DE-FRGs were shown in Figure 2B. Further analysis indicated that there were 20 DE-FRGs with good prognostic value ( Figure 2C).
The correlation was close and high between these DE-FRGs ( Figure 2D). Then, LASSO and multivariate Cox regression analyses were applied to construct a prognostic model. A fourgene signature (FeSig) was identified ( Figures 2E-G). The risk score = 0.49*expression of BID + 0.70*expression of TAZ + 0.09*expression of MT1G + 0.45*expression of SLC7A11.
The Diagnostic and Prognostic Value of the Four-Gene Signature PCA and t-SNE analysis proved that FeSig could significantly divide patients into different risk groups ( Figures 3A, B). In the TCGA-KIRC cohort, patients were divided into two groups (high-risk group and low-risk group) according to the median risk score ( Figure 3C). As shown in Figure 3D, patients with a high-risk score had a higher probability of death earlier than those with a low risk score. Further Cox regression analysis revealed that the risk score was an independent prognostic factor ( Figures 3E, F). Time-dependent ROC curves indicated that the risk score had good predictive performance in both the TCGA-KIRC ( Figure 3G) cohort and the ICGC-RECA-EU ( Figure 3H) cohort. Moreover, survival analysis also uncovered that the highrisk group predicted poor OS ( Figures 3I, J). These findings suggested that FeSig had good diagnostic and prognostic value.

FeSig Is Closely Correlated With
Immune-Related Pathways GSEA, GSVA, and KEGG analyses were performed to identify potential downstream signaling pathways of FeSig. As shown in Figure 4A, FeSig was associated with immune cell-related pathways. Analogously, GSVA also verified that there were many differentially enriched immune cell-related pathways between the high-risk group and the low-risk group ( Figures  4B, C). Then, 314 DEGs were identified via DESeq2 and edgeR package in the R programme ( Figures 4D, E). Further KEGG pathway enrichment analysis indicated that these DEGs were mainly enriched in immune-related pathways ( Figure 4F, Ras signaling pathway, PPAR signaling pathway, and IL-17 signaling pathway). To further study the correlation between FeSig and immune status, we quantified the enrichment scores of diverse immune cell subpopulations, related functions, or pathways with ssGSEA. As shown in Figure 5A, nine immune cell subpopulations (CD8+ T cells, macrophages, mast cells, neutrophils, T helper cells, Tfhs, Th1 cells, Th2 cells, and TILs) were clearly different between the high-risk group and the lowrisk group. For immune-related functions, seven pathways (CCR, checkpoint, cytolytic activity, inflammation promotion, parainflammation, T cell coinhibition, and T cell costimulation) had higher scores in the high-risk group, and only the type II IFN response had higher scores in the low-risk group ( Figure 5B). These findings revealed that FeSig is significantly associated with the regulation of tumor immunity. and immune checkpoint pathways. As shown in Figure 6A, we found that PD-1, CTLA4, LAG3, and TIGIT were upregulated in the high-risk group and TIM-3 was downregulated in the high-risk group. Moreover, correlation analysis indicated that the riskScore was significantly positively correlated with PD-1, CTLA4, LAG3, and TIGIT ( Figures 6B-H). These findings suggested that FeSig may be a biomarker of the response to immunotherapy. Therefore, the predictive value of FeSig was tested in an immunotherapy dataset of ccRCC. Survival analysis verified that there was a significant difference  (log-rank p = 0.038) between the high-and low-risk groups of patients with anti-PD-1 monotherapy. Patients in the highrisk group had longer OS than those in the low-risk group ( Figure 6I). The ROC curve of ICI response ( Figure 6J) revealed that the diagnostic value of FeSig (AUC = 0.603) was better than that of PD-1 (AUC = 0.544).

Identification of the TAZ/WNT10B Axis as an Immune Checkpoint Regulatory Pathway in ccRCC
The WGCNA algorithm was performed to find further downstream targets offour-gene signature (FeSig). As shown in Figure 7A, b = 5 was selected as the soft threshold power in this study. The DEGs  between high-risk group and low-risk group were divided into four modules according to the expression level of each gene ( Figure 7B). Combined with the clinical information, we found that the MEturquoise module was significantly positively associated with grade (R = 0.23, p = 2e-07), TNM stage (R = 0.29, p = 3e-11), risk (R = 0.66, p = 2e-66), and ImmuneScore (R = 0.2, p = 6e-06) ( Figure 7C). According to cor.geneTraitSignificance >0.2 and cor.geneModuleMembership >0.6, CPNE7, WNT10B, ADAMTS14, and RUFY4 were regarded as downstream hub genes of FesSig ( Table 1, Figures 7D-H). Moreover, we analyzed the correlation among four-gene signature (FeSig) (BID, TAZ, MT1G, and SLC7A11), downstream hub genes (CPNE7, WNT10B, ADAMTS14, and RUFY4), and immune checkpoints (PD-1, CTLA4, LAG3, and TIGIT) to further discover potential molecular mechanisms of tumor immunity. According to the correlation value, BID, TAZ, CPNE7, WNT10B, and RUFY4 were selected for subsequent analysis ( Figure 7I). Survival analyses of OS and DFS proved that WNT10B had better prognostic prediction performance than CPNE7 and RUFY4 ( Figures 8A-F). Similarly, TAZ had better prognostic value than BID ( Figures 8G, H). Therefore, TAZ and WNT10B were selected for further study. In addition, the expression of TAZ was clearly positively correlated with WNT10B ( Figure 8I, R = 0.66, p = 2e-67). Interestingly, GSEA verified that the low TAZ group was enriched in WNT signaling pathway ( Figures 8J, K). These findings indicated that TAZ might regulate tumor immunity through mediating WNT10B (WNT signaling pathway) in ccRCC.

TAZ/WNT10B Was Closely Correlated With TNM/Grade Stage and UpRegulated in ccRCC Tissue
We further analyzed the correlation between the expression level of TAZ/WNT10B and TNM/Grade stage and verified their mRNA expression level in ccRCC tissues. As shown in Figures  9A, B, TAZ and WNT10B were both elevated and positively correlated with TNM/Grade stage in ccRCC. Moreover, qRT-PCR assay also indicated that TAZ and WNT10B were upregulated in ccRCC tissues ( Figure 9C).

DISCUSSION
Recently, immunotherapy has emerged as a promising approach for cancer treatment, which mainly includes non-specific immunostimulation, immune checkpoint inhibitors (ICIs), tumor vaccines, and adoptive cellular immunotherapy (28). In metastatic RCC (mRCC), immune checkpoint inhibitors have changed the treatment paradigm because most patients with newly diagnosed mRCC are now treated with these medicines. It has been reported that immune checkpoint inhibitors have provided significant clinical benefit for mRCC patients in multiple clinical trials (29). However, there are still some patients without good effects due to adverse reactions and drug resistance to ICIs (30). Therefore, it is necessary to study the molecular mechanisms of drug resistance and reduce the adverse reactions to ICIs. At present, some studies have reported that the induction of ferroptosis enhances the effect of ICIs (16). Interestingly, many ferroptosis-related genes are abnormally expressed in ccRCC and might be potential therapeutic targets. However, the correlation between ferroptosis and immune checkpoints in ccRCC is still unclear.
In this study, a ferroptosis-related risk model (FeSig) was constructed through Cox regression analysis. We found that FeSig had good diagnostic/prognostic value and was closely associated with immune checkpoints. Moreover, patients with high-risk scores had better OS than those with low-risk scores  after PD-1 inhibitor treatment, which indicated that FeSig was a potential prognostic biomarker for immunotherapy. Then, the TAZ/WNT10B axis was identified as a potential regulatory pathway of immune checkpoints via WGCNA and correlation analysis.
TAZ (Tafazzin) is a transcriptional regulator and plays a vital role in tumorigenesis and tumor progression of most solid tumors (31). TAZ stimulates cell proliferation by regulating DNA duplication and mitosis (32). It has been reported that TAZ maintains plasticity in cell-ECM adhesion and favors  cytoskeletal remodeling to promote tumor metastasis (33). TAZ is also regarded as an important regulator of ferroptosis. Two research groups recently reported that TAZ promoted tumor cell ferroptosis by regulating members of the NOX family (34,35).
Furthermore, many studies have verified that TAZ is closely involved in tumor immunity and the microenvironment. TAZ mediates the expression of tumor-secreted factors to drive the differentiation and recruitment of immune suppressive cells, such as tumor-associated macrophages (TAMs) (36) and regulatory T cells (Tregs) (37). In addition, TAZ promotes immune evasion of tumor cells by regulating the expression of immune checkpoints. Feng et al. indicated that tumor cellderived lactate activated the TAZ/PD-L1 axis to enhance tumor evasion from the immune response (38). Similarly, a recent study also reported that the TAZ/YAP/TEAD signaling pathway increased PD-L1 promoter activity and induced immune evasion of tumors (39). In our study, we found that the expression of TAZ was highly positively associated with PD-1, which suggested that TAZ was a potential therapeutic target of immunotherapy in ccRCC.
To study the specific molecular mechanisms of TAZ, WGCNA and GSEA were performed and indicated that WNT10B is a potential downstream target of TAZ. WNT10B (Wnt Family Member 10B) is a regulator encoding secreted proteins and activating the Wnt signaling cascade (40). It has been reported that the Wnt signaling pathway is closely involved in regulating immune checkpoints. Notably, the expression of PD-L1 has been proven to be regulated by MYC, which is a well-documented target of the Wnt signaling pathway (41). Moreover, inhibiting the Wnt/bcatenin axis can promote antitumor immunity by suppressing PD-L1 expression (42). In this study, we also found that TAZ was significantly positively correlated with WNT10B and that high WNT10B predicted poor OS/DFS in ccRCC. These findings suggested that the TAZ/WNT10B axis might regulate tumor immunity by activating the WNT signaling pathway.
In conclusion, the TAZ/WNT10B axis (ferroptosis-related pathway) is regarded as a tumor immune-related regulatory pathway via integrated bioinformatics analysis. Further analysis revealed that immune checkpoints are potential targets of the TAZ/ WNT10B pathway. Therefore, the TAZ/WNT10B axis is expected Ferroptosis-Related Pathway Regulates Immune-Checkpoints to be a novel therapeutic target of immunotherapy in ccRCC. However, the specific mechanisms still require further research.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Human Research Ethics Committee of Huazhong University of Science and Technology (HUST). The patients/ participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
XZ designed this study. TX and SG performed data collection and analysis. TX, SG, HR, and JL performed the majority of the experiments. TX and SG wrote the manuscript and contributed to preparing and making figures and tables. YL, DL, and JT collected the clinical samples and managed the clinical data. JT and JS reviewed the relevant literature. HY and XZ provided conceptual advice and critically reviewed the manuscript. All authors contributed to the article and approved the submitted version.