Construction and Validation of a Novel Immune Checkpoint-Related Model in Clear Cell Renal Cell Carcinoma

Background With the highest mortality and metastasis rate, kidney renal clear cell carcinoma (KIRC) is one of the most common urological malignant tumors and not sensitive to chemotherapy and radiotherapy. Immunotherapy, which proves to be effective and a big progression, such as PD-1/PD-L1 inhibitors, is not sensitive to all KIRC patients. To predict prognosis and immunotherapy response, a novel immune checkpoint gene- (ICG-) related model is essential in clinics. Methods From the public database-downloaded dataset, a novel ICG-related model for predicting prognosis and immunotherapy response in KIRC patients was built up and verified with R packages and Cox regression analysis. The Kaplan-Meier curve was plotted. Results 39 ICGs were identified to have different expression in KIRC patients and enriched in immune-related biological pathways and activities. Three ICGs (CTLA4, TNFSF14, and HHLA2) were screened to generate KIRC-ICG model. The KIRC-ICG model was verified to be effective. With conducting KIRC-SYS model, KIRC-ICGscore was verified to be an independent factor regardless of age, gender, stage, grade, and TNM stage. Compared to the ICG-low subgroup, the ICG-high subgroup had more immune activities. KIRC-ICGscore was significantly positively correlated with the expression of Treg markers. KIRC-ICG model could also be reliable to predict immunotherapy response. Conclusion The KIRC-ICG model was reliable to predict prognosis and immunotherapy response for KIRC patients and could be an independent factor regardless of clinical characteristics.


Introduction
Taking up 80% of renal carcinoma, kidney renal clear cell carcinoma (KIRC) is one of the most common urological malignancies [1]. Patients with renal clear cell carcinoma have the highest mortality and metastasis rate [2][3][4] and are not sensitive to chemotherapy and radiotherapy [5]. Combined with traditional cancer treatment, tumor immunotherapy has high efficacy and safety, which can significantly improve the survival rate of patients [6].
Despite significant progression in immunotherapy [7], there are still a proportion of KIRC patients who are not sensitive to these inhibitors due to immunosuppression [8].
Tumor immunosuppression refers to the inhibition of the host immune response to tumor antigens under the action of some factors, resulting in tumor cells escaping from the surveillance of the body's immune system, rapidly dividing and proliferating in the body, and accelerating the deterioration of tumors [9,10]. These factors include a series of molecules expressed on immune cells that regulate the extent of immune activation, called immune checkpoints. For example, PD-1, a type I transmembrane glycoprotein expressed mainly on the surface of T cells, inhibits the activation and proliferation of T cells by binding to PD-L1 [11]. Nevertheless, due to the diversity of immune checkpoints, different types of cancer can express different immune checkpoints. To predict prognosis and immunotherapy response, a new immune checkpoint gene-(ICG-) related model is essential in clinics for KIRC patients.
To address this issue, we identified ICG-related differential genes (DEGs) and constructed a new KIRC-ICG model. Meanwhile, we also verified the KIRC-ICG model performance in prognosis and immunotherapy response effectiveness. It would provide a new vision for improving accuracy in immunotherapy in KIRC patients.

Data Collection.
The KIRC data (RNA-seq data and clinical data) from the TCGA database were downloaded from the University of California Santa Cruz's official webpage (UCSC database) (https://xenabrowser.net/datapages/). TCGA-KIRC dataset included RNA transcript of 531 KIRC patients and 72 paracarcinoma tissues, in which we obtained 515 KIRC patients with survival time greater than 30 days. To supplement the calculations, we also obtained the RNAseq data from normal people in The Genotype-Tissue Expression (GTEx) databases. In addition, we got the RECA-EU dataset with expression profiling of KIRC patients as the validation cohort from the International Cancer Genome Consortium (ICGC) database. The complete dataset of 56 renal cell carcinoma patients with immune therapy information was downloaded from the IMvigor210 (http:// research-pub.gene.com/IMvigor210CoreBiologies) [12]. 47 immune checkpoint genes (ICGs) were collected from published literature (Table 1) [13].
2.2. Identification of ICGs. The "limma" package was applied to filter immune checkpoint genes in both TCGA-KIRC and RECA-EU datasets by setting p < 0:05 and jlog 2 FoldChangej > 1 as the threshold. The GeneMANIA was used to explore the protein network in these ICGs from multiple perspectives [14].

Construction of Prognostic Risk
Model. Univariate Cox regression analysis was used to identify ICGs which were related to overall survival (OS). Next, through a multivariate Cox analysis (stepwise model), we further filtered these ICGs affecting prognosis and built up the prognostic risk model with the lowest Akaike information criterion (AIC) value. AIC is a standard for assessing the complexity of statistical models and measuring the goodness of fit. Based on the coefficient and expression value of ICGs, the risk score of each sample was estimated. The equation was as follows: To evaluate the model, the "survivalROC" package was used to plot time-dependent ROC curves. Based on the KIRC-ICGscore, we classified KIRC patients (n = 515) into the ICG-high (n = 258) and ICG-low subgroups (n = 257).

Analysis of Tumor Immune
Microenvironment (TIME) in KIRC-ICG Subgroups The "GSVA" and "Estimate" packages were used to analyze the differences in the TIME status between the two KIRC-ICG subgroups [15,16]. The Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard .edu/) was applied to predict the effect of immunosuppressive therapy on the two KIRC-ICG subgroups [17].

Immune Clustering Based on ssGSEA
The ssGSEA algorithm was based on the 29 immuneassociated gene sets which represented diverse immune cell types, functions, and pathways [18]. We applied the ssGSEA algorithm via R packages "GSVA" to comprehensively assess the immunologic characteristics of every sample included in the TCGA-KIRC set and performed an unsupervised hierarchical clustering algorithm to divide these samples into Immune_H and Immune_L clusters based on the ssGSEA results.
With the median value of the KIRC-ICGscore, we classified KIRC patients (n = 515) into the ICG-high (n = 258) and ICG-low subgroups (n = 257). The distribution of the KIRC-ICGscore and survival status between the two KIRC-ICG subgroups were shown in Figure 2(c). There were more death samples with higher KIRC-ICGscore (Figure 2(d)). The Kaplan-Meier analysis confirmed a lower OS rate in the ICG-high subgroup (p < 0:05, Figure 2(e)). The AUCs of time-dependent ROC curve at 1, 3, and 5 years were 0.718, 0.696, and 0.751, which indicated that the KIRC-ICG model was reliable for predicting the prognosis of KIRC patients ( Figure 2(f)). To illustrate the underlying mechanism of the     survival difference between the ICG-high and ICG-low subgroups, GSEA was applied. The results related to KEGG and GO pathways showed that the ICG-high subgroup had more immune activities (Figures 3(a)-3(f)).

Clinical Application of the KIRC-ICG Model.
To evaluate the KIRC-ICGscore in the clinic, six characteristics (age, gender, stage, grade, TNM stage, and KIRC-ICGscore) were enclosed to build up the systematic prognostic risk model named KIRC-SYS model by conducting univariate Cox regression and multivariate Cox regression analyses. As shown in Figures 4(a)-4(b), the KIRC-ICGscore was identified as an independent risk factor for OS. And the AUC value for the KIRC-ICGscore was 0.771 (Figure 4(c)), which was superior to other clinicopathological parameters, espe-cially the TNM stage. All these features indicated that the KIRC-ICGscore was systematically dependable in predicting the prognosis of KIRC patients. Moreover, we drew a distribution heat map of clinicopathological characteristics and three ICGs between the two KIRC-ICG subgroups (Figure 4(d)) and discovered that TNM stage, stage, grade, and survival status (fustat) were different between the two KIRC-ICG subgroups (all p < 0:05). With the KIRC-ICGscore as a continuous variable, we further demonstrated that the high KIRC-ICGscores were linked to poorly differentiated tumors, advanced TNM stages, and death ( Figure 5). It further revealed that the high KIRC-ICGscores might be associated with the progression of KIRC. We performed the Kaplan-Meier curves for OS in subsets of different clinical characteristics such as age,   (Figures 6(a)-6(g), respectively). We observed the ICG-high subgroup was consistently relevant to adverse prognosis in different subsets except for N1, which might mean that the KIRC-ICG model could independently predict the survival outcome of KIRC patients regardless of clinical features.

The Impact of the KIRC-ICG Model on Immune Therapy
Response. To check the impact of the KIRC-ICG model on immune therapy, we applied the TIDE algorithm. With correlating gene expression profiles of cytotoxic T lymphocyte (CTL) markers (CD8A, CD8B, GZMA, GZMB, and PRF1) with T-cell characteristics, the TIDE algorithm predicts immunotherapy of tumor patients. The high TIDE score implies that the patient is not responding well to immune therapy due to T-cell dysfunction or more exclusion of T-cell infiltration. The results revealed that the ICG-high subgroup was characterized by significantly higher TIDE scores (Figure 8(a)), although the KIRC-ICGscore had a weak correlation with the TIDE score (r = 0:33, p < 0:01, Figure 8(b)), which mean that the ICG-high subgroup would have a worse response on immune therapy. By analyzing IMvigor210 data of 64 KIRC patients which were classified as progressive disease (PD), stable disease (SD), partial response (PR), and complete response (CR) according to the treatment response of atezolizumab (PD-L1 inhibitor), we got the result that the CR/PR group had significantly lower KIRC-ICGscore than the SD/PD group (Figure 8(c)). The Kaplan-Meier curve confirmed that the ICG-high subgroup was correlated with worse OS after atezolizumab treatment (Figure 8(d)).

Construction and Validation of Immue
Clustering of KIRC Based on the ssGSEA algorithm, KIRC patients were analyzed using 29 immune gene sets. The hierarchical clustering algorithm classified KIRC patients (n = 515) into the Immunity_H (n = 258) and Immunity_L clusters (n = 257, Figure 9(a)) based on ssGSEA results. The ESTIMATE algorithm was used to score the TIME of each patient. The results showed that Immunity_H cluster had higher immune and ESTIMATE scores but lower tumor purity status (p < 0:05, Figure 9(b)). Then, GSVA enrichment analysis revealed that immune-related pathways such as T-cell receptor signaling pathway, cytokine-cytokine receptor interaction, and primary immunodeficiency were significantly enriched in the Immunity_H cluster (Figure 9(c)). The above results indicated that the Immunity_H cluster had higher immune activities and more immune cell infiltration.
The Kaplan-Meier analysis demonstrated that patients in the Immune_H cluster had a poorer prognosis than Immune_L cluster (Figure 9(d)). Meanwhile, cytokines (IL-4 and IL-10) and checkpoints (LAG3, CTLA4, TIGIT, and PDCD1) related to immune suppression were highly expressed in the Immunity_H cluster. The FOXP3 expression was also upregulated in the Immunity_H cluster, which may be attributed to the infiltration of Tregs (Figure 9(e)).

Discussion
Conventional radiation therapy and cytotoxic chemotherapy have little efficacy in advanced KIRC [5]. Immunotherapy is proven to be a big progression, such as PD-1/PD-L1 inhibitors. PD-1/PD-L1 is one of the immune checkpoint genes that could impair the killing ability of T cells [20]. However, not all patients with PD-L1 overexpression can benefit from PD-L1 inhibitors [21,22]. Novel ICGs have been discovered [23]. In our study, we collected 47 ICGs from published literature. Inside the 47 ICGs, cytotoxic T lymphocyte-associated antigen 4 (CTLA-4) belongs to the same immunoglobulin supergene family as CD28. It is mainly expressed on activated T cells and regulatory T cells. With ligands of B7-1 and B7-2 on APCs, CTLA-4 is one of the most important T cell costimulatory receptors [24], and it is considered as a negative regulator of T cells [25]. High expression of CTLA4 was reported to associate with poor overall survival (OS) and might serve as prognostic biomarkers in KIRC patients [26].
Tumor necrosis factor superfamily member 14 (TNFSF14) is mainly expressed on activated T cells, NK cells, and immature dendritic cells (DCs) [27], and its two primary    [28]. TNFSF14-LTβR interaction is necessary for lymphoid structure development and maintenance [29]. However, it is expressed on lymphocytes, NK cells, and smooth muscle; HVEM functions is an important T-cell costimulatory factor [30,31]. TNFSF14-HVEM can induce strong antitumor immunity in a therapeutic context [32]. TNFSF14 was reported to be a prognostic factor as one of the TNF-related signatures and could assess the efficacy of immunotherapy [33].
Human endogenous retrovirus-H long terminal repeatassociating protein 2 (HHLA2) is a new member of the B7 family and could mediate costimulation by interacting with TMIGD2 [34]. Patients with HHLA2 overexpression exhibited a worse clinical prognosis in KIRC [35]. HHLA2 overexpression is also found in other carcinomas like intrahepatic cholangiocarcinoma [36], stomach cancer [37], and colorectal cancer [38] and is associated with poor prognosis. However, there were also reports that HHLA2 was highly expressed in KIRC and predicted a favorable survival outcome in KIRC [39,40], which was consistent with our research. According to the report, HHLA2 has a costimulatory receptor TMIGD2 (transmembrane and immunoglobulin domain containing 2) and a coinhibitory receptor KIR3DL3 (killer cell Ig-like receptor, three Ig domains, and long cytoplasmic tail), which endows it with both immunostimulant and immunosuppression functions in cancer development [41]. These opposing functions emphasized that the expression and distribution of HHLA2 and its receptor may determine the immune response in the tumor microenvironment, thereby influencing prognosis.
Together with three ICGs (CTLA4, TNFSF14, and HHLA2), here, we successfully established the predicting KIRC-ICG model for the first time. Based on the IMvi-gor210 cohort, we proved that patients in the ICG-high subgroup were less sensitive to ICGs. Meanwhile, we applied the TIDE algorithm to speculate the reason why the ICG-high subgroup was less sensitive to ICIs from the perspective of the CTL. The result showed the TIDE score in the ICGhigh subgroup that was higher than the ICG-low subgroup, which indicated that the patients are not responding well to immune therapy due to more T-cell dysfunction or more exclusion of T-cell infiltration. Notably, the KIRC-ICGscore only weakly correlated with TIDE scores. We hypothesized that in addition to the CTL, there may be other immune cells that influence the immunotherapy response by interacting with CTL. Even though there are limitations on the bioinformatic model for prediction prognosis and immunotherapy response, it could give us a different perspective on the function of the corresponding risk factors.
Contrary to the general trend of tumor immunity, patients with poorer prognosis had higher immune activities and more immune cell infiltration (including Tregs) in our study. To further verify this opposite immune characteristic, the ssGSEA algorithm was applied. By the ssGSEA algorithm, KIRC patients were classified into the Immunity_H and Immunity_L clusters. And the Immune_H cluster had shorter OS and more Treg cell infiltration than those in the Immune_L cluster, which was the same as the immune characteristic of the KIRC-ICG model. Some reports have proved that high tumoral immune activity does not necessarily indicate a better prognosis in certain cancer types [42,43]. The above trend in cancer immunology is uncommon, and it may be influenced by cancer-type-specific factors. According to the report, high immune activity in uveal melanoma might induce epithelial-mesenchymal transition/endothelial-mesenchymal transition and loosen the blood barrier structure, further leading to a poor prognosis [44]. However, the underlying mechanisms of the special immune characteristic of KIRC are unclear and needed to be further explored.

Conclusion
In conclusion, we created a new KIRC-ICG model which was reliable to predict prognosis and immunotherapy response for KIRC patients and could be an independent factor regardless of clinical characteristics. However, this study still has some limitations. To further validate the model, bench experiments are needed in the future.