Multidimensional Analyses of Tumor Immune Microenvironment Reveal the Possible Rationality of Immunotherapy and Identify High Immunotherapy Response Subtypes for Renal Papillary Cell Carcinoma

Abstract

Kidney renal papillary cell carcinoma (KIRP), the second most common subtype of renal cell carcinoma, still lacks effective treatment regimens for individualized immunotherapy because of the heterogeneity of its elusive immune microenvironment. Therefore, we aimed to comprehensively evaluate the immune microenvironment of KIRP by using the computational biology strategy to analyze the expression profile data of 289 KIRP patients obtained from The Cancer Genome Atlas database. Based on multidimensional, multi-omics bioinformatics analysis, we found that the tumor of patients with KIRP exhibited “hot” tumor characteristics but the CD8+ T cells in the tumor tissues did not limit tumor progression. Thus, patients with KIRP may realize higher clinical benefits by receiving treatment that can reverse CD8+ T-cell exhaustion. Among them, C1 and C3 immune subtypes could realize the best efficacy of reversing CD8+ T-cell exhaustion. Moreover, CCL5 and FASLG expression may be related to the formation of the immunosuppressive microenvironment in the tumors of patients with KIRP. In conclusion, the immune microenvironment landscape presented in this study provides a novel insight for further experimental and clinical exploration of tailored immunotherapy for patients with KIRP.

Introduction

Renal cell carcinoma (RCC) is one of the most common types of kidney cancer in humans. In the past two decades, the incidence of kidney cancer has been increasing (accounting for 2%–3% of all new tumor cases) and RCC accounting for approximately 85% of kidney cancer cases (1). Kidney renal papillary cell carcinoma (KIRP) is the second most common type of RCC, accounting for 10%–20% of all RCC cases (2). Although considerable progress has been made in the diagnosis and treatment of RCC, the clinical outcome of RCC is still not satisfactory (3, 4). As an important subtype of kidney cancer, the preferred treatment for local and locally advanced KIRP is surgical resection. Although multiple drugs have been used for treating advanced and metastatic KIRP, including anti-vascular endothelial factor drugs and mTOR inhibitors, effective and individualized immunotherapeutic programs for treating KIRP have still not been developed (2).

Tumor-infiltrating immune cells (TIICs) form an ecosystem that regulates cancer progression in the tumor microenvironment and have exhibited potential prognostic value (5). The most studied TIICs are tumor-infiltrating lymphocytes, among which cytotoxic CD8+ T cells can prevent tumor growth (6). However, antigen-epitope-specific CD8+ T cells present in the tumor in situ or the peripheral blood are often inactivated first (7, 8). The combination of the CTLA4–B7 complex on the surface of T cells inhibits the activation of T cells, leading to the inhibition of T-cell proliferation (9). A previous study has shown that TIM-3 and PD-1 have a potential synergistic effect and induce the dysfunction of tumor-infiltrating CD8+ T cells (10). Therefore, these lead to the exhaustion of CD8+ T cells in the tumor microenvironment and, thus, the loss of antitumor effects. At present, targeted drugs for PD-1 and CTLA4 have been approved treating various cancers including KIRC and have achieved good results (11–15). The size and number of intraperitoneal and retroperitoneal metastases were considerably reduced in a patient with type 2 recurrent metastatic KIRP after 5 months of treatment with nivolumab monoclonal antibody (16), suggesting the potential feasibility of restoring the antitumor ability of CD8+ T cells to treat KIRP. However, a comprehensive evaluation of the tumor immune microenvironment in patients with KIRP has still not been performed, which has severely limited the further exploration of the feasibility of immunotherapy in these patients.

Therefore, in this study, we comprehensively evaluated the tumor immune microenvironment of patients with KIRP and identified key genes and T-cell subsets that are closely related to the tumor microenvironment of patients with KIRP and can be used as immunotherapeutic targets or markers. Furthermore, we identified subtypes with different immunotherapy responses of KIRP patients through an unsupervised clustering method and further demonstrated the prospect of immunotherapy based on reversing CD8+ T-cell exhaustion in KIRP patients.

Materials and Methods

Data Sources and Differential Analysis of mRNAs

The RNA-seq matrix files (count format) and clinical information of 289 KIRP samples and 32 peritumoral kidney samples were downloaded from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/) (Tables 1 and 2). The mRNA-seq matrix files were extracted for the next differential analysis. The R package “edgeR” was applied to screen differentially expressed mRNAs between normal and tumor tissues. Next, the p-value was calculated by the FDR-corrected method. The mRNAs with fold change >2 (| log2 fold-change | >1) and p < 0.05 were filtrated as differentially expressed genes.

Evaluation of Immune Infiltration

The CIBERSORT method is used to evaluate the proportion of immune cells based on the standardized gene expression files. This method has been validated in gene expression file studies measured by gene chips (17). We used the CIBERSORT to calculate the proportion of 22 immune cells (B cells naïve, B cells memory, plasma cells, T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), T cells gamma delta, NK cells resting, NK cells activated, monocytes, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, eosinophils, and neutrophils). Samples with p < 0.05 indicated that the proportion of immune cells calculated by CIBERSORT was correct. In addition to CIBERSORT, we used the tumor immune estimation resource (TIMER) database (https://cistrome.shinyapps.io/timer/) to assess the abundance of six immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and myeloid dendritic cells) for the TCGA and GEO cohorts (18).

Construction of Weighted Co-Expression Network

The regulation network method is widely applied to analyze the gene expression data. When compared with the node-based method, the regulation network method pays attention to both the difference in and the correlation between the gene expression profiles (19). The weighted gene co-expression network analysis (WGCNA) was superior to several other methods in constructing the co-expression networks (20). To determine the gene groups with similar expression patterns, we applied the previously reported R package “WGCNA” to construct a weighted co-expression network for the TCGA cohort (21). First, the outlier samples were filtered to reduce the differences. The soft threshold method was selected to construct the correlation network so that the adjacency matrix becomes a continuous value between 0 and 1. The network thus constructed was closer to the real biological network state. Then, the BlockwiseModules function was applied to build scale-free networks, after which the module division analysis and dynamic tree-cutting algorithm were applied to group the genes with similar expression patterns and define the modules (22). In order to obtain a large number of genes for each module, we selected a dividing line (0.25) to merge some of the modules with a high correlation.

Identification of Important Modules Related to Immune Cells

To determine the important modules related to T cells in the TCGA cohort, we calculated the value of module eigengenes (MEs), gene significance (GS), and module significance (MS). MEs are the main components of gene principal component analysis in the module. By calculating the correlation between MEs and clinical information, vital relevant modules could be determined. GS is a log10 transformation of p-value in the linear regression between the gene expression and clinical information, which represents the correlation between genes and samples. MS is the average GS in the module. The module with the highest correlation with the CD8+ T cells was selected as the module to be analyzed. The enrichment pathways of the key modules were obtained using the Metascape (http://metascape.org/gp/index.html).

Construction of the PPI Network and Screening of Its Hub Genes

For the TCGA cohort, the protein–protein interaction (PPI) network in the hub module was constructed using default parameters of the STRING database (23) and then visualized using Cytoscape. The hub genes in the network were calculated using the MCC algorithm in CytoHubba. When compared with the 10 other previously reported methods, MCC can more accurately identify the central objects in the network (24).

Screening of Hub Genes in the WGCNA Network

For the TCGA cohort, the node with weights >0.4 between the two nodes in the key module was input into the Cytoscape, after which the hub genes in the co-expression network were calculated using the MCC algorithm of CytoHubba. The genes highly correlated with certain clinical features in the modules were identified as the more important genes in the module; therefore, we defined genes that simultaneously satisfied the module membership (MM) >0.9 and gene significance (GS) >0.6 as the hub genes in the WGCNA network.

Identification of Real Hub Genes in the Key Module

We identified the genes obtained from the intersection of hub genes of the PPI network and the WGCNA network as the hub genes of the key module in the TCGA cohort. The selection of the intersection and drawing of the Wayne diagram was performed using the online tool VENNY (Https://bioinfogp.cnb.csic.es/tools/venny/index.html) (25). After obtaining the true hub genes, we used the R package “corrplot” to draw the correlation heat map between these hub genes.

Immunohistochemistry and Real-Time Quantitative PCR Analyses

In total, 30 KIRP tissues and 10 adjacent normal tissues were collected from 30 patients who underwent partial nephrectomies or nephrectomies at The First Hospital of China Medical University. All patients provided their written informed consent. Total RNA was extracted using the TRIzol reagent (Waltham, Massachusetts, USA). The PrimeScript RT Kit (Takara Bio, Inc., Dalian, China) was used to reversely transcribe RNA into cDNA in accordance with the manufacturer’s protocols. The SYBR Green PCR Kit (Takara Bio, Inc., Dalian, China) was used to conduct real-time fluorescence quantitative PCR with the ABI7500 Fluorescent Quantitative PCR Machine (Applied Biosystems, Lincoln Centre Drive Foster City, CA, USA). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the internal control. The primer sequences in this experiment were as follows: CCL5 (forward: 5′-CCAGCAGTCGTCTTTGTCAC-3′; reverse: 5′-CTCTGGGTTGGCACACACTT-3′), FASLG (forward: 5′-TGCCTTGGTAGGATTGGGC-3′; reverse: 5′-GCTGGTAGACTCTCGGAGTTC-3′), EOMES (forward: 5′-GCCATGCTTAGTGACACCGA-3′; reverse: 5′-GGACTGGAGGTAGTACCGC-3′), PDCD1 (forward: 5′-CCAGGATGGTTCTTAGACTCCC-3′; reverse: 5′-TTTAGCACGAAGCTCTCCGAT-3′), and GAPDH (forward: 5′-ACAACTTTGGTATCGTGGAAGG-3′ reverse: 5′-GCCATCACGCCACAGTTTC-3′). The mRNA expression levels were analyzed using the 2^-ΔΔCT method. The immunohistochemical data of the hub genes were obtained from The Human Protein Atlas (HPA) database (http://www.proteinatlas.org) (26). The method of the immunohistochemical analysis of the HPA could be found in the HPA website (https://www.proteinatlas.org/about/assays+annotation#ih).

Overall Survival Analysis of KIRP Patients in the TCGA Cohort

Survival-related hub genes were determined by the Kaplan–Meier (KM) survival analysis on the UALCAN (http://ualcan.path.uab.edu/analysis.html) (27). KIRP patients were categorized into the high expression (with gene expression levels above the upper quartile) group and the low/medium expression (with gene expression levels below the upper quartile) group. The log-rank p < 0.05 was considered to be statistically significant.

Correlation Analysis of Hub Genes and Immune Cells

TIMER (https://cistrome.shinyapps.io/timer/) was used for correlational analysis between the hub genes and immune cells in the TCGA cohort. The correlation analysis between the immune cells and the clinical prognosis was also performed on the TIMER. The R software was used to draw box plots of the proportion of CD8+ T cells and Treg cells in different clinical stages.

Identification of the Immune Subtypes by Consensus Clustering

We used the R package “ConsensusClusterPlus” to perform consensus clustering and the screening of the immune subtypes of KIRP in the TCGA cohort based on the 12 hub genes of the CD8+ T-cell-related module. In order to present the immune landscapes in different immune subtypes, we used the R package “ESTIMATE” to calculate the immune score, stromal score, and tumor purity of each tumor sample. The package “GSVA” was used to evaluate the single-sample gene-set enrichment analysis (ssGSEA) score based on 29 immune gene sets. The Kaplan–Meier (KM) survival curves of the immune subtypes were performed by using the R package “Survival.”

Prediction of Immunotherapy Responses in KIRP Patients

“TIDE” is a computational method that calculates the score of T-cell exclusion and T-cell dysfunction to predict the responses of immunotherapy (28). A lower TIDE score predicted a higher immunotherapy response. We classified KIRP patients into a true response group and a false response group by the median of TIDE scores in the TCGA cohort.

Acquisition and Processing of Validation Datasets

Three expression matrix files of KIRP, GSE2748 (n = 34), GSE7023 (n = 35), and GSE26574 (n = 34) were extracted from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). After merging the three expression matrix files, the batch effects were adjusted by the R package “sva.” The merged files were then used as the validation cohort.

Identification of Frequently Mutated Genes and Gene-Set Enrichment Analysis in Different Immunotherapy Response Groups

The somatic mutation files (maf format) of 273 KIRP patients of the TCGA cohort were downloaded from TCGA (https://cancergenome.nih.gov/). The R package “Maftools” was used to draw the waterfall plots in order to visualize the genes with mutation frequency >5% in the high and low/medium immunotherapy response groups. The r package “clusterProfiler” was used to perform GSEA analysis to evaluate the biological mechanisms between the high and low/medium immunotherapy response groups. We visualized the top five enrichment results of the high immunotherapy response group.

Statistical Analyses

Spearman correlation was applied to calculate the correlation coefficients between the hub genes. Wilcoxon test and Kruskal–Wallis test were applied to separately conduct the group comparisons of two groups and more than two groups. Overall survival curves were generated using the Kaplan–Meier method, and the group comparisons were performed with the log-rank test.

Results

Construction of a Weighted Co-Expression Network of Differentially Expressed Genes and Evaluation of Immune Infiltration

Under the condition of | log2 (fold-change) |> 1 and p < 0.05, 5,132 differentially expressed genes were screened (3,130 upregulated and 2,002 downregulated) in the TCGA cohort (Supplementary Figure 1A). The 5,132 differential genes in 289 tumor samples were used to construct a co-expression network. The optimal soft threshold β = 5 was used for the next calculation (Supplementary Figure 1B). The modules with ME of <0.25 in the cluster were merged, and 32 modules were finally obtained in the TCGA cohort (Supplementary Figures 1C, D). In the results of CIBERSORT, 137 tumor samples with a p < 0.05 were selected (Supplementary Figure 1E).

Correlation Between the Purple Module and T Cells

In the WGCNA analysis results, we found a high correlation between the purple module (169 genes) and CD8+ T cells (r = 0.67, p = 5e-38) (Figure 1A). We also plotted the relationship diagram of GS and MM of the purple module (Figure 1B). These results indicate that the purple module was significantly related to CD8+ T cells.

Figure 1

Analysis of the Function and Pathway of Genes Identified in the Purple Module

The genes of the purple module were related to a variety of immune functions. The first three significantly enriched GO biological processes in the purple module were about T-cell activation (Log10(P) = −26.47), response to interferon-gamma (Log10(P) = −16.90), and the cytokine-mediated signaling pathway (Log10(P) = −16.10). The first three KEGG pathways were antigen processing and presentation (Log10(P) = −11.27), cytokine–cytokine receptor interaction (Log10(P) = −11.19), and natural killer cell-mediated cytotoxicity (Log10(P) = −8.66) (Figure 1C, Supplementary Table 1). In addition, by using the MCODE algorithm of the Metascape platform, six closely linked network components were calculated and their functionally analyzed (Figure 1D, Supplementary Table 2).

Screening of Hub Genes of the Purple Module

A PPI network with 115 nodes and 661 edges was generated and visualized, and the genes in the top 30 scores were identified as hub genes of this network (Figure 2A). The nodes with a weight of >0.4 (81 nodes, 1,530 edges) between two nodes in the weighted co-expression network were visualized, and the genes in the top 30 scores were identified as hub genes of the weighted co-expression network (Figure 2B). In the purple module, genes that were simultaneously satisfied MM >0.9 and GS >0.6 were selected. Then, they were intersected to obtain 12 genes (GZMA, FASLG, CD2, PDCD1, CCL5, CD8A, CD3D, EOMES, NKG7, CD3E, CD8B, and CTLA4) (Figure 2C). We believed that these 12 genes were the true hub genes in the purple module and they were strongly correlated (Figure 2D).

Figure 2

Differential Expression of Hub Genes and Survival Analyses

As shown in Figure 3A, in the TCGA cohort, the expression of these 12 genes in tumor tissues was higher than that in normal tissues (p <0.05). The immunohistochemical data from HPA showed that EOMES and PDCD1 exhibited medium staining in the cytoplasm and cytomembrane of the tumor cells in the renal cancer. Normal kidney tissue staining of EOMES exhibited high staining in the cytoplasm and cytomembrane of cells in tubules, but as for cells in glomeruli, the staining was not detected. The staining of PDCD1 was not observed in normal kidney tissues (Figure 3B). In our samples, the mRNA levels of CCL5, FASLG, EOMES, and PDCD1 in tumor tissues were significantly higher than those in normal tissues (Figure 3C). Survival analysis revealed that the overexpression of eight genes was significantly associated with poor prognosis. As shown in Figure 3D, the p values of these genes are as follows: CCL5 (p = 0.015), CD2 (p = 0.025), CD8A (p = 0.004), CD8B (p = 0.039), EOMES (p = 0.045), FASLG (p = 0.015), GZMA (p = 0.012), and PDCD1 (p = 0.01).

Figure 3

Hub Genes Were Significantly Correlated With the Fraction of CD8+ T Cells

We validated the significant correlation between hub genes and the fraction of CD8+ T cells in the TCGA cohort by using the TIMER online platform (Figure 4A). The high infiltration status of B cells and CD8+ T cells was associated with poor prognosis in patients with KIRP (Figure 4B), and the number of CD8+ T cells and Treg cells increased with tumor progression (Figures 4C, D).

Figure 4

Immune Microenvironment Landscapes of the Five Immune Subtypes

Based on the 12 hub genes, the KIRP samples of the TCGA cohort were clustered into five subtypes (C1: 71 samples, C2: 52 samples, C3: 29 samples, C4: 87 samples, and C5: 44 samples) (Figures 5A–C). The prognostic outcomes of the five subtypes showed significant statistical differences (Figure 5D). The abundance of six types of immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and myeloid dendritic cells) was compared among the five subtypes (Figure 5E). The number of patients at different tumor stages in the immune subtypes is shown in Figure 5F. We also compared the ssGSEA score, immune score, stromal score, and tumor purity of the five immune subtypes (Figures 6A, B). Among the immune subtypes, C1 and C3 showed the highest abundance of the six immune cell types with the highest ssGSEA, immune, and stromal scores. C4 and C5 showed a moderate abundance of the six immune cell types along with moderate ssGSEA, immune, and stromal scores. Notably, C2 showed the lowest abundance of the six immune cell types with the lowest ssGSEA, immune, and stromal scores. On the other hand, tumor purity showed the opposite trend among these immune subtypes. Interestingly, in C1 and C3, which were rich in immune components, the percentage of patients at high tumor stages (III/IV) in C1 was 18.2%, and that in C3 was 44.8%.

Figure 5

Figure 6

Predictions of Immunotherapeutic Responses Among the Five Immune Subtypes

Among the five immune subtypes, C1 and C3 had the lowest TIDE and highest T-cell dysfunction scores. C4 and C5 had moderate TIDE and T-cell dysfunction scores. C2 showed the highest TIDE and lowest T-cell dysfunction scores (Figure 7A). The details of TIDE and T-cell dysfunction scores of each KIRP sample are presented in Supplementary Table 3. A sample with the TIDE score greater than 1.15 was defined as a false responder to immunotherapy, and a sample with the TIDE score less than 1.15 was defined as a true responder to immunotherapy. The percentage of true responders to immunotherapy decreased in the following order: C3 (96.6%) > C1 (81.7%) > C4 (44.8%) > C5 (29.5%) > C2 (7.7%) (Figure 7B). The expression of four immune checkpoint genes (LAG3, PD-L1, PD-1, and CTLA4) and two hub genes (CCL5 and FASLG) decreased in the following order: C3 > C1 > C4 > C5 > C2 (Figure 7C). Based on these findings, C1 and C3 were considered the high immunotherapy response group, C4 and C5 were considered the median immunotherapy response group, and C2 was considered the low immunotherapy response group.

Figure 7

Validation of the Immune Subtypes of KIRP

The GEO cohort was clustered into five subtypes based on the 12 hub genes to validate the immune subtypes of KIRP as identified in the TCGA cohort (Figures 8A–C). C1 and C3 had the highest abundance of immune cells, especially CD8+ T cells and myeloid dendritic cells. This observation was consistent with the results of the TCGA cohort (Figure 8D). In terms of immune-related gene-set enrichment score, immune score, stromal score, and tumor purity, the analysis results of the GEO cohort showed the same trend as that of the TCGA cohort (Figures 9A, B). In terms of the expression of two hub genes (CCL5 and FASLG) and four immune checkpoint genes (LAG3, PD-L1, PD-1, and CTLA4), the analysis results of the GEO cohort also exhibited the same trend as that of the TCGA cohort (Figure 9C). These comparable results between the TCGA and GEO cohorts indicated that the immune subtypes of KIRP identified in this study existed in KIRP.

Figure 8

Figure 9

Comparing the Immune Subtypes With Histological and Molecular Subtypes

Some studies have previously proposed the existence of two different histological subtypes and four different molecular subtypes in KIRP: type 1 with pale cytoplasm and small cells and type 2 with eosinophilic cytoplasm and large cells. Type 2 KIRP presented at a higher tumor stage more often than type 1 (29, 30). Molecular subtype KIRP.C1 was enriched in type 1 KIRP and associated with the MET mutation and gain of chromosome 7. Molecular subtype KIRP.C2a was enriched in type 2 KIRP and associated with DNA methylation cluster 2. Molecular subtype KIRP.C2b mainly consisted of type 2 and unclassified KIRP and was associated with SETD2 mutation and DNA methylation cluster 1. Molecular subtype KIRP.C2c-CIMP was enriched in type 2 KIRP and associated with hypermethylation of the CDKN2A promoter and mutation of FH (30). In the TCGA cohort, we compared the immune subtypes identified with the histological and molecular subtypes. In terms of histological subtypes, C3 was mainly composed of type 2 and C2 was mainly composed of type 1. The composition ratio of the histological subtypes was comparable among C1, C4, and C5 (Figure 10A). C3 was mainly presented in type 2 (Figure 10B). With regard to molecular subtypes, KIRP.C2a accounted for a large proportion of C3, whereas KIRP.C1 was dominant in C1, C2, C4, and C5 (Figure 10C). C3 was mainly presented in KIRP.C2a and KIRP.C2c-CIMP. C2 mainly appeared in KIRP.C1 and KIPR.C2b (Figure 10D).

Figure 10

The Mutation Frequency of TTN Was Lower in the High Immunotherapy Response Group Than in the Low/Medium Immunotherapy Response Group

In the TCGA cohort, the genes with a mutation frequency greater than 5% were different in the high and low/medium immunotherapy response groups. MUC16, KMT2C, MET, TTN, and MUC4 were common between the two groups. OBSCN, ARID1A, FAT1, USH2A, CENPF, HELZ2, and WDFY3 were specific to the high immunotherapy response group. SETD2, KIAA1109, CUBN, KMT2D, MACF1, PCLO, DNAH8, KDM6A, LRP2, PBRM1, PCF11, PKHD1, and SYNE1 were specific to the low/medium immunotherapy response group. The difference in the mutation frequency of TTN was most significant between these two groups. The high response group was 7% and the low/medium group was 17% (Figures 11A, B). GSEA performed on these two groups revealed that multiple immune-related functions and pathways such as eosinophil chemotaxis, eosinophil migration, IgG binding, interleukin-12 secretion, MHC class II protein complex, allograft rejection, graft-versus-host disease, the intestinal immune network for IgA production, systemic lupus erythematosus, and type I diabetes mellitus were significantly enriched in the high immunotherapy response group (Figures 11C, D).

Figure 11

Discussion

Tumor immunotherapy targeting immune checkpoints has shown encouraging therapeutic effects. Because of the differences in individual responses of patients, screening of the beneficiary population for specific tumor subtypes and identifying new targets and biomarkers to develop single-therapy drugs and combined therapy programs still need to be intensively performed. Based on the computational biology method and the in-depth analysis of the multidimensional data of 289 patients with KIRP in the TCGA cohort and 103 patients with KIRP in the GEO cohort, our research shows the following. 1) There is a certain number of potentially functional immune cells in the tumors of patients with KIRP but the tumor-infiltrating CD8+ T cells are exhausted, causing them to fail to exert their antitumor effects. 2) The C1 and C3 immune subtypes that we identified may have the highest clinical benefit in reversing CD8+ T-cell exhaustion. Our comprehensive evaluation of the tumor immune microenvironment of patients with KIRP suggests that these patients may be suitable for immunotherapy that reverses CD8+ T-cell exhaustion.

Based on the comprehensive evaluation of the tumor immune microenvironment of patients with KIRP, we identified 12 hub genes that were closely related to this tumor microenvironment. These 12 hub genes were highly expressed in these patients and showed a strong correlation with CD8+ T cells. The roles of these 12 hub genes in tumors can be divided into two major groups: those with immunostimulatory functions (CD2, CD8A, CD8B, CD3D, CD3E, EOMES, and GZMA) and those with immunosuppressive functions (PDCD1, CTLA4, FASLG, and CCL5). Some of these genes are cancer suppressor genes. CD2, CD8A, CD8B, CD3D, and CD3E are membrane proteins on the surface of T cells and play an important role in the immune recognition and activation of T cells (31–33). EOMES has an important effect on the development of NK cells and the differentiation of CD8+ T cells (34, 35). GZMA, a cytotoxic protein secreted by NK cells and cytotoxic T cells, can induce caspase-dependent cell death (36). At present, the role of NKG7 in cancer has not been elucidated. In addition, some of these genes are involved in the formation of an immunosuppressive microenvironment in the tumor. CTLA4 is an inhibitory receptor found on the surface of T cells. It can downregulate the activity of T cells after binding to CD80 and CD86 on antigen-presenting cells (37, 38) (Figure 12). PDCD1 belongs to the cell-surface immunoglobulin superfamily. After interacting with PD-L1, the proliferation and activity of tumor-specific T cells are inhibited (39) (Figure 12). FASLG is a member of the tumor necrosis factor superfamily. Soluble FASLG can kill effector immune cells and cause immune tolerance (40, 41) (Figure 9). In adoptive cell immunotherapy, when genetically engineered Fas receptor-negative T cells are used for cancer treatment, Fas receptor-negative T cells show longer durability and stronger antitumor immunity both in the periphery blood and in tumors of tumor-bearing animals (42). The chemokine CCL5 can recruit Treg cells to enhance immune tolerance (43) (Figure 12). In our study, CCL5 and Treg cells also showed a strong correlation (Supplementary Figure 2). Additionally, CCL5 produced by immune cells can promote tumor growth and proliferation by regulating macrophage production of metalloproteinases or inducing epithelial–mesenchymal transformation in tumor cells (44, 45). There is evidence that drugs blocking the CCR5–CCL5 axis or reducing CCL5 production can be used in the clinical setting (46). It is worth noting that in addition to the common immune checkpoint targets such as PD-1, we found that CCL5 and FASLG may play an important role in the formation of the KIRP tumor immunosuppressive microenvironment; however, studies on the role of CCL5 and FASLG expression in KIRP have not been performed. Therefore, whether CCL5 and FASLG and their receptors can be used as new targets for KIRP immunotherapy remains to be further studied.

Figure 12

The main finding of this study was to identify that patients with KIRP are suitable for receiving immunotherapy that can reverse CD8+ T-cell exhaustion. In our study, we found a strong co-expression relationship between some genes with immunosuppressive effects and those with immune-promoting effects in the purple module, indicating that the greater the activation potential and antitumor potential of CD8+ T cells, the greater the possibility of inhibition. In addition, our results showed that in the tumor of patients with KIRP, the numbers of Treg cells with immunosuppressive effects and CD8+ T cells continue to increase with tumor progression. This finding is consistent with the observation that the number of CD8+ T cells in RCC is positively correlated with the tumor grade but because of the immunosuppressive environment in renal cell carcinoma, TILs with lower proliferative activity cannot exert their function effectively (47, 48). Survival analysis showed that CD8+ T cells were significantly associated with poor prognosis of patients with KIRP. Based on the above findings, we can infer that KIRP belongs to the “hot” tumor category. Although there is a certain number of CD8+ T cells around KIRP tumor cells, because of an immunosuppressive environment and the expression of immunosuppressive receptors by CD8+ T cells, CD8+ T cells are exhausted, resulting in the inability of CD8+ T cells to exert their antitumor effects. Therefore, reversing exhausted CD8+ T cells and restoring their antitumor ability could be a reasonable choice in KIRP immunotherapy.

To find potentially suitable patients who could benefit the most from the treatment of reversing CD8+ T-cell exhaustion, we clustered a TCGA cohort of KIRP into five immune subtypes based on 12 hub genes associated with CD8+ T cells and validated these immune subtypes in the GEO cohort. C1 and C3 had the highest immune component, T-cell dysfunction score, and percentages of true responders to immunotherapy (C1: 81.7%, C3: 96.6%). C2 had the lowest immune component and T-cell dysfunction score. In terms of percentages of true esponders to immunotherapy, C2 had the lowest one (7.7%). C4 and C5 had moderate immune components, T-cell dysfunction scores, and percentages of true responders to immunotherapy (C4: 44.8%, C5: 29.5%). Therefore, among the five immune subtypes, high immunotherapy response subtypes C1 and C3 could obtain the best efficacy of reversing CD8+ T-cell exhaustion therapy. The low immunotherapy response subtype C2 was not suitable for reversing CD8+ T-cell exhaustion therapy. It is worth noting that although C1 and C3 are both subtypes of high immunotherapy response, the number of patients at higher tumor stages (III/IV) in C1 was 18.2% and that in C3 was 44.8%. This finding suggests that immunotherapeutic interventions for early-stage patients of type C1 may be effective in preventing patients from transitioning to late stages.

In conclusion, based on an in-depth analysis of the multiomic and multidimensional data of KIRP from the largest sample available currently, we found that patients with KIRP whose immune microenvironment exhibits “hot” tumor characteristics are suitable for receiving therapy that causes the reversal of CD8+ T-cell exhaustion, and their C1 and C3 immune subtypes may achieve the best therapeutic effect. Genes such as CCL5 and FASLG may play a crucial role in the formation of the KIRP immunosuppressive microenvironment. Given the limited high-quality clinical data of KIRP, further experimental and clinical studies are required to confirm the above findings and explore the corresponding immunotherapy regimens.

Funding

This work was supported by the national natural science foundation of China (grants 81672523, 81472404, 81472403, 81272834, and 31000572); the 2018 support plan for innovative talents in colleges and universities of Liaoning province; the 2018 “million talents project” funded project of Liaoning province (grant 33013); and the 2019 Key R & D projects of Shenyang (grant 19-112-4-102).

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  LjungbergBBensalahKCanfieldSDabestaniSHofmannFHoraMet al. EAU Guidelines on Renal Cell Carcinoma: 2014 Update. Eur Urol (2015) 67(5):913–24. doi: 10.1016/j.eururo.2015.01.005

  • CrossRef
  • Google Scholar

• 2

  CourthodGTucciMDi MaioMScagliottiGV. Papillary Renal Cell Carcinoma: A Review of the Current Therapeutic Landscape. Crit Rev Oncol Hematol (2015) 96(1):100–12. doi: 10.1016/j.critrevonc.2015.05.008

  • CrossRef
  • Google Scholar

• 3

  BarataPCRiniBI. Treatment of Renal Cell Carcinoma: Current Status and Future Directions. CA Cancer J Clin (2017) 67(6):507–24. doi: 10.3322/caac.21411

  • CrossRef
  • Google Scholar

• 4

  GarciaJARiniBI. Recent Progress in the Management of Advanced Renal Cell Carcinoma. CA Cancer J Clin (2007) 57(2):112–25. doi: 10.3322/canjclin.57.2.112

  • CrossRef
  • Google Scholar

• 5

  GrivennikovSIGretenFRKarinM. Immunity, Inflammation, and Cancer. Cell (2010) 140(6):883–99. doi: 10.1016/j.cell.2010.01.025

  • CrossRef
  • Google Scholar

• 6

  VeselyMDKershawMHSchreiberRDSmythMJ. Natural Innate and Adaptive Immunity to Cancer. Annu Rev Immunol (2011) 29:235–71. doi: 10.1146/annurev-immunol-031210-101324

  • CrossRef
  • Google Scholar

• 7

  AlbersAEFerrisRLKimGGChikamatsuKDeLeoABWhitesideTL. Immune Responses to P53 in Patients With Cancer: Enrichment in Tetramer+ P53 Peptide-Specific T Cells and Regulatory T Cells at Tumor Sites. Cancer Immunol Immunother (2005) 54(11):1072–81. doi: 10.1007/s00262-005-0670-9

  • CrossRef
  • Google Scholar

• 8

  KimJWWieckowskiETaylorDDReichertTEWatkinsSWhitesideTL. Fas Ligand-Positive Membranous Vesicles Isolated From Sera of Patients With Oral Cancer Induce Apoptosis of Activated T Lymphocytes. Clin Cancer Res (2005) 11(3):1010–20.

  • Google Scholar

• 9

  FifeBTBluestoneJA. Control of Peripheral T-Cell Tolerance and Autoimmunity via the CTLA-4 and PD-1 Pathways. Immunol Rev (2008) 224:166–82. doi: 10.1111/j.1600-065X.2008.00662.x

  • CrossRef
  • Google Scholar

• 10

  FourcadeJSunZBenallaouaMGuillaumePLuescherIFSanderCet al. Upregulation of Tim-3 and PD-1 Expression Is Associated With Tumor Antigen-Specific CD8+ T Cell Dysfunction in Melanoma Patients. J Exp Med (2010) 207(10):2175–86. doi: 10.1084/jem.20100637

  • CrossRef
  • Google Scholar

• 11

  AnsellSMLesokhinAMBorrelloIHalwaniAScottECGutierrezMet al. PD-1 Blockade With Nivolumab in Relapsed or Refractory Hodgkin's Lymphoma. N Engl J Med (2015) 372(4):311–9. doi: 10.1056/NEJMoa1411087

  • CrossRef
  • Google Scholar

• 12

  BorghaeiHPaz-AresLHornLSpigelDRSteinsMReadyNEet al. Nivolumab Versus Docetaxel in Advanced Nonsquamous Non-Small-Cell Lung Cancer. N Engl J Med (2015) 373(17):1627–39. doi: 10.1056/NEJMoa1507643

  • CrossRef
  • Google Scholar

• 13

  FerrisRLBlumenscheinGJr.FayetteJGuigayJColevasADLicitraLet al. Nivolumab for Recurrent Squamous-Cell Carcinoma of the Head and Neck. N Engl J Med (2016) 375(19):1856–67. doi: 10.1056/NEJMoa1602252

  • CrossRef
  • Google Scholar

• 14

  HodiFSO'DaySJMcDermottDFWeberRWSosmanJAHaanenJBet al. Improved Survival With Ipilimumab in Patients With Metastatic Melanoma. N Engl J Med (2010) 363(8):711–23. doi: 10.1056/NEJMoa1003466

  • CrossRef
  • Google Scholar

• 15

  MotzerRJEscudierBMcDermottDFGeorgeSHammersHJSrinivasSet al. Nivolumab Versus Everolimus in Advanced Renal-Cell Carcinoma. N Engl J Med (2015) 373(19):1803–13. doi: 10.1056/NEJMoa1510665

  • CrossRef
  • Google Scholar

• 16

  Adrianzen HerreraDAFleisigSBGartrellBA. Impressive and Durable Response to Nivolumab in a Patient With Metastatic Type 2 Papillary Renal Cell Carcinoma: On-Label But Without Evidence. Invest New Drugs (2017) 35(5):665–8. doi: 10.1007/s10637-017-0469-5

  • CrossRef
  • Google Scholar

• 17

  NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

  • CrossRef
  • Google Scholar

• 18

  LiTFanJWangBTraughNChenQLiuJSet al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77(21):e108–10. doi: 10.1158/0008-5472.CAN-17-0307

  • CrossRef
  • Google Scholar

• 19

  LiuZP. Identifying Network-Based Biomarkers of Complex Diseases From High-Throughput Data. Biomark Med (2016) 10(6):633–50. doi: 10.2217/bmm-2015-0035

  • CrossRef
  • Google Scholar

• 20

  AllenJDXieYChenMGirardLXiaoG. Comparing Statistical Methods for Constructing Large Scale Gene Networks. PloS One (2012) 7(1):e29348. doi: 10.1371/journal.pone.0029348

  • CrossRef
  • Google Scholar

• 21

  LangfelderPHorvathS. WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

  • CrossRef
  • Google Scholar

• 22

  LangfelderPZhangBHorvathS. Defining Clusters From a Hierarchical Cluster Tree: The Dynamic Tree Cut Package for R. Bioinformatics (2008) 24(5):719–20. doi: 10.1093/bioinformatics/btm563

  • CrossRef
  • Google Scholar

• 23

  SzklarczykDMorrisJHCookHKuhnMWyderSSimonovicMet al. The STRING Database in 2017: Quality-Controlled Protein-Protein Association Networks, Made Broadly Accessible. Nucleic Acids Res (2017) 45(D1):D362–8. doi: 10.1093/nar/gkw937

  • CrossRef
  • Google Scholar

• 24

  ChinCHChenSHWuHHHoCWKoMTLinCY. Cytohubba: Identifying Hub Objects and Sub-Networks From Complex Interactome. BMC Syst Biol (2014) 8(Suppl 4):S11. doi: 10.1186/1752-0509-8-S4-S11

  • CrossRef
  • Google Scholar

• 25

  ChenHBoutrosPC. VennDiagram: A Package for the Generation of Highly-Customizable Venn and Euler Diagrams in R. BMC Bioinf (2011) 12:35. doi: 10.1186/1471-2105-12-35

  • CrossRef
  • Google Scholar

• 26

  UhlenMFagerbergLHallstromBMLindskogCOksvoldPMardinogluAet al. Proteomics. Tissue-Based Map of the Human Proteome. Science (2015) 347(6220):1260419. doi: 10.1126/science.1260419

  • CrossRef
  • Google Scholar

• 27

  ChandrashekarDSBashelBBalasubramanyaSAHCreightonCJPonce-RodriguezIChakravarthiBet al. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia (2017) 19(8):649–58. doi: 10.1016/j.neo.2017.05.002

  • CrossRef
  • Google Scholar

• 28

  JiangPGuSPanDFuJSahuAHuXet al. Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

  • CrossRef
  • Google Scholar

• 29

  DelahuntBEbleJN. Papillary Renal Cell Carcinoma: A Clinicopathologic and Immunohistochemical Study of 105 Tumors. Mod Pathol (1997) 10(6):537–44.

  • Google Scholar

• 30

  Cancer Genome Atlas Research NLinehanWMSpellmanPTRickettsCJCreightonCJFeiSSet al. Comprehensive Molecular Characterization of Papillary Renal-Cell Carcinoma. N Engl J Med (2016) 374(2):135–45. doi: 10.1056/NEJMoa1505917

  • CrossRef
  • Google Scholar

• 31

  BiererBEPetersonAGorgaJCHerrmannSHBurakoffSJ. Synergistic T Cell Activation via the Physiological Ligands for CD2 and the T Cell Receptor. J Exp Med (1988) 168(3):1145–56. doi: 10.1084/jem.168.3.1145

  • CrossRef
  • Google Scholar

• 32

  LeongSRSukumaranSHristopoulosMTotpalKStaintonSLuEet al. An Anti-CD3/anti-CLL-1 Bispecific Antibody for the Treatment of Acute Myeloid Leukemia. Blood (2017) 129(5):609–18. doi: 10.1182/blood-2016-08-735365

  • CrossRef
  • Google Scholar

• 33

  RosensteinYRatnofskySBurakoffSJHerrmannSH. Direct Evidence for Binding of CD8 to HLA Class I Antigens. J Exp Med (1989) 169(1):149–60. doi: 10.1084/jem.169.1.149

  • CrossRef
  • Google Scholar

• 34

  GordonSMChaixJRuppLJWuJMaderaSSunJCet al. The Transcription Factors T-Bet and Eomes Control Key Checkpoints of Natural Killer Cell Maturation. Immunity (2012) 36(1):55–67. doi: 10.1016/j.immuni.2011.11.016

  • CrossRef
  • Google Scholar

• 35

  PearceELMullenACMartinsGAKrawczykCMHutchinsASZediakVPet al. Control of Effector CD8+ T Cell Function by the Transcription Factor Eomesodermin. Science (2003) 302(5647):1041–3. doi: 10.1126/science.1090148

  • CrossRef
  • Google Scholar

• 36

  CullenSPBrunetMMartinSJ. Granzymes in Cancer and Immunity. Cell Death Differ (2010) 17(4):616–23. doi: 10.1038/cdd.2009.206

  • CrossRef
  • Google Scholar

• 37

  KormanAJPeggsKSAllisonJP. Checkpoint Blockade in Cancer Immunotherapy. Adv Immunol (2006) 90:297–339. doi: 10.1016/S0065-2776(06)90008-X

  • CrossRef
  • Google Scholar

• 38

  MeleroIHervas-StubbsSGlennieMPardollDMChenL. Immunostimulatory Monoclonal Antibodies for Cancer Therapy. Nat Rev Cancer (2007) 7(2):95–106. doi: 10.1038/nrc2051

  • CrossRef
  • Google Scholar

• 39

  ButteMJKeirMEPhamduyTBSharpeAHFreemanGJ. Programmed Death-1 Ligand 1 Interacts Specifically With the B7-1 Costimulatory Molecule to Inhibit T Cell Responses. Immunity (2007) 27(1):111–22. doi: 10.1016/j.immuni.2007.05.016

  • CrossRef
  • Google Scholar

• 40

  AndreolaGRivoltiniLCastelliCHuberVPeregoPDehoPet al. Induction of Lymphocyte Apoptosis by Tumor Cell Secretion of FasL-Bearing Microvesicles. J Exp Med (2002) 195(10):1303–16. doi: 10.1084/jem.20011624

  • CrossRef
  • Google Scholar

• 41

  KimREmiMTanabeKUchidaYTogeT. The Role of Fas Ligand and Transforming Growth Factor Beta in Tumor Progression: Molecular Mechanisms of Immune Privilege via Fas-Mediated Apoptosis and Potential Targets for Cancer Therapy. Cancer (2004) 100(11):2281–91. doi: 10.1002/cncr.20270

  • CrossRef
  • Google Scholar

• 42

  YamamotoTNLeePHVodnalaSKGurusamyDKishtonRJYuZet al. T Cells Genetically Engineered to Overcome Death Signaling Enhance Adoptive Cancer Immunotherapy. J Clin Invest (2019) 129(4):1551–65. doi: 10.1172/JCI121491

  • CrossRef
  • Google Scholar

• 43

  YouYLiYLiMLeiMWuMQuYet al. Ovarian Cancer Stem Cells Promote Tumour Immune Privilege and Invasion via CCL5 and Regulatory T Cells. Clin Exp Immunol (2018) 191(1):60–73. doi: 10.1111/cei.13044

  • CrossRef
  • Google Scholar

• 44

  BarashiNWeissIDWaldOWaldHBeiderKAbrahamMet al. Inflammation-Induced Hepatocellular Carcinoma Is Dependent on CCR5 in Mice. Hepatology (2013) 58(3):1021–30. doi: 10.1002/hep.26403

  • CrossRef
  • Google Scholar

• 45

  LiouGYDopplerHNecelaBKrishnaMCrawfordHCRaimondoMet al. Macrophage-Secreted Cytokines Drive Pancreatic Acinar-to-Ductal Metaplasia Through NF-KappaB and MMPs. J Cell Biol (2013) 202(3):563–77. doi: 10.1083/jcb.201301001

  • CrossRef
  • Google Scholar

• 46

  AldinucciDColombattiA. The Inflammatory Chemokine CCL5 and Cancer Progression. Mediat Inflamm (2014) 2014:292376. doi: 10.1155/2014/292376

  • CrossRef
  • Google Scholar

• 47

  LahnMFischPKohlerGKunzmannRHentrichIJesuiterHet al. Pro-Inflammatory and T Cell Inhibitory Cytokines Are Secreted at High Levels in Tumor Cell Cultures of Human Renal Cell Carcinoma. Eur Urol (1999) 35(1):70–80. doi: 10.1159/000019821

  • CrossRef
  • Google Scholar

• 48

  NakanoOSatoMNaitoYSuzukiKOrikasaSAizawaMet al. Proliferative Activity of Intratumoral CD8(+) T-Lymphocytes as a Prognostic Factor in Human Renal Cell Carcinoma: Clinicopathologic Demonstration of Antitumor Immunity. Cancer Res (2001) 61(13):5132–6.

  • Google Scholar