Identication of Tumor Microenvironment-Related Genes in Kidney Renal Clear Cell Carcinoma Patients via Bioinformatic Analysis

Background Tumor microenvironment has been implicated in the development and progression of cancers. However, the prognostic signicance of tumor microenvironment-related genes in kidney renal clear cell carcinoma (KIRC) remains unclear. Methods In this study, we obtained and analyzed gene expression proles from The Cancer Genome Atlas database. Stromal and immune scores were calculated based on the ESTIMATE algorithm. Results In the discovery series of 537 patients, we identied a list of differentially expressed genes which was signicantly associated with prognosis in KIRC patients. Protein-protein interaction networks and functional enrichment analysis were both performed, indicating that these identied genes were related to the immune response. Conclusions The could KIRC: cell carcinoma; ESTIMATE: Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data; The Cancer Genome Atlas; DEGs: expressed genes; 10;


Introduction
Cancer remains a primary cause of death globally [1], with the 8th most common form of cancer being renal cell carcinoma (RCC), making it a major public health threat [2]. Roughly 78,000 persons die of RCC worldwide each year, and this gure is increasing annually [3][4]. Of all renal cell carcinomas, approximately 85% are classi ed as kidney renal clear cell carcinoma (KIRC) [5]. Due to the slow development and protection of surrounding tissues, no early clinical symptoms of KIRC are revealed until the tumor volume is large enough or it enters advanced stages and distant metastasis [6]. While there have been several recent advances in surgical techniques and chemotherapeutic regimens available for treating KIRC, more advanced KIRC cases are not amenable to curative treatment efforts, with a 3-year overall survival rate of under 5% among KIRC patients [7]. Therefore, better understanding the local tumor microenvironment and the molecular mechanisms of KIRC tumors is essential, and the tumor biomarkers suitable for early detection should be identi ed.
Tumor microenvironment, as the tumor cellular milieu, is composed of myriad different immune, mesenchymal, and endothelial cells as well as a variety of in ammatory and extracellular matrix proteins [8][9][10]. The tumor microenvironment [11][12] would play a profound role in regulating gene expression of tumor tissues, and it is also related to tumor development, progression, chemoresistance, and metastasis [13]. Given its ability to in uence different outcomes in a context-dependent manner, the tumor microenvironment is likely to contain certain molecular biomarkers suitable for diagnostic or prognostic detection of particular cancers. Meanwhile, increasing reports have applied the ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) algorithm to breast and colon cancers [14][15], highlighting the utility of such algorithms driven by large datasets. However, at present there have been insu cient investigations regarding how the tumor microenvironment in uences KIRC outcomes. As such, we analyzed microenvironment-related gene expression as a means of assessing their prognostic value in patients with KIRC in this present study.

Data preparation
Gene expression data for KIRC patients was downloaded from The Cancer Genome Atlas (TCGA) data portal (http://cancergenome.nih.gov/), along with corresponding patient demographic details such as age, gender, outcomes, and survival. The ESTIMATE algorithm was then used in order to calculate both stromal and immune scores for the downloaded datasets [11]. The Ethics Committee of the Institutional Review Board of People's Hospital of Beilun district, Cixi People's Hospital of Zhejiang and Ningbo Yinzhou Second Hospital all approved this study. This study was conducted in a manner consistent with TCGA publication guidelines [16]. The exclusion criteria were as follows: (1) samples without differentially expressed genes (DEGs) and detailed clinical characteristics; (2) samples without prognostic data; (3) samples in which follow-up time were < 1 month. In total, 537 KIRC patients were enrolled in this study. DEGs were identi ed using an R package, with the limma package used for data analysis [17]. Fold change (FC) values for individual genes were determined, and DEGs with FC > 1.5 and p < 0.05 were deemed signi cant.
Heatmaps and enrichment analyses R packages were used for generation of heatmaps and for clustering analyses. The online Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics tool was used for gene ontology (GO) analyses, based on their molecular functions (MF) and biological processes (BP) (https://david.ncifcrf.gov/) [18]. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway was also performed for identi ed DEGs, with p-value < 0.05 as the cut-off for signi cance.

Overall survival curves and PPI network construction
The link between patient OS and the expression of particular genes was analyzed via a Kaplan-Meier analysis, with a log-rank test used to gauge signi cance based on a p < 0.001 signi cance threshold. SPSS 25.0 (SPSS Inc., Chicago, IL, USA) and R v3.5.1 were used for statistical analyses (http:// www.rproject.org/) [19]. A protein-protein interaction (PPI) network was produced for identi ed DEGs through the STRING database [20], with Cytoscape software used for visualization purposes [21]. Select protein pairs from this network with > 10 nodes were the outputs of this analysis.

Characteristics of KIRC patients
A total of 537 KIRC patients that had initially been diagnosed between 1989 and 2019 were enrolled in this study, with corresponding gene expression pro les and clinical data downloaded from TCGA. Table 1 list the detailed clinical characteristics, which includes age, gender, stage, disease grade, size of original (primary) tumor (T stage), neoplasm metastasis (M stage) and lymph node involvement (N stage). Of the enrolled patients, 64.4% were male, and 50.5% were order than 60 years. The most common tumor grades were G3 and G4 (53.1%). Similar, stage I and II were also listed as a major tumor stage (64.1%).

Relationship between ESTIMATE scores and overall survival in KIRC patients
The potential correlations of clinical information with stromal/immune scores were displayed in Fig. 1.
We found that patients > 60 years old had signi cantly lower stromal scores than younger patients; while males had higher immune scores than did females ( Fig. 1a and b, p < 0.05). With respect to tumor grade, G4 cases had the highest average immune scores, followed by G3, G2, and then G1 cases (Fig. 1c, p < 0.001). With respect to tumor staging, immune scores were highest for Stage IV diseases, followed by Stage II, Stage III, and then Stage I (Fig. 1d, p < 0.001). Correlations between patient OS and stromal/immune scores were also assessed via separating patients into either high or low score groups based on the median scores and comparing survival between these groups ( Fig. 1e and f). We found that patients with lower immune scores had a longer median OS relative to those with higher immune scores (p < 0.05 in log-rank test). Consistently, the median overall survival of cases with the low score group of immune scores is longer than the cases in the high score group, although it was not statistically signi cant (p = 0.258 in log-rank test).
Comparison of stromal/immune scores with gene expression in KIRC All 537 KIRC cases were obtained from TCGA database in order to show the association between global gene expression pro les and stromal/immune scores. We observed distinct gene expression pro les when comparing patients in the high and low stromal and immune score groups (Fig. 2), A total of 259 genes were upregulated and 152 genes were downregulated when comparing the high stromal score group to the low stromal score group (FC > 1.5, p < 0.05), whereas 512 genes were upregulated and 147 downregulated when comparing the high immune score group to the low immune score group (FC > 1.5, p < 0.05). In addition, Fig. 2c, d showed that 48 and 47 genes were commonly up-and down-regulated respectively, which would be used for the subsequent analysis.

Functional analysis and protein-protein interactions of identi ed DEGs
To elucidate the biological function of identi ed DEGs, GO enrichment analysis and KEGG pathway were performed in this study. A functional enrichment clustering analyses of these DEGs indicated that they were strongly associated with the immune response. The GO terms were mainly enriched in immune response, cytokine binding and chemotaxis (Fig. 4). In addition, the KEGG pathways were signi cantly enriched in hematopoietic cell lineage, cytokine-cytokine receptor interaction, NF-kappa B signaling pathway and primary immunode ciency ( Figure S1), which were also associated with immune response. Further, PPI networks were also obtained using STRING tool. As depicted in Figure S2, IL10, CCL19, CD79A, IGLL5, CD19, IGJ and TNFRSF17 were the remarkable nodes in this network, as they had the most connections with other members of the module.

Discussion
Recently, the incidence and mortality of KIRC have risen sharply in the world [22]. According to literature reports, there were about 63990 newly diagnosed renal cancer patients in 2017 in the United States [23], most of which were KIRC, placing a signi cant burden on affected patients, as well as their families and on the healthcare system as a whole. Although the improvement of basic research and clinical treatment in KIRC has improved the survival in KIRC patients in recent years, there are still 25%-30% of patients with tumor metastasis, and nearly one-third of patients have tumor recurrence [24]. Patients with recurrent or metastatic KIRC have an extremely poor prognosis. The identi cation of reliable diagnostic and prognostics KIRC biomarkers would allow for better treatment planning and a signi cant improvement in overall patient survival. As such, it is of great signi cance to further understand the molecular mechanisms of occurrence, development and metastasis of KIRC and to explore new biomarkers.
Numerous experiments [25] have shown that stromal cells in the tumor microenvironment interact with tumor cells. Tumor cells in turn can directly in uence the differentiation and activation of local mesenchymal cells through direct contact and via secreted factors, transforming them into tumorassociated mesenchymal cells which then secret proteins that can aid tumor progression. Given the importance of the tumor microenvironment to tumor progression, it represents a potentially ideal target worthy of study as a means of determining the prognosis of cancer patients based on tumor microenvironment-related gene expression. Therefore, in this study we assessed tumor microenvironment-related genes that were strongly correlated with KIRC patient OS, and conducted pathway and functional enrichment analyses of these genes, in addition to grouping them into proteinprotein interaction networks. We found that several DEGs were signi cantly associated with KIRC patient prognosis, with all of these genes being associated with immune and in ammatory responses. These results thus highlighted to potential value of tumor microenvironment-related gene expression as a means of predicting clinical outcomes in individuals suffering from KIRC.
From the protein-protein interaction network in this study, IL-10 and CCL19 are highly interconnected nodes. IL-10 has currently been reported to be closely associated with proliferation, immune evasion, and prognosis in a variety of tumors, such as cervical cancer [26], lung cancer [27], colorectal cancer [28], noncardia gastric cancer [29] and so on. Wittke F, et al [30] also reported that IL-10 was an immunosuppressive factor and independent predictor in patients with metastatic renal cell carcinoma.
Exactly how IL-10 functions in the development of renal clear cell tumors remains uncertain, with multiple mechanisms having been proposed to date. (1) Direct effects of cytokines. IL-10 is able to act directly on T cells, impairing Th1 cytokine production and enhancing Th2 cell differentiation, thus shifting the balance between Th1 and Th2 cells. (3) directly binding to CCR7 receptor and activating the receptor pathway to inhibit tumor proliferation, invasion and metastasis. In addition, our results further suggested that CD79A, IGLL5, CD19, IGJ and TNFRSF17, which have not previously been linked with KIRC prognosis, could serve as potential biomarkers for KIRC, and further molecular studies are needed to con rm these predictions.
Several limitations of the study also exist. First, we recapitulated our ndings in a single published dataset, which may impact the stability of the study and differ from the general population; a larger and multi-centric cohort study would be performed to replicate and validate our ndings. Second, the mechanisms behind the prognostic values of the tumor microenvironment-related genes were not investigated; experimental studies should be performed to provide information about their functions with regards to KIRC.

Conclusions
In conclusion, our study reported the following new ndings: (1) we identi ed some tumor microenvironment-related genes from functional enrichment analysis based on the immune/stromal scores; (2) These microenvironment-related genes could serve as the predictor of KIRC prognosis; (3) Some of the previously ignored genes have the potential to become additional biomarkers for KIRC. Finally, the study underscores the role of tumor microenvironment-related genes in KIRC patients and provides a new insight of the potential relationship between tumor microenvironment and KIRC prognosis. Availability of data and material All data are fully available without restriction.

Competing interests
The authors have declared that no competing interests exist. (e) KIRC cases were divided into two groups based on their immune scores: median survival of the low score group is longer than high score group, p = 0.044. (f) KIRC cases were divided into two groups based on their stromal scores: median survival of the low score group is longer than high score group; however, it is not statistically different, p = 0.258.