Integrative Analysis of Immune-Related Genes in the Tumor Microenvironment of Renal Clear Cell Carcinoma and Renal Papillary Cell Carcinoma

Kidney cancer encompasses a range of primary cancers, such as clear cell renal cell carcinoma (ccRCC) and papillary renal cell carcinoma (pRCC). Our knowledge about the tumor microenvironment (TME) of kidney cancer is still limited. Therefore, we comprehensively assessed the TME of kidney cancers (including ccRCC and pRCC) using the ESTIAMTE, and CIBERSORT algorithms, and conducted distinct functional and correlation analyses with data from The Cancer Genome Atlas (TCGA), International Cancer Genome Consortium (ICGC), Gene Expression Omnibus (GEO), Connectivity map and CellMiner database. Next, we identified two immune-related hub genes, IGLL5 and IL2RA, which play essential roles in the TME as well as on survival in ccRCC and pRCC. Furthermore, ccRCC and pRCC samples from our medical center were collected to verify the clinical application value of these two immune-related genes. A specific enrichment analysis of immune cells related to IGLL5 and IL2RA was also conducted in two types of renal cell cancer. Based on selected genes, we predicted the drug response and uncovered novel drug candidate for RCC treatment. Considering the unfavorable outcomes of kidney cancer and emerging interest in TME-targeted treatments, our results may offer insights into immune-related molecular mechanisms and possible targets to control the kidney cancer.


INTRODUCTION
Kidney or renal cancer is the sixth most common malignant cancer in men and the ninth most common in women. In 2020, it was estimated that 76,080 patients would be diagnosed with kidney cancer and 13,780 patients would die from the disease in the United States (Siegel et al., 2021). Renal cell carcinoma (RCC) is the most common form and accounts for 85% of all histological types, with ccRCC comprising 70-80% and pRCC 15-20% of all RCC (Linehan et al., 2016;Barata and Rini, 2017;Vuong et al., 2019). Although pRCC has a better prognosis than ccRCC, the overall prognosis for ccRCC and pRCC remains limited (Steffens et al., 2012). Therefore, effective therapy and accurate biomarkers need to be determined. TME, as an integral part of tumors, is a cellular environment consisting of tumor cells and other non-malignant cells, including surrounding immune cells, lymphocytes, fibroblasts, stromal cells, and blood vessels (Wu and Dai, 2017). TME acts as a complex ecosystem, which could support tumor growth as well as metastasis while attenuating immunosurveillance. A wealth of new information has emerged which reveals how the functionality of TME determines its integral and indispensable role in various cancers (Hanahan and Coussens, 2012;Pitt et al., 2016). Within the TME, various types of immune and nonimmune cell are found. Some nontumor cells, such as stromal cells, fibroblasts, may accelerate cancer cell proliferation and stimulate cancer progression (Hanahan and Coussens, 2012). However, although several immune cells, such as T cells as well as B cells also exist in TME, the exact contribution of these immune cells to the prognosis of patients is not clear (Grivennikov et al., 2010). With a variety of factors that diverse cells secrete, these lead to a chronic inflammatory, immunosuppressive environment (Osipov et al., 2019). Within that environment, cancer cells are able to adapt and grow with less possibility of detection and eradication by immunosurveillance.
Currently, RCC is considered an immunogenic and vascularized tumor. Several studies have found that immune cells could infiltrate into the TME of RCC, but these immune cells block the effective anti-tumor responses (Díaz-Montero et al., 2020). Due to the immunosuppressed state of RCC, attention has focused on immune checkpoint inhibitors. However, a large number of patients with RCC do not benefit from these immune-based treatments, and therapies are related with drug resistance because of heterogeneous and adaptive TME (Ravaud et al., 2016). In order to eradicate tumor cells, effector immune cells must be relieved from the immunosuppressive networks that constitute the TME. Moreover, with the improved understanding of interactions between immune cells, and between cancer cells and immune cells, the optimizing of our approach to treating cancers will significantly progress (Joyce and Pollard, 2009). Therefore, it is indispensable to elucidate the cross-talk process between various immune cells and the TME.
In this study, we identified the immune-related genes in the TME of ccRCC as well as pRCC, and explored their potential mechanisms, biological functions, and predict drug response as well as potential drugs for ccRCC and pRCC. These immunegenes could serve as reliable prognostic biomarkers as well as therapeutic targets.

Data Collection
The gene expression data and clinical information of ccRCC (539 ccRCC and 72 normal samples) and pRCC (289 pRCC and 32 normal samples) were downloaded directly from TCGA (https:// portal.gdc.cancer.gov/). The gene expression data and clinical information of 91 and 61 RCC samples were obtained from ICGC (https://dcc.icgc.org/) and GEO (including GSE40912 and GSE2748) (https://www.ncbi.nlm.nih.gov/geo/), separately. We eliminated samples in which clinical information was lost or unknown.

Identification and Analysis of Differentially Expressed Genes
Differentially expressed genes (DEGs) were identified using the "limma" package and visualized using the "pheatmap" package. The |logFC| > 1, p-value < 0.05, and false discovery rate (FDR) < 0.05 were used to screen for DEGs. The correlation between immune-related genes and clinical factors was performed using the "ggplot2" and "corrplot" package.

Functional Analysis
Protein-protein interaction (PPI) network construction was performed using the STRING online database (http://string-db. org/), and the network was visualized using Cytoscape (Szklarczyk et al., 2019). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses was conducted by "clusterProfiler" package. Gene set enrichment analysis (GSEA) was performed by the "gsva" R package. The Connectivity Map (Cmap) and CellMiner database were uesd to predict small molecular compounds and individual response to drugs.

Tumor Microenvironment Analysis
We used ESTIMATE algorithm (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) algorithm to measure stromal and immune scores in ccRCC and pRCC with data from TCGA database (https:// portal.gdc.cancer.gov/) (Yoshihara et al., 2013). The prognostic-related genes were selected after performing the Univariate Cox analysis with p < 0.1. Tumor-infiltrating immune cells were calculated by the CIBERSORT algorithm. Venn diagrams as well as violin plots were visualized with the 'VennDiagram' and 'vioplot' packages.

Immunohistochemistry
Ten pairs of ccRCC, pRCC and adjacent normal tissues were collected from Shandong provincial Hospital affiliated to Shandong First University from April 2021 to June 2021. The study was approved by the Ethics Committee of Shandong Provincial Hospital (Approval no. SWYX: NO. 2021-118). IHC was performed according to previously published method . All samples were incubated with rabbit polyclonal anti-IL2RA (ab245687) and anti-IGLL5 (NBP2-14574) antibodies overnight at 4°C and then washed. Two pathologists independently assessed the IHC slides.

Statistical Analysis
Mann-Whitney test was utilized to measure gene expression between tumor and non-tumor tissues. The Kaplan-Meier curves with a log-rank test were adopted in the survival analysis. Statistical analysis in this study was performed with the R statistical package (R version 4.0.1).

Correlation of ESTIMATE Scores in ccRCC and pRCC
We first determined the association between ESTIMATE scores (including immune scores and stromal scores), prognosis as well as clinical characteristics in ccRCC and pRCC. As shown in Figures 1A,B, ccRCC patients with high immune scores had poor prognosis, but stromal scores did not have apparently association with survival. However, in pRCC patients, a decrease in overall survival was correlated with high immune scores and stromal scores ( Figures 1C,D). In addition, for ccRCC, immune scores were positively associated with tumor grade, tumor stage, T stage, and M stage. Patients with advanced T stage had relatively higher stromal scores ( Figures 1E-K). For pRCC, higher immune scores were associated with T stage, while patients with advanced T stage and N stage have high stromal scores ( Figures 1L-N).

Identification of Differentially Expressed Genes
Based on immune scores and stromal scores, we determined the DEGs by analyzing the expression data of all ccRCC and pRCC samples. In ccRCC, 1066 genes were upregulated and 204 genes were downregulated in samples with high immune score. Similarly, we also identified 515 genes that were highly expressed and 221 genes with low expression in ccRCC samples with high stromal score (logFC >1; p < 0.05; Figures  2A,B). In pRCC, 1155 genes had high expression levels and 174 genes had low expression in patients with high immune scores, and 1395 genes were upregulated as well as 141 downregulated in patients with high stromal scores ( Figures 2C,D). Using Venn plot, we determined 218 genes with high expression and 68 genes with low expression in ccRCC patients with high scores (Figures 2E,F). pRCC patients with high scores presented 961 genes that were upregulated and 38 genes that were downregulated, respectively ( Figures 2G,H).
On the basis of these selected DEGs, we performed GO and KEGG enrichment analyze in two types of kidney cancers. The GO analysis showed that DEGs identified in ccRCC and pRCC cases were primarily enriched in immune-related functions, such as B cell-mediated immunity and lymphocyte-mediated immunity ( Figures 2J-L). We also listed the top 10 KEGG pathways for both ccRCC and pRCC. As shown in Figures 2I-K, various common immune-related signaling pathways were enriched. Furthermore, DEGs from pRCC samples were also enriched in several pathways related to T-cell receptor and B-cell receptor signaling pathways ( Figure 2K).

PPI Network and Prognostic Analysis of DEGs
Then these identified DEGs were used to construct the PPI networks. We ranked the DEGs genes by their nodes (Figures . Moreover, we further investigated prognostic genes by performing the univariate Cox analysis among all DEGs in ccRCC and pRCC, respectively. This revealed that 69 genes were related to ccRCC prognosis, and 110 genes were associated with the prognosis of pRCC patients (p < 0.1, Figures 3C,D).
Then, we interrogated the top 15 prognosis-related hub genes in ccRCC and pRCC utilizing Venn algorithm ( Figures 3E,F). Among 30 selected prognosis-related hub genes in ccRCC and pRCC, two common genes, IGLL5 and IL2RA, were identified, which could potentially play essential role in TME as well as survival in ccRCC and pRCC. Correlation and GSEA Enrichment Analyses of IGLL5 and IL2RA To further assess the expression profile and clinical values of IGLL5 and IL2RA, we extracted the corresponding expression data, as well as clinical information. Both IGLL5 and IL2RA showed high expression in ccRCC samples (p < 0.05, Figures  4A,B). Furthermore, a decrease in overall survival was associated with high expression of IGLL5 and IL2RA (p < 0.05, Figures  4C,D). IGLL5 and IL2RA expression were also positively associated with the stage of ccRCC patients, including high tumor stage, tumor grade, T stage, and more lymph node metastases (p < 0.05, Figures 4E-K).
The same approach was used to analyze two genes in pRCC patients. The differential expression analysis confirmed that IL2RA expressed highly in tumor as well (p < 0.05, Figure 4M), but significant differences in expression of IGLL5 between normal and tumor patients was not observed (p > 0.05, Figure 4M). Based on these results, we chose IL2RA for further analysis. In pRCC patients, poor survival was associated with high expression of IL2RA (p < 0.05, Figure 4N). Moreover, high IL2RA expression was associated with various clinical characteristics, including T and N stage indicative of disease progression. (p < 0.05, Figures 4O,P).
Next, we explored the potential regulation pathways through GSEA enrichment analysis in ccRCC and pRCC samples. Several immune-related pathways, including B cell and T cell receptor signaling pathways, chemokine signaling, and natural killer cell mediated cytotoxicity, which were closely correlated with IGLL5 and IL2RA expression (NES >1, p < 0.05, FRD <0.25, Supplementary Figures S2A-C).

Composition of Immune Cells and Correlation Analysis
To further verify the essential role of major immune cells in ccRCC and pRCC, we explored immune cell infiltration and correlation between two genes and major immune cells. The profile of the principle immune cells showed that adaptive immune cells (such as native B cells, memory B cells, native and memory CD4 + cell, CD8 + T cells) accounted for a substantial part of the immune cells in ccRCC as well as pRCC ( Figures 5A,B). Next, we investigated differences in the fraction of major types of immune cells according to expression levels of IGLL5 and IL2RA. For ccRCC, high IGLL5 expression was directly proportional to most adaptive immune cells (including native B cell, CD4 + memory T cells, CD8 + T cells, T follicular helper cells, and regulatory T cells) and was inversely proportion to several innate immune cells (such as NK cell, monocytes, M2 macrophages, activated dendritic cells, and mast cells) (p < 0.05, Figure 5C). Moreover, the difference analysis between high and low expression level of IL2RA in ccRCC illustrated that the high IL2RA expression set presented a significantly larger fraction of native B cells, CD4 + memory T cells, and M2 macrophages, than CD8 + T cells, NK cells, and neutrophils (p < 0.05, Figure 5D). The same approach was used to assess IL2RA in pRCC, which revealed that high IL2RA expression set had higher fraction of native B cells, neutrophils, and a lower proportion of memory B cells, activated CD4 + T memory cells, activated NK cells and resting mast cells (p < 0.05, Figure 5E).
We also detected the correlation of IGLL5 and IL2RA with immune cell abundance. In ccRCC, the correlation analyses showed that IGLL5 and most adaptive immune cells had a clear linear positive correlation and had a negative correlation with the infiltration of native immune cells (p < 0.05, Supplementary Figure S3). For IL2RA, there was a positive relationship with native B cells, activated CD4 + memory T cells, and M0 and M2 macrophages. However, an inverse correlation was observed between IL2RA and CD8 + T cells, T follicular helper cells, and NK cells (p < 0.05, Supplementary Figure S4). In pRCC, IL2RA showed a positive correlativity with activated CD4 + T memory cells, native B cells, M1 macrophages and neutrophils but had negative relationship with resting CD4 + T memory cells, memory B cells, and mast cells (p < 0.05, Supplementary Figure S5). Finally, on the basis of above analyses, we further investigated the intersected immune cells associated with both IGLL5 and IL2RA in ccRCC and pRCC (

Prognostic Value of IGLL5 and IL2RA in Validation Cohort
We utilized RCC data from the ICGC and GEO database to validate the prognostic value as well. In RCC patients from the ICGC database, a decrease in overall survival was significantly related to the high expression of IGLL5 and IL2RA (p < 0.05, Figures 5F,G). We then explored the prognostic value of IGLL5 and IL2RA using data from GSE40912 (n 32 samples) and GSE2748 (n 28 samples) datasets. The results indicated that high expression of IGLL5 was significantly associated with poor Frontiers in Molecular Biosciences | www.frontiersin.org November 2021 | Volume 8 | Article 760031 7 prognosis (p < 0.05, Figures 5H,I). Patients with high expression of IL2RA had poor survival, but the difference was not statistically significant (p 0.112, Figure 5J) Furthermore, the robustness of IGLL5 as well as IL2RA as biomarkers was verified using primary ccRCC and pRCC samples from Shandong provincial Hospital affiliated to Shandong First University. For ccRCC, IHC images indicated that normal renal tissue had weak staining for IGLL5 and IL2RA, but relatively strong staining patterns for IGLL5 and IL2RA were observed in the cytoplasm of tumor tissues (Figures 6A,B). Moreover, similar staining patterns for IL2RA were observed in pRCC tissues and adjacent normal tissues ( Figure 6C). However, IGLL5 staining was weakly positive by IHC, and no significant differences was detected between pRCC tissues and adjacent normal tissues ( Figure 6D). These unique IHC staining patterns further confirmed the above results and illustrate that the two genes could be used to predict clinical outcome and distinguish cancerous tissue from normal tissue.

Pharmacological Analysis
We asked whether selected regulators have ability to predict drug response and potential treatment options for ccRCC and pRCC. Through CellMiner database and Connectivity Map (Cmap), the top 16 relevant correlations between IGLL5 and IL2RA and drug activity were shown in Figure 7. The response to all anticancer  Frontiers in Molecular Biosciences | www.frontiersin.org November 2021 | Volume 8 | Article 760031 8 drugs was negatively correlated to the expression level of IL2RA, including panobinostat, pralatrexate, AT-13387, belinostat, and ponatinib. However, except for irofulven, other anticancer drugs were positively associated with the expression of IGLL5 (p < 0.05, Supplementary Table S1). Moreover, potential drugs for ccRCC and pRCC were predicted based on the top 15 prognosis-related hub genes ( Figures 3E,F). The prediction results revealed that 325 and 47 candidate drugs exhibited promising clinical applications for ccRCC and pRCC, respectively (Score ≤ −0.75, Supplementary  Table S2, S3). The molecular structure of top 10 potentially active drugs for the treatment of two types of RCC was shown in Supplementary Figure S6.

DISCUSSION
RCC accounts for more than 90% of kidney cancers, which consists of over 10 histological subtypes. Among these, ccRCC and pRCC are the main subtypes with high incidence and high mortality (Comprehensive molecular c, 2013; Chen et al., 2016;Linehan et al., 2016;Hsieh et al., 2017). Given limited options of effective therapy and modest effects of immunotherapy, there is an urgent need for novel biomarkers for future clinical application and evaluation (Barata and Rini, 2017).
Several studies have established that TME is a potential regulator of cancer progression and a source of therapeutic targets (Quail and Joyce, 2013;Vitale et al., 2019;Ho et al., 2020). In complex TME, immune cells and stromal cells play a significant role in cancer development. Stromal cells are genetically stable, and are considered as promising therapeutic targets for reducing tumor recurrence and drug resistance (Quail and Joyce, 2013). In TME, immune cells are divided into tumor antagonizing and promoting-immune cells. With the ability to regulate tumor progression in various stages (Lei et al., 2020). For example, inflammation could lead to CD8 + T cell dysfunction and consequently promote cancer progression (Wherry, 2011). Furthermore, in breast TME, NK cells are immature and dysfunctional NK in hepatocellular TME is related with cancer development (Krneta et al., 2016;Zhang et al., 2017). Garaud et al. found that infiltration of B cells in breast TME was correlated with a favorable prognosis (Garaud et al., 2019). However, several investigations have revealed that B cells may facilitate tumorigenesis by suppressing CD8 + T cells and producing cytokines, which contribute to angiogenesis (de Visser et al., 2005;Inoue et al., 2006;Tsou et al., 2016). Currently, the understanding of knowledge of TME of kidney cancers is restricted to specific tumor types and there is a lack of comprehensive analysis. Here, we used several analyses to explore the TME landscape of ccRCC and pRCC. Firstly, we calculated the immune scores as well as stromal scores and assessed the clinical relevance by integrating of clinical information, which revealed that the survival of ccRCC and pRCC patients was significantly correlated with high immune scores and/or stromal scores. Moreover, various clinical characteristics, such as tumor grade, tumor stage, T stage, M stage were associated with immune scores and stromal scores. Then, by investigating the associated biological functions and pathways via GO enrichment and KEGG pathway analyses, we uncovered diverse immunerelated functions and pathways, such as B cell mediated immunity, lymphocyte mediated immunity; cytokine-cytokine receptor interaction, chemokine signaling pathways, T cell receptor and B cell receptor signaling pathways. These results are consistent with previous studies and confirmed the essential role of immune cells in ccRCC and pRCC (Xu et al., 2019;Díaz-Montero et al., 2020;Lei et al., 2020;Wan et al., 2020;Zeng et al., 2020). Subsequently, according to STRING database, univariate Cox analysis, and Venn algorithm, we identified the two prognostic hub genes, IL2RA and IGLL5 in ccRCC and pRCC. To further exploit the potential role of these genes, we validated each individual gene in TCGA ccRCC and pRCC cohorts. The results revealed that IL2RA is a DEG and is correlated with prognosis of ccRCC and pRCC. However, IGLL5 was significantly expressed in ccRCC, but not in pRCC. In the ICGC, GEO40912, and GSE2748 validation cohort, we also observed a significantly differential survival trend. However, K-M analysis of OS indicated that high expression of IL2RA had a negative impact in GSE2748 dataset, but the difference was not significant. Moreover, the clinical application of the two genes was also validated using primary KIRC and KIRP samples from our hospital, and the results further confirmed the prognostic value of the two genes in clinical application. Previous research on IGLL5 and IL2RA focused predominantly on hematological malignancies. Fusion or mutation of IGLL5 was identified in Burkitt lymphoma, multiple lymphadenopathy, diffuse large B cell lymphoma, multiple myeloma, and chronic lymphocytic leukemia (Kasar et al., 2015;White et al., 2018;Panea et al., 2019;Li et al., 2020;Mosquera Orgueira et al., 2021). It has been reported that IGLL5 may act as a risk factor and was related with metastasis and poor prognosis of multiple myeloma (White et al., 2018;Barwick et al., 2019;D'Agostino et al., 2020). IGLL5 is also correlated with glioblastoma survival and breast cancer, which could serve as a biomarker for relapse-free survival in breast cancer with more than 85% accuracy (Ascierto et al., 2012;Qin et al., 2020). IL2RA also exhibits prognostic and diagnostic value in hematopoietic cancers and solid cancers. Elevated expression of IL2RA is associated with poor prognosis and may act as a biomarker in acute myeloid leukemia (Du et al., 2019;Nguyen et al., 2019). Up-regulated IL2RA was also found in tumor of lung, prostate, breast and melanoma, and thus, could serve as a marker of prognosis and immunoreaction (Kuhn and Dou, 2005).
Based on these results, we explored underlying mechanisms by GSEA analysis and also uncovered various immune-related pathways, further confirming the role of IGLL5 and IL2RA in immune regulation. Afterwards, CIBERSORT indicated that adaptive immune cells account for most of the immune cell infiltrations in ccRCC as well as pRCC. According to difference analysis and correlation analysis, we finally identified the immune cells commonly associated with IGLL5 and IL2RA in ccRCC and pRCC, such as CD8 + T cells, Tregs cells and T follicular helper cells, which was consistent with prior studies (Roth et al., 2018;Cañete et al., 2019;Ohue et al., 2019). Considering the unsatisfactory outcomes of therapy, we predicted small molecular compounds and individual response to drugs. Some commonly used drugs (such as oxaliplatin, ponatinib) have been shown to be influenced by IGLL5 and IL2RA, and numerous small molecular compounds have promising potential in future treatments. Currently, drugs targeting IL2RA consist of recombinant IL2, anti-IL2RA antibodies (such as basiliximab, daclizumab), anti-IL2RA radioimmunoconjuates (such as CHT-25, 90 Y-daclizumab), anti-IL2RA immunotoxins (such as LMB-2, RFT5-SMPT-dgA) and anti-IL2RA antibody drug conjugates (Flynn and Hartley, 2017). However, the efficacy and safety of some drugs is still unclear. Furthermore, IGLL5-targeted drugs have not been developed. Thus, further investigation of drugs targeting IGLL5 and IL2RA is of the essence.
In addition, there are other potential factors which also affect clinical outcomes in kidney cancers. Firstly, the spatial distribution of TME-regulated T cells may physically exclude these immune cells from the vicinity of cancer cells. Other essential considerations are that the stages of primary tumor and potentially different constitution of TME at metastatic sites (Fridman et al., 2012). To improve the outcome of ccRCC and pRCC, accurate and efficient biomarkers are indispensable. We investigated the TME and immune-related genes to identify associated immune cells and molecular mechanism as well as clinical value. Our findings suggested that IGLL5 and IL2RA could be used to estimate prognosis and may serve as potential therapeutic targets. The limited number of patients in GEO validation cohort may partially lead to bias. Therefore, in future investigations, expression data and clinical information from our medical center will be collected and further experiments using in vivo and in vitro models will be implemented to confirm the significant role of the two genes in the TME of kidney cancer.

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 authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Shandong Provincial Hospital Ethics Committee. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.