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

DOI: https://doi.org/10.21203/rs.3.rs-471905/v1

Abstract

Background

Tumor microenvironment has been implicated in the development and progression of cancers. However, the prognostic significance of tumor microenvironment-related genes in kidney renal clear cell carcinoma (KIRC) remains unclear.

Methods

In this study, we obtained and analyzed gene expression profiles 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 identified a list of differentially expressed genes which was significantly associated with prognosis in KIRC patients. Protein-protein interaction networks and functional enrichment analysis were both performed, indicating that these identified genes were related to the immune response.

Conclusions

The tumor microenvironment-related genes could serve as the potential biomarkers for KIRC.

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 figure is increasing annually [34]. Of all renal cell carcinomas, approximately 85% are classified 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 identified.

Tumor microenvironment, as the tumor cellular milieu, is composed of myriad different immune, mesenchymal, and endothelial cells as well as a variety of inflammatory and extracellular matrix proteins [810]. The tumor microenvironment [1112] 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 influence 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 [1415], highlighting the utility of such algorithms driven by large datasets. However, at present there have been insufficient investigations regarding how the tumor microenvironment influences 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.

Materials And Methods

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 identified 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 significant.

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 identified DEGs, with p-value < 0.05 as the cut-off for significance.

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 significance based on a p < 0.001 significance threshold. SPSS 25.0 (SPSS Inc., Chicago, IL, USA) and R v3.5.1 were used for statistical analyses (http:// www.r-project.org/) [19]. A protein-protein interaction (PPI) network was produced for identified 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.

Results

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 profiles 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%).

Table 1

Baseline characteristics of the patients

Variables

All patients (n = 537)

No.

%

Gender

   

Female

191

35.6

Male

346

64.4

Age at diagnosis

   

> 60

271

50.5

≤ 60

266

49.5

Grade

   

G1 + G2

244

45.4

G3 + G4

285

53.1

NA

8

1.5

Stage

   

I + II

326

60.7

III + IV

208

38.7

NA

3

0.6

T stage

   

T1 + T2

344

64.1

T3 + T4

193

35.9

Metastasis

   

M0

426

79.3

M1

79

14.7

MX

30

5.6

NA

2

0.4

Lymph node status

   

N0

240

44.7

N1

17

3.2

NX

280

52.1

NA, no available

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 significantly 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 significant (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 profiles and stromal/immune scores. We observed distinct gene expression profiles 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.

Individual DEGs strongly related to OS in KIRC patients

We next sought to explore how KIRC patient OS was related with DEGs using the Kaplan-Meier curves and log-rank tests. Figure 3 showed that 23 DEGs were significantly related to OS in log-rank test (p < 0.001). Among them, ADGRV1, CWH43, FREM1, GPAT3, HMGCS2, HSD11B2, LDHD, OGDHL, PPARGC1A, SLC22A6, SLC22A8 and SLC22A12 were all positively associated with OS; while APCDD1L, AQP9, CPA4, FCRL5, FDCSP, GJB6, MIXL1, OBP2A, PAEP, SLN and TNFSF13B were negatively associated with OS.

Functional analysis and protein-protein interactions of identified DEGs

To elucidate the biological function of identified 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 significantly enriched in hematopoietic cell lineage, cytokine-cytokine receptor interaction, NF-kappa B signaling pathway and primary immunodeficiency (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 significant 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 identification of reliable diagnostic and prognostics KIRC biomarkers would allow for better treatment planning and a significant improvement in overall patient survival. As such, it is of great significance 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 influence the differentiation and activation of local mesenchymal cells through direct contact and via secreted factors, transforming them into tumor-associated 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 protein-protein interaction networks. We found that several DEGs were significantly associated with KIRC patient prognosis, with all of these genes being associated with immune and inflammatory 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. (2) The inhibitory effect of cytokines on apoptosis. Some cytokines, such as IL-10, IL-6 and GM-CSF, could inhibit apoptosis and further promote KIRC development. (3) Counteraction of KIRC tumor tissue. KIRC cells are able to induce monocytes to produce IL-10 and PGE2, while PGE2 can induce monocytes and macrophages to produce IL-10, which in turn causes a Th1-Th2 shift by regulating the balance of Th1/Th2 cell subsets, thus inducing the development of KIRC. Meanwhile, as a class of small molecular secreted proteins that can chemotactic the directional movement of immune cells, CCL19 also plays a key role in inhibiting tumor proliferation, migration, and invasion [31]. By studying a lung cancer model [32], Hillinger S, et al. reported that intratumoral injection of CCL19 significantly reduced the large size of tumors, and also found that the antitumor effect of CCL19 in the model was related to T cell-mediated immune responses. Lu J, et al [33] also found that CCL19 expression in colon cancer was significantly lower than that in adjacent normal tissues in clinical specimens, and patients with CCL19 positivity in tumors had a higher 5-year survival rate. CCL19 may play an anti-tumor effect through the following pathways [3435]: (1) chemotactic T cells are involved in killing tumor cells; (2) inducing dendritic cells to enter tumor tissues; (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 confirm these predictions.

Several limitations of the study also exist. First, we recapitulated our findings 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 findings. 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 findings: (1) we identified 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.

Abbreviations

RCC: Renal cell carcinoma; KIRC: Kidney renal clear cell carcinoma; ESTIMATE: Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data; TCGA: The Cancer Genome Atlas; DEGs: Differentially expressed genes; OS: Overall survival; IL-10: Interleukin 10; FCs: Fold changes; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: Protein-protein interaction; GO: Gene Ontology; MF: Molecular functions; BP: Biological process; DAVID: Database for Annotation, Visualization and Integrated Discovery.

Declarations

Ethics approval and consent to participate

This study was approved by Ethics Committee of the Institutional Review Board of Cixi People’s Hospital of Zhejiang Ningbo and Ningbo Urology and Nephrology Hospital, Ningbo Yinzhou No 2. Hospital, Ningbo, China.

Consent for publication

Not applicable.

Availability of data and material

All data are fully available without restriction.

Competing interests

The authors have declared that no competing interests exist.

Funding

None

Author contributions

RJZ, YSS and MMW conceived and designed the project, acquired the data, and analyzed the data. YLT and MLH interpreted the data. RJZ and YSS wrote the paper. All authors read and approved the final manuscript.

Acknowledgements

Not applicable.

References

  1. Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R.L., Torre, L.A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 68, 394-424.
  2. Chow, W.H., Dong, L.M., and Devesa, S.S. (2010). Epidemiology and risk factors for kidney cancer. Nat Rev Urol 7, 245-257.
  3. Lessi, F., Mazzanti, C.M., Tomei, S., Cristofano, C.D., Minervini, A., and Menicagli, M., et al. (2014). VHL and HIF-1α: gene variations and prognosis in early-stage clear cell renal cell carcinoma. Med Oncol 31,
  4. Zbar, B., Klausner, R., and Linehan, W.M. (2003). Studying cancer families to identify kidney cancer genes. Annu Rev Med 54, 217-233.
  5. Posadas, E.M., and Figlin, R.A. (2012). Systemic therapy in renal cell carcinoma: advancing paradigms. Oncology (Williston Park) 26, 290-301.
  6. Sejima, T., Iwamoto, H., Masago, T., Morizane, S., Hinata, N., and Yao, A., et al. (2013). Oncological and functional outcomes after radical nephrectomy for renal cell carcinoma: a comprehensive analysis of prognostic factors. Int J Urol 20, 382-389.
  7. Kuehn, H.S., Ouyang, W., Lo, B., Deenick, E.K., Niemela, J.E., and Avery, D.T., et al. (2014). Immune dysregulation in human subjects with heterozygous germline mutations in CTLA 4. Science 345, 1623-1627.
  8. Kobayashi, K., Omori, K., and Murata, T. (2018). Role of prostaglandins in tumor microenvironment. Cancer Metastasis Rev 37, 347-354.
  9. Park, S.A., and Surh, Y.J. (2017). Modulation of tumor microenvironment by chemopreventive natural products. Ann N Y Acad Sci 1401, 65-74.
  10. Abadjian, M.Z., Edwards, W.B., and Anderson, C.J. (2017). Imaging the tumor microenvironment. AdvExp Med Biol 1036, 229-257.
  11. Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., and Torres-Garcia, W., et al. (1999). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 4,
  12. Curry, J.M., Sprandio, J., Cognetti, D., Luginbuhl, A., Bar-ad, V., and Pribitkin, E., et al. (2014). Tumor microenvironment in head and neck squamous cell carcinoma. Semin Oncol 41, 217-234.
  13. Zhou, Z., and Lu, Z.R. (2017). Molecular imaging of the tumor microenvironment. Adv Drug Deliv Rev 113, 24-48.
  14. Priedigkeit, N., Watters, R.J., Lucas, P.C., Basudan, A., Bhargava, R., and Horne, W., et al. (2017). Exome-capture RNA sequencing of decade-old breast cancers and matched decalcified bone metastases. JCI Insight 2,
  15. Alonso, M.H., Aussó, S., Lopez-Doriga, A., Cordero, D., Guinó, E., and Solé, X., et al. (2017). Comprehensive analysis of copy number aberrations in microsatellite stable colon cancer in view of stromal component. Br J Cancer 117, 421-431.
  16. Publication Guidelines. http://cancergenome.nih.gov/publications/publicationguidelines.
  17. Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., and Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43,
  18. Jiao, X., Sherman, B.T., Huang, D.W., Stephens, R., Baseler, M.W., and Lane, H.C., et al. (2012). DAVID-WS: a stateful web service to facilitate gene/protein list analysis. Bioinformatics 28, 1805-1806.
  19. Team RC. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. 2014.
  20. Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., and Huerta‚ÄźCepas, J., et al. (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 43, 447-452.
  21. Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., and Ramage, D., et al. (2003). Cytoscape: a software environment for integratedmodels of biomolecular interaction networks. Genome Res 13, 2498-2504.
  22. Srigley, J.R., Delahunt, B., Eble, J.N., Egevad, L., Epstein, J.I., and Grignon, D., et al. (2013). ISUP renal tumor pane. the international society of urological pathology (ISUP) vancouver classification of renal neoplasia. Am J SurgPathol 37, 1469-1489.
  23. Barata, P.C., and Rini, B.I. (2017). Treatment of renal cell carcinoma: Current status and future directions. CA Cancer J Clin 67, 507-524.
  24. Cohen, H.T., and McGovern, F.J. (2005). Renal-cell carcinoma. N Engl J Med 353, 2477-2490.
  25. Turley, S.J., Cremasco, V., and Astarita, J.L. (2015). Immunological hallmarks of stromal cells in the tumour microenvironment.Nat Rev Immunol 15, 669-682.
  26. Du, G.H., Wang, J.K., Richards, J.R., and Wang, J.J. (2019). Genetic polymorphisms in tumor necrosis factor alpha and interleukin-10 are associated with an increased risk of cervical cancer. IntImmunopharmacol 66, 154-161.
  27. Hsu, T.I., Wang, Y.C., Hung, C.Y., Yu, C.H., Su, W.C., and Chang, W.C., et al. (2016). Positive feedback regulation between IL10 and EGFR promotes lung cancer formation. Oncotarget 7, 20840-54.
  28. Gulubova, M., Aleksandrova, E., and Vlaykova, T. (2018). Promoter polymorphisms in TGFB1 and IL10 genes influence tumor dendritic cells infiltration, development and prognosis of colorectal cancer. J Gene Med 20,
  29. Kim, J., Cho, Y.A., Choi, I.J., Lee, Y.S., Kim, S.Y., and Shin, A., et al. (2012). Effects of interleukin-10 polymorphisms, Helicobacter pylori infection, and smoking on the risk of noncardia gastric cancer. PLoS One 7,
  30. Wittke, F., Hoffmann, R., Buer, J., Dallmann, I., Oevermann, K., and Sel, S., et al. (1999). Interleukin 10 (IL-10): an immunosuppressive factor and independent predictor in patients with metastatic renal cell carcinoma. Br J Cancer 79, 1182-1184.
  31. Franciszkiewicz, K., Boissonnas, A., Boutet, M., Combadière, C., and Mami-Chouaib, F. (2012). Role of chemokines and chemokine receptors in shaping the effector phase of the antitumor immune response. Cancer Res 72, 6325-6332.
  32. Hillinger, S., Yang, S.C., Zhu, L., Huang, M., Duckett, R., and Atianzar, K., et al. (2003). EBV-induced molecule 1 ligand chemokine (ELC/CCL19) promotes IFN-gamma-dependent antitumor responses in a lung cancer model. J Immunol 171, 6457-6465.
  33. Lu, J., Zhao, J., Feng, H., Wang, P., Zhang, Z., and Zong, Y., et al. (2014). Antitumor efficacy of CC motif chemokine ligand 19 in colorectal cancer. Dig Dis Sci 59, 2153-2162.
  34. Okada, N., Gao, J.Q., Sasaki, A., Niwa, M., Okada, Y., and Nakayama, T., et al. (2004). Anti-tumor activity of chemokine is affected by both kinds of tumors and the activation state of the host's immune system: implications for chemokine-based cancer immunotherapy. Biochem Biophys Res Commun 317, 68-76.
  35. Braun, S.E., Chen, K., Foster, R.G., Kim, C.H., Hromas, R., and Kaplan MH., et al. (2000). The CC chemokine CK beta-11/MIP-3 beta/ELC/Exodus 3 mediates tumor rejection of murine breast cancer cells through NK cells. J Immunol 164, 4025-4031.