DOI: https://doi.org/10.21203/rs.3.rs-471905/v1
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.
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 [3–4]. 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 [8–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 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 [14–15], 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.
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.
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.
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.
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%).
| 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 | ||
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).
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.
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.
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.
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 [34–35]: (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.
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.
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.
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.