Prognostic value and immune infiltration of novel signatures in clear cell renal cell carcinoma microenvironment

Elevated immune and stromal score significantly correlated with advanced clinicopathological indicators and poor prognosis

Transcriptional expression profiles and phenotype data was download and integrated in 533 ccRCC patients from TCGA cohort. 64.7% (n=345) patients were male and 35.5% (n=188) were female. T1-T2 stage patients accounted for 64.1% (n=342) of the total number. N0 and N1 patients accounted for 45% (n=240) and 3% (n=16), respectively. M0 and M1 patients accounted for 79.2% (n=422) and 14.8% (n=79). In addition, after ESTIMATE algorithm was processed, stromal scores and immune scores were obtained, ranging from -2,716.84 to 4,773.7 and -1,158.94 to 3,076.7, respectively. Estimate score was significantly associated with higher ISUP grade and AJCC stage (Figure 1A, p<0.001, Figure 1B, p=0.0005). The highest Estimate score was found in the most progressive clinicopathological stage, G4 and stage IV. Immune score indicated significant prognostic implications, associated with elevated ISUP grade and AJCC stage (Supplementary Figure 1A–1B, p<0.0001). Stromal score significantly correlated with advanced ISUP grade (Supplementary Figure 1C, p=0.0463), while showed no association with AJCC stage (Supplementary Figure 1D, p=0.0674).

Figure 1. Association between immune/stromal/Estimate score and prognosis in TCGA after ESTIMATE algorithm processed. (A–B) Estimate score was significantly associated with higher ISUP grade and AJCC stage (p<0.001). The highest Estimate score was found in the most progressive clinicopathological stage, G4 and stage IV. (C) Survival curves indicated that elevated immune score significantly correlated with poor overall survival in 533 ccRCC patients (p=0.001, 1165 vs. 1217 days). (D) Increased stromal score significantly associated with shorter OS (p=0.002, 1117.5 vs. 1230 days). (E) Significant Estimate score also predict significant OS for ccRCC patients (p=0.003, 1172.5 vs. 1223.5 days).

To detect potential correlation between immune/stromal/Estimate score and survival benefits, we divided 533 patients into high and low score groups. Survival curves indicated that elevated immune score significantly correlated with poor overall survival (Figure 1C; p=0.001, 1165 vs. 1217 days). Increased stromal score significantly associated with shorter OS (Figure 1D; p=0.002, 1117.5 vs. 1230 days). Significant Estimate score also predict significant OS for ccRCC patients (Figure 1E; p=0.003, 1172.5 vs. 1223.5 days).

Differential expressed genes with immune and stromal score in ccRCC

To explore differential expressed genes (DEGs) profiles with immune and/or stromal scores, we performed transcriptional microarray analysis of 533 ccRCC cases from TCGA cohort. Based on immune score comparison, 162 genes were up-regulated and 747 genes down-regulated in the high score than the low score group after propensity analysis using limma package algorithm (Figure 2A). Similarly, for high stromal score compared with low score, 261 up-regulated genes and 1198 down-regulated genes were obtained (Figure 2B). A total of 77 DEGs were commonly upregulated in the high scores groups (Figure 2C), and 787 genes were synchronously downregulated using Venn algorithm (Figure 2D).

Figure 2. Differential expressed genes with immune and stromal score and related functional annotations in ccRCC. (A) Based on immune score comparison, 162 genes were up-regulated and 747 genes down-regulated in the high score than the low score group after propensity analysis using limma package algorithm. (B) Similarly, for high stromal score compared with low score, 261 up-regulated genes and 1198 down-regulated genes were obtained. (C–D) A total of 77 DEGs were commonly upregulated in the high scores groups, and 787 genes were synchronously downregulated using Venn algorithm. (E) functional enrichment analysis including GO: BP, GO: CC, GO: MF and KEGG pathways, was performed in 864 commonly DEGs. (F) Cluster analysis and heat map including 77 up-regulated DEGs suggested distinct mRNA expression profiles of DEGs in 533 ccRCC samples.

In addition, functional enrichment analysis including GO: BP, GO: CC, GO: MF and KEGG pathways, was performed in 864 commonly DEGs in Figure 2E. After –Log (FDR) sorting, we listed the top 10 function annotations of each part. As illustrated in Supplementary Figure 2, DEGs were mostly enriched in immune defense, plasma membrane, cytokine binding and cytokine-cytokine receptor interaction. Cluster analysis and heat map including 77 up-regulated DEGs suggested distinct mRNA expression profiles of DEGs in 533 ccRCC samples (Figure 2F).

Significant modular analysis based on PPI network

PPI network was constructed using a total of 77 commonly up-regulated DEGs in Figure 3A. MCODE, plug-in of Cytoscape, was used to detect most significant co-regulated modular. It indicated four closely related subgroups displayed in different color, and the most significant modular including AGPAT9, AQP7, HMGCS2, KLF15, MLXIPL, PPARGC1A, was marked in yellow. Survival curves of other nodes were illustrated in Supplementary Figure 3. It suggested that decreased SLC27A2, G6PC, MGAM, TRPM3, PKHD1, MYL3, MAPT, SLC22A6, TRHDE, TMEM174, SLC22A8, OGDHL, SCGN, SLC51B, SLC22A12, REN, PAH, GABRG1, SLC13A2, SST, KCNJ11 significantly correlated with poor OS, while elevated TUBB4A and RGS7 expression significantly predicted poor prognosis (p<0.05).

Figure 3. Significant modular analysis and function enrichment analysis based on PPI network. (A) PPI network was constructed using a total of 77 commonly up-regulated DEGs. MCODE, plug-in of Cytoscape, was used to detect most significant co-regulated modular. The most significant modular including AGPAT9, AQP7, HMGCS2, KLF15, MLXIPL and PPARGC1A, was marked in yellow. (B) functional annotations indicated that 77 DEGs were mostly involved in carbohydrate digestion and absorption, fatty acid transmembrane transport activity, PPAP signaling pathway, response to methionine, insulin resistance, water channel activity, enamel mineralization, negative regulation of mitochondrial fission, etc.

As shown in Figure 3B, functional annotations indicated that 77 DEGs were mostly involved in carbohydrate digestion and absorption, fatty acid transmembrane transport activity, PPAP signaling pathway, response to methionine, insulin resistance, water channel activity, enamel mineralization, negative regulation of mitochondrial fission, etc.

Survival analysis of significant DEGs in ccRCC from TCGA database

After integrating mRNA expression profile of six significant hub genes (AGPAT9, AQP7, HMGCS2, KLF15, MLXIPL, PPARGC1A) and clinical information, univariate regression analysis of overall survival were performed in TCGA cohort. As shown in Table 1, stromal score (ref. low), Immune score (ref. low), pT stage (ref. T1-T2), pN stage (ref. N0), pM stage (ref. M0), AJCC stage (ref. I-II), ISUP grade (ref. 1-2), AGPAT9, AQP7, HMGCS2, KLF15, MLXIPL and PPARGC1A expression (ref. low) were demonstrated as independent prognostic indicators for ccRCC patients (p<0.05). Multivariate Cox analysis showed that poor OS was significantly associated with pM stage (ref. M0; HR=2.807, p<0.001), ISUP grade (ref. 1-2; HR=1.765, p=0.029), MLXIPL expression (ref. low; HR=2.537, p=0.005) and PPARGC1A expression (ref. low; HR=0.468, p=0.009).

As demonstrated in Figure 4, among 6 significant hub genes, significantly decreased AGPAT9, AQP7, HMGCS2, KLF15, PPARGC1A mRNA expressions were found in ccRCC tissues compared with adjacent normal tissues, while MLXIPL mRNA expression was significantly elevated in tumor samples compared with normal samples. Kaplan-Meier method indicated that decreased AGPAT9, AQP7, HMGCS2, KLF15, PPARGC1A mRNA expression significantly correlated with poor OS (p<0.001), and elevated MLXIPL mRNA expression was significantly associated with shorter OS for ccRCC patients (p=0.012).

Figure 4. Survival analysis of significant DEGs in 533 ccRCC from TCGA database. Among 6 significant hub genes, significantly decreased AGPAT9, AQP7, HMGCS2, KLF15, PPARGC1A mRNA expressions were found in ccRCC tissues compared with adjacent normal tissues, while MLXIPL mRNA expression was significantly elevated in tumor samples compared with normal samples. Kaplan-Meier method indicated that decreased AGPAT9, AQP7, HMGCS2, KLF15, PPARGC1A mRNA expression significantly correlated with poor OS (p<0.001), and elevated MLXIPL mRNA expression was significantly associated with shorter OS for ccRCC patients (p=0.012).

Prognostic validation of MLXIPL and PPARGC1A in FUSCC cohort

To validate AQP9 mRNA expression profile in ccRCC tissues, we performed RT-qPCR using 380 paired tumor and normal samples with available clinical follow-up data from a real-world cohort. It revealed dramatically increased MLXIPL and decreased PPARGC1A mRNA expression in ccRCC samples than normal tissues (Figure 5A–5B). Survival curves suggested that patients with elevated MLXIPL and decreased PPARGC1A mRNA levels significantly correlated with poorer PFS and OS (p<0.001; Figure 5C–5F).

Figure 5. Prognostic validation of MLXIPL and PPARGC1A in FUSCC cohort. (A–B) To validate AQP9 mRNA expression profile in ccRCC tissues, we performed RT-qPCR using 380 paired tumor and normal samples with available clinical follow-up data from a real-world cohort. It revealed dramatically increased MLXIPL and decreased PPARGC1A mRNA expression in ccRCC samples than normal tissues. (C–F) Survival curves suggested that patients with elevated MLXIPL and decreased PPARGC1A mRNA levels significantly correlated with poorer PFS and OS (p<0.001).

Cox regression analysis and ROC curves

Multivariate Cox regression analysis of PFS and OS were performed in FUSCC cohort using mRNA expression profile of MLXIPL and PPARGC1A and clinicopathological information. As shown in Table 2, multivariate Cox analysis showed that poor PFS and OS were significantly associated with pT stage (ref. T1-T2), pN stage (ref. N0), pM stage (ref. M0), AJCC stage (ref. I-II), ISUP grade (ref. 1-2) and gene (MLXIPL or PPARGC1A) expression (ref. low) for ccRCC patients of FUSCC cohort (p<0.05).

After integrating all the significant clinicopathological parameters and gene expression profiles in the multivariate Cox regression models of FUSCC cohort, we generated the formula: Integrated score_(MLXIPL) = 2.105×pT stage (ref. T1-T2) + 1.922×pN stage (ref. N0) + 1.771×pM stage (ref. M0) + 2.413×AJCC stage (ref. I-II) + 1.934×ISUP grade (ref. 1-2) + 1.963×MLXIPL expression (ref. low) for PFS; Integrated score_(MLXIPL) = 1.840×pT stage (ref. T1-T2) + 1.832×pN stage (ref. N0) + 2.001×pM stage (ref. M0) + 3.434×AJCC stage (ref. I-II) + 1.764×ISUP grade (ref. 1-2) + 1.545×MLXIPL expression (ref. low) for OS. For PPARGC1A, we generated formula: Integrated score_(PPARGC1A) = 1.931×pT stage (ref. T1-T2) + 2.029×pN stage (ref. N0) + 1.641×pM stage (ref. M0) + 2.721×AJCC stage (ref. I-II) + 1.823×ISUP grade (ref. 1-2) + 0.524×PPARGC1A expression (ref. low) for PFS; Integrated score_(PPARGC1A) = 1.862×pT stage (ref. T1-T2) + 1.821×pN stage (ref. N0) + 1.912×pM stage (ref. M0) + 3.511×AJCC stage (ref. I-II) + 1.749×ISUP grade (ref. 1-2) + 0.665×PPARGC1A expression (ref. low) for OS. The AUC index of MLXIPL and PPARGC1A for the FUSCC-PFS and FUSCC-OS were 0.765, 0.768 and 0.778, 0.799, respectively (p<0.001; Figure 6A–6B). External validation was implemented in TCGA cohort. The AUC index of MLXIPL and PPARGC1A for the TCGA-PFS and TCGA-OS were 0.753, 0.750 and 0.748, 0.737, respectively (p<0.001; Figure 6C–6D). Survival curves suggested that integrated scores of MLXIPL and PPARGC1A expression significantly correlated prognosis in FUSCC cohort, and were validated significant in predicting prognosis in TCGA cohort (Supplementary Figure 4).

Figure 6. ROC curves were generated to validate the ability of the logistic model to predict prognosis. After integrating all the significant clinicopathological parameters and gene expression profiles in the multivariate Cox regression models of FUSCC cohort, we generated the formulas for MLXIPL and PPARGC1A to predict prognosis in FUSCC cohort, and validated prognostic ability in TCGA cohort.

Immune infiltration of MLXIPL and PPARGC1A

After identifying prognostic value of MLXIPL and PPARGC1A, we performed correlation analysis between MLXIPL and PPARGC1A and immune infiltration level for ccRCC in Figure 7. Scatter plots were generated with partial Spearman’s correlation and statistical significance. MLXIPL and PPARGC1A expression were significantly associated purity (correlation=0.207 and 0.287, respectively). In addition, elevated MLXIPL and PPARGC1A significantly correlated with B cell, CD8⁺ T cell, macrophage, neotrophil, and dendritic cell infiltration (p<0.05), prompting a general decline in immune infiltration level.

Figure 7. Immune infiltration of MLXIPL and PPARGC1A. After identifying prognostic value of MLXIPL and PPARGC1A, we performed correlation analysis between MLXIPL and PPARGC1A and immune infiltration level for ccRCC. Scatter plots were generated with partial Spearman's correlation and statistical significance. MLXIPL and PPARGC1A expression were significantly associated purity (correlation=0.207 and 0.287, respectively). In addition, elevated MLXIPL and PPARGC1A significantly correlated with B cell, CD8+ T cell, macrophage, neotrophil, and dendritic cell infiltration (p<0.05), prompting a general decline in immune infiltration level.

In Table 3, Spearman’s correlation and estimated statistical significance between MLXIPL, PPARGC1A expression and immune cell signature infiltration were displayed in detail. Correlation analysis between MLXIPL and PPARGC1A and immune cell infiltrations in ccRCC and normal samples were assessed in TCGA cohort in Supplementary Table 1. Partial correlation and correlation adjusted by tumor purity were also provided. Important signatures of a variety of immune cells include CD8⁺ T cell, T cell (general), B cell, Monocyte, tumor-associated macrophage (TAM), M1 Macrophage, M2, Macrophage, Neutrophils, Natural killer cell, Dendritic cell, Th1, Th2, Tfh, Th17, Treg, T cell exhaustion, were illustrated. (*p< 0.05; **p< 0.01; ***p< 0.001; ****p< 0.0001).

Ethics statement

All of the study designs and test procedures were performed in accordance with the Helsinki Declaration II. Study protocols were obtained by Fudan University Shanghai Cancer Center (FUSCC) (Shanghai, China) included in this research.

Raw biological microarray data

KIRC patients with available RNA-sequence data from TCGA database (https://tcga-data.nci.nih.gov/tcga/) were consecutively recruited in analyses [33]. The gene expression profile was measured experimentally using the Illumina HiSeq 2000 RNA Sequencing platform by the University of North Carolina TCGA genome characterization center. Level 3 data was downloaded from TCGA data coordination center, with available clinicopathological and survival data. ESTIMATE algorithm was used to calculate immune and stromal scores using "estimate" package (http://r-forge.r-project.org; repos=rforge, dependencies=TRUE) [9].

Patients and transcriptional expression profile

Clinicopathological parameters including ISUP grade and AJCC stage in 533 ccRCC patients from TCGA were analyzed and displayed according the immune, stromal and ESTIMATE score. One-way ANOVA test were utilized to measure statistically significance. X-tile software was utilized to take the cut-off value of immune score, stromal score and ESTIMATE score, in concordance of which overall participants were divided to two groups, respectively [34]. Survival comparison between distinct three scores identified as binary variables (high vs. low) was analyzed in 533 ccRCC patients. The primary end point for patients was progression-free survival (PFS), and overall survival (OS) was the secondary end point, which was evaluated from the date of first therapy to the date of death or last follow-up. The follow-up duration was estimated using the Kaplan-Meier method with 95% confidence intervals (95% CI) and log-rank test in distinct curves. All hypothetical tests were two-sided and P-values less than 0.05 were considered significant in all tests.

A total of 380 ccRCC patients from the Department of Urology of Fudan University Shanghai Cancer Center (FUSCC, Shanghai, China) from April 2009 to July 2018 were consecutively recruited in analyses. Tissue samples, including ccRCC and normal tissues, were collected during surgery and available from FUSCC tissue bank.

Identification, normalization and elucidation of DEGs

Preprocessing and normalization of raw biological data were the first step to process DNA microarray. This process removes bias from the microarray data to ensure its uniformity and integrity. Next, background correction, propensity analysis, normalization and visualization output of probe data were performed by robust multi-array average analysis algorithm in limma package [35]. Fold change > 1.5 and adj. p < 0.05 were set as the cut-offs to screen for differentially expressed genes (DEGs). Therefore, DEGs were identified based on |Log₂FC (fold change)|<0.5849 as statistically significance.

DAVID (http://david.ncifcrf.gov, Version 6.8) online database was performed to explore the role of development-related signaling pathways in ccRCC [36]. P-value<0.05 was considered statistically significant. Function annotations including biological processes (BP), molecular functions (MF), and cellular component (CC) were extracted from Gene Ontology (GO) enrichment analysis [37] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [38]. Hierarchical partitioning was performed using transcriptional expression profiles of selected positively-regulated DEGs in a heat map. Color gradients suggest high (red) or low (blue) expression level.

Protein-protein interaction (PPI) network and functional annotations

In this study, Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 10.0) online database was used to predict PPI network of significantly positive DEGs and analyze the degree of interactions between proteins [39]. Statistically significant edger was considered as interaction combined score>0.4 (medium confidence). Cytoscape (version 3.5), an open-access bioinformatics software platform providing the possibility of molecular maps, was utilized to visualize interactive network data [40]. Molecular Complex Detection (MCODE) (version1.4.2) is a plug-in for Cytoscape used for clustering a given network based on topology to find densely connected regions [41], which is able to identify the most significant module in the PPI networks with selection as follows: MCODE scores>5, degree cut-off=2, node score cut-off=0.2, Max depth=100 and k-score=2.

ClueGO is a Cytoscape plug-in that visualizes the non-redundant biological terms for large clusters of genes in a functionally grouped network [42]. GO: biological process and KEGG pathways analysis of selected hub genes were enrolled and visualized using ClueGO (version 2.5.3) and CluePedia (version 1.5.3), a functional extension of ClueGO, plug-in of Cytoscape [43].

Real-time quantitative PCR analysis

Total RNA sequence was extracted using TRIzol^® reagent (Invitrogen Life Technologies, USA) from 380 paired tumor and para-carcinoma normal samples according to manufacturer’s protocol. Then, total RNA from each sample was reverse transcribed to cDNA using the PrimeScriptTM RT reagent Kit (Takara Bio Inc., Japan). Primers were diluted in ddH₂O with SYBR Green PCR Master Mix (Applied Biosystems, Japan) according to the manufacturer’s instructions. The forward and reverse PCR primers for MLXIPL (ChREBP) are 5′- AAAACTGGGTTCTGGGTGTTC -3′ and 5′- AGGGAGTTCAGGACAGTTGG -3′, respectively. The forward and reverse primers for PPARGC1A are 5′- TGAACTGAGGGACAGTGATTTC -3′ and 5′- CCCAAGGGTAGCTCAGTTTATC -3′, respectively. Transcriptional expression was determined and normalized to β-Actin expression, and then analyzed by the -ΔΔCt method. Relative expression in ccRCC was represented using the ratio in Tumor/Normal tissues (T/N) or Normal/Tumor tissues (N/T).

Hub genes selection and statistics analysis in two cohorts

The most co-regulated hub genes penal were strived from MCODE. Clinical and pathological parameters and transcriptional expression profiles of hub genes in 533 ccRCC patients from TCGA cohort and 380 patients from FUSCC cohort were analyzed and displayed. Expression of hub genes was respectively identified as binary variables (high vs. low) referring to median expression taking the cut-off value of each hub genes. Partial likelihood test from Cox regression analysis was developed to address the influence of independent factors on PFS and OS. The follow-up duration was estimated using the Kaplan-Meier method with 95% confidence intervals (95%CI) and log-rank test in separate curves. Student’s t tests were utilized to compare differential hub genes expression between tumor and normal tissues (*p<0.05, **p<0.01, ***p<0.001, ****p<0.0001). Univariate and multivariate analysis were performed with Cox logistic regression models to find independent variables, including stromal score (ref. Low), immune score (ref. Low), pT stage (ref. T1-T2), pN stage (ref. N0), pM stage (ref. M0), AJCC stage (ref. I-II), ISUP grade (ref. 1-2), and each hub genes expression (ref. Low).

Integrated score was identified as sum of the weight of each significant hub gene and significant clinicopathological prognostic indicators. Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) and GEPIA (http://gepia.cancer-pku.cn/detail.php) was used to perform comprehensive correlation analysis between tumor-infiltrating immune cells signatures and selected hub genes. All hypothetical tests were two-sided and p-values less than 0.05 were considered significant in all tests. All of these statistical analyses were performed in R or corresponding R packages survival and survminer.

Ethics approval

The Ethics approval and consent to participate of the current study was approved and consented by the ethics committee of Fudan University Shanghai Cancer center.