- Open Access
Identification of an at-risk subpopulation with high immune infiltration based on the peroxisome pathway and TIM3 in colorectal cancer
BMC Cancer volume 22, Article number: 44 (2022)
Peroxisomes are pivotal metabolic organelles that exist in almost all eukaryote cells. A reduction in numbers and enzymatic activities of peroxisomes was found in colon adenocarcinomas. However, the role of peroxisomes or the peroxisome pathway in colorectal cancer (CRC) is not defined.
In the current study, a peroxisome score was calculated to indicate the activity of the peroxisome pathway using gene set variant analysis based on transcriptomic datasets. CIBERSORTx was chosen to infer enriched immune cells for tumors among subgroups. The SubMap algorithm was applied to predict its sensitivity to immunotherapy.
The patients with a relatively low peroxisome score and high level of T-cell immunoglobulin and mucin domain 3 (TIM-3) presented the worse overall survival than others. Moreover, low peroxisome scores were associated with high infiltration of lymphocytes and poor prognosis in those CRC patients. Thus, a PERLowTIM3High CRC risk subpopulation was identified and characterized by high immune infiltration. The results also showed that CD8 T cells and macrophages highly infiltrated tumors of the PERLowTIM3High group, regardless of consortium molecular subtype and microsatellite instability status. This subgroup had the highest tumor mutational burden and overexpression of immune checkpoint genes. Further, the PERLowTIM3High group showed a higher probability of responding to programmed cell death protein-1-based immunotherapy. In addition, genes involved in peroxisomal metabolic processes in CRC were also investigated since peroxisome is a rather pleiotropic and highly metabolic organelle in cell. The results indicated that only those genes involved in fatty acid alpha oxidation could be used to stratify CRC patients as similar as peroxisome pathway genes.
We revealed the favorable prognostic value of the peroxisome pathway in CRC and provided a new CRC stratification based on peroxisomes and TIM3, which might be helpful for CRC diagnostics and personalized treatment.
Colorectal cancer (CRC), one of the most prevalent malignancies, is the second leading cause of cancer-related death worldwide . Although the use of chemotherapy and targeted therapy has dramatically prolonged survival time for patients with unresectable CRCs, the efficacy of these therapies can be limited by drug resistance and side effects . In the past few years, a better knowledge of the complex interactions between tumor cells and the immune system has led to the establishment of novel immunotherapies [3, 4]. However, only a small subset of metastatic CRCs with the microsatellite instability (MSI) phenotype benefit from immune checkpoint inhibitors (ICIs) that target co-inhibitory receptors [5, 6]. One major immunotherapeutic obstacle of CRC is the immunosuppressive tumor microenvironment (TME), a highly complex and heterogenous contexture [7, 8]. For CRC tumors with pre-existing strong immune infiltration, restoring the suppressive effects of the TME on immune cells is a potential strategy to improve the efficacy of immunotherapy . Thus, TME-based classifications are valuable for personalized immunotherapy in CRC.
Peroxisomes are membrane-bound organelles first discovered by Johannes Rhodin in 1954. Peroxisomes are involved in multiple metabolic processes, including the metabolism of fatty acid oxidation and D-amino acids, biosynthesis of ether lipids and bile acids, and turnover of reactive oxygen species (ROS) . Apart from being a highly metabolic organelle, the peroxisome has emerged as a pivotal regulator of inflammation and immune responses in these few years [11,12,13]. Furthermore, a recent study suggested that enhanced peroxisome signaling might be involved in reprograming the lipid metabolism of T cells in urological cancers . The attention on the role of peroxisomes or peroxisomal lipid metabolism in cancer, a disease characterized by metabolic reprogramming, is continuously increasing. From earlier publications, an overall reduction in peroxisomal proteins and enzymatic activities has been found in neoplastic tissue, including that of the colon [15,16,17]. Moreover, a more recent study suggested that NUDT7, a peroxisomal gene, is a potent tumor suppressor during KRASG12D-driven CRC development . Although these publications pointed to a potential loss of peroxisome function during CRC tumorigenesis, our current view on the role of peroxisomes or the peroxisome pathway in CRC remains fragmentary and limited.
In the current study, the prognostic role of peroxisome pathway was evaluated in CRC based on the genomic data of public datasets. In addition, an at-risk CRC subpopulation was identified and characterized by massive immune infiltration. Then, we explored the relationship between the newly characterized subpopulation and immunotherapy efficiency. Finally, we further investigated whether the genes specifically involved in peroxisomal metabolic processes can stratify CRC patients individually.
Data collection and processing
This study used two public patient datasets for discovery and validation. The CRC cohort from the Cancer Genome Atlas (TCGA) was used for discovery. We downloaded the RNA-Seq expression data (FPKM values) of 593 primary colon or rectum tumors with overall survival (OS) data using the GDCquery of the R package TCGAbiolinks . We then transformed the FPKM values into transcripts per kilobase million (TPM) values for subsequent analysis. We used a consensus measurement of tumor purity to remove samples with low tumor purity (< 60%) . We updated the clinical information from TCGA Pan-Cancer Clinical Data Resource . A total of 504 TCGA CRC patients enrolled in this study possessed gene somatic mutation data (MAF files). We analyzed and summarized the mutation data using maftools . A cohort from the Gene Expression Omnibus (GEO), GSE39582 (n = 556), was used for validation. We downloaded the normalized expression data and clinical information for the GSE39582 dataset from the GEO repository. We obtained the MSI status and consortium molecular subtype (CMS) classifications from the Colorectal Cancer Subtyping Consortium synapse data portal (https://www.synapse.org/#!Synapse:syn2623706/files/). We further combined MSI-low and microsatellite stable (MSS) cases and defined them as MSS and defined MSI-high cases as MSI for this analysis.
Gene set expression analyses
Gene set variant analysis (GSVA) is a computational approach that has been utilized to measure the overall activity of biological pathways [23,24,25]. With a pathway-based approach, the expression of pathway member genes can be summarized thus reducing complexity, and yielding more reliable results . We applied it to the transcriptomic data samples to estimate the HALLMARK_PEROXISOME score from the Molecular Signatures Database (MSigDB) Hallmark collection using the GSVA R package (version 3.10) . The median value of the GSVA HALLMARK_PEROXISOME score was used as the cut-off to categorize the CRC patients into either a peroxisome-high (Per-High) or peroxisome-low (Per-Low) group. We used the gene set enrichment analysis (GSEA) method to investigate the relationship between HALLMARK_PEROXISOME and other Hallmark gene sets using the GSEABase package (http://www.bioconductor.org/packages/release/bioc/html/GSEABase.html) . Gene set related to β-oxidation of very-long-chain fatty acids was obtained from Reactome Pathway database (https://reactome.org/) . Gene sets associated with fatty acid alpha oxidation and ROS were from Gene Ontology (http://geneontology.org/) [29, 30]. Gene sets related to biosynthesis of bile acids and ether lipid metabolism were from Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/kegg/kegg1.html) [31, 32].
Estimation of immune infiltration
The approach utilized for estimating immune infiltration was CIBERSORTx (https://cibersortx.stanford.edu/), which can predict the abundance and proportion of 22 types of tumor-infiltrating lymphocytes based on the expression of 547 genes for each tissue sample . We additionally applied the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumours using Expression data) algorithm to assess the enrichment of immune and stromal cells in the TME and to infer purity .
Prediction of the immunotherapy response for each subgroup
To predict the sensitivity of each subgroup to immunotherapy, we used Class mapping analysis (SubMap from Gene Pattern, https://www.genepattern.org/) to compare the similarity of the gene expression profiles of our subgroups to those of melanoma patients treated with checkpoint blockades against programmed cell death protein-1 (PD1) or cytotoxic T lymphocyte antigen-4 (CTLA4) .
We applied the optimal cut-off for T-cell immunoglobulin and mucin domain 3 (TIM3) expression with regard to the associated hazard of death events in a survfit model based on the cutp function of the survMisc R package (version 0.5.5; https://CRAN.R-project.org/package=survMisc). We chose the highest log-rank test score as the optimal cut-off for classifying patients into high- and low-expression groups with different risks. Kaplan–Meier analysis was used to estimate OS. We conducted multivariate Cox proportional hazards regression models to investigate the prognostic value of group III using the survival R package (version 2.42.3; https://CRAN.R-project.org/package=survival). To compare the contingency table variables, we used the chi-square test or Fisher’s exact test. We used the Mann–Whitney test for two-group comparison and the Kruskal–Wallis test for three-group comparisons. All statistical results with a p-value < 0.05 were considered significant.
A low peroxisome pathway score is associated with a worse clinical outcome and high immune cell infiltration in CRC patients
The peroxisome score was calculated by GSVA for CRC patients enrolled in the study and patients were categorized into either a Per-High or Per-Low group by the median of peroxisome scores. Intriguingly, the Per-Low group showed significantly worse OS than the Per-High group in TCGA CRC dataset (p = 0.024; Fig. 1a), findings that were validated in the GSE39582 dataset (p = 0.017; Fig. 1b). We next performed GSEA on 50 hallmark gene sets in TCGA cohort and found lipid-associated gene sets, CHOLESTEROL_HOMEOSTASIS and BILE_ACID_METABOLISM, which were significantly enriched in the Per-High group (Fig. 1c), which was validated in the GSE39582 cohort (Additional file 2: Fig. S1a). Furthermore, we found that immune-related gene sets were highly enriched in the Per-Low group, both in TCGA discovery dataset and the GEO-GSE39582 validation dataset (Fig. 1c; Additional file 2: Fig. S1a; Additional file 1: Table S1–2).
To further investigate the association between the peroxisome pathway and TME in CRC, we applied the ESTIMATE algorithm to TCGA cohort. We observed that the peroxisome pathway activity was negatively correlated with immune and stromal cell infiltration but was positively correlated with tumor purity (Fig. 1d). The Per-Low group had significantly higher immune and stromal scores but a lower level of tumor purity than the Per-High group (Fig. 1e). We also used a deconvolution method, CIBERSORTx, to assess the total immune infiltrates. The Per-Low group demonstrated higher immune cell infiltration than the Per-High group (Fig. 1f). A negative correlation between peroxisome scores and total immune infiltrates was observed (Fig. 1g). These results were validated with the GSE39582 dataset (Additional file 2: Fig. S1b-d). Taken together, our findings revealed that patients with low peroxisome scores had poor clinical outcomes despite having a high degree of immune infiltrates.
Identification and validation of an at-risk CRC subpopulation stratified by peroxisome score and TIM3 expression
To further investigate immunosuppressive mechanisms that could contribute to the poor prognosis observed in the Per-Low group, we analyzed the expression of immune checkpoints (ICs), as these genes are associated with the failure of T cells to efficiently eliminate tumor cells. The Per-Low group showed higher levels of all inhibitory checkpoint genes that we examined (Additional file 2: Fig. S2), including PDCD1 (PD1), PDL1, CTLA4, LAG3, TIGIT, TIM3, OX40, GITR, TNFRSF9 (4-1BB), ICOS, IDO1, CD40, CD70, and CD27. We observed that the peroxisome score was negatively correlated with the respective expression levels of these checkpoint genes in TCGA cohort (Fig. 2a). Among these ICs, TIM3 had the highest negative correlation with the peroxisome score. Additionally, we used an optimal cut-off for TIM3 expression to classify CRC patients into high and low TIM3 groups, followed by OS analysis. We found that for TCGA colorectal tumors, a high level of TIM3 expression was associated with a poor prognosis (Additional file 2: Fig. S3), which is consistent with the findings of previous studies [36, 37]. Intriguingly, Figs. 2b-e show that TIM3 prognostic behaviors were dependent on peroxisome scores. The median peroxisome score was used to classify patients into the Per-High and Per-Low groups. We applied an optimal cut-off based on the log-rank test for each independent group and subsequently defined four groups as follows: group I (PERHighTIM3High), group II (PERHighTIM3Low), group III (PERLowTIM3High), and group IV (PERLowTIM3Low) (Fig. 2b and c). In the Per-Low group, we found that TIM3 served as a prognostic marker of poor outcome; however, this was not seen in the Per-High group (Fig. 2d and e). Group III (PERLowTIM3High), with the highest proportion of MSI-status tumors, showed the most unfavorable outcomes compared with those in group IV and the Per-High group (group I + II), based on an analysis of OS (Fig. 3a and c). In addition, we adopted the same method to validate the existence of risk group III in the GSE39582 validation dataset. Group III consistently had the worst OS, as expected (Fig. 3b and d). Our results showed that group III patients had the worst OS; therefore, it would be of interest to explore differences in clinical and molecular characteristics across three subpopulations.
We gathered clinical and molecular data of 1149 CRC patients from TCGA and GSE39582. Characteristics of all patients are shown in Table 1. Advanced age and stages, a higher mutation frequency of BRAF, a CpG-island methylator phenotype (CIMP)-high status, and the MSI phenotype were associated with group III. Furthermore, we investigated the proportions of different CMS subtypes across subgroups and found that group III largely overlapped with CMS1 (MSI immune subtype; strong immune infiltration and activation) and CMS4 (mesenchymal; TGF-β activation; stromal invasion and angiogenesis).
To assess the independent poor prognostic value of group III, we performed univariate and multivariate cox regression analysis to investigate the association between traditional clinical and molecular characteristics and OS for the entire cohort. We found three independent indicators for poor OS, stage IV (Stage IV: HR = 6.07, 95% CI: 3.25–11.35, p < 0.001; Stage I as the reference), mutated KRAS (HR = 1.29, 95% CI: 1.00–1.66, p = 0.048) and group III (Group III: HR = 1.41, 95% CI: 1.05–1.90, p = 0.023; Group I + II as the reference). Whereas group III (PERLowTIM3High) had the highest proportions of CMS1 and CMS4 subtypes, of interest, our multivariable OS model showed that group III, but not CMS4, was an independent indicator for predicting poor OS (Table 2).
The PERLowTIM3High group has high levels of tumor infiltrated CD8+ T cells and macrophages
To explore the immune heterogeneity across subgroups, we used CIBERSORTx to estimate total immune infiltrates. Group III had the highest enrichment level of immune-infiltrated cells (Fig. 4a; Additional file 2: Fig. S4a). We then used the ESTIMATE-based approach to infer tumor purity; compared with group I + II and IV, group III had the lowest tumor purity (Fig. 4b; Additional file 2: Fig. S4b).
To further investigate the proportions of immune infiltrates in group III, we aggregated CIBERSORTx results of relative cell-types into six major cell types (B cells, CD4 T cells, CD8 T cells, dendritic cells, macrophages, and neutrophils). We observed that CD8 T cells and macrophages were highly enriched in group III both in TCGA and GSE39582 datasets (Fig. 4c; Additional file 2: Fig. S4c). CRC MSS tumors tend to have poor immunogenicity and lack responses to ICIs ; our results showed that group III MSS tumors were also highly infiltrated by CD8 T cells and macrophages (Fig. 4d and g; Additional file 2: Fig. S4d and S4g). As mentioned, group III had high proportions of CMS1 and CMS4, and these types were related to high immune infiltration. Next, we compared tumor infiltrating levels of CD8 T cells and macrophages between CMS1 + CMS4 patients among the three subgroups; CMS1 + CMS4 tumors in group III had the highest enrichment levels (Fig. 4e and h; Additional file 2: Fig. S4e and S4h). Moreover, other CMS subtypes (CMS2, CMS3, and indeterminate) tumors in group III had the highest infiltration of CD8 T cells and macrophages among the three subgroups (Fig. 4f and i; Additional file 2: Fig. S4f and S4i). Together, these results indicate that group III tumors were highly infiltrated by CD8 T cells and macrophages, which was not a simple reflection of high proportions of MSI tumors or CMS1/CMS4 phenotypes.
The PERLowTIM3High group has the highest tumor mutational burden (TMB) and is likely to respond to immunotherapy
To further explore whether these groups have different mutational profiles, we analyzed somatic mutation data from TCGA. Figure 5a shows the 20 most frequently mutated genes for each subgroup. We observed that the mutational load increased significantly in this group (Fig. 5b).
Given that group III was associated with the highest immune infiltration level and TMB, we next sought to assess the possibility of efficient immunotherapy for the three subgroups. We found that group III had the highest expression of ICs in both TCGA and GSE39582 cohorts (Fig. 5c-d), including PDL1, TIM3, LAG3, TIGIT, and VISTA. We then used subclass mapping to measure the similarity of expression profiles between the three CRC subgroups and a public cohort for melanoma patients treated with checkpoint blockades against CTLA4 or PD-1 . When comparing the expression profiles of the TCGA CRC subgroups and this external melanoma cohort, we observed that group III had a significant correlation with the PD-1 response group (Bonferroni corrected p = 0.012; Fig. 5e), a finding that was validated in the GSE39582 dataset (Bonferroni corrected p = 0.012; Fig. 5f). These results suggested that group III patients were likely to benefit from immunotherapy and specifically PD1-based strategies.
Genes involved in fatty acid α-oxidation stratified CRC patients similarly with the peroxisome pathway
Since peroxisomes are highly dynamic and pleiotropic metabolic organelles, we further investigated whether the genes specifically involved in peroxisomal metabolic processes can stratify CRC patients individually. Thus, we attempted to stratify TCGA CRC patients by using gene sets associated with specific peroxisome metabolic functions including REACTOME_BETA_OXIDATION_OF_VERY_LONG_CHAIN_FATTY_ACIDS (BOVLCFA), GOBP_FATTY_ACID_ALPHA_OXIDATION (FAAO), KEGG_PRIMARY_BILE_ACID_BIOSYNTHESIS (PBAB), KEGG_ETHER_LIPID_METABOLISM (ELM) and GOBP_CELLULAR_RESPONSE_TO_REACTIVE_OXYGEN_SPECIES (CRTROS) (Additional file 2: Fig. S5). For TCGA CRC patients with a relatively low BOVLCFA score or FAAO score, high TIM3 expression indicated worse OS (Additional file 2: Fig. S5a and S5b). However, our multivariable OS models showed that FAAOLowTIM3High, but not BOVLCFALowTIM3High, was an independent predictor for poor OS (Additional file 1: Table S3–4). Thus, among these five gene sets, FAAO was found as the only gene set that could stratify CRC patients similarly as HALLMARK_PEROXISOME (Additional file 2: Fig. S5b); TCGA CRC patients with a relatively low FAAO score and high TIM3 expression (FAAOLowTIM3High) had the worst OS, which was validated in GSE39582 (Fig. 6a and b). Similarly, FAAOLowTIM3High group tumors had the lowest tumor purity and were highly infiltrated by immune cells, especially CD8 T cells and macrophages (Fig. 6c-e). We also observed that in TCGA, FAAOLowTIM3High group tumors had significantly increased TMB and expression of ICs (Fig. 6f and g). Then, we applied the SubMap to compare the expression profiles of TCGA CRC subgroups stratified by FAAO and TIM3 with the external melanoma cohort mentioned previously herein. We found that the FAAOLowTIM3High group had a significant correlation with the PD-1 response group, which was validated in GSE39582 (Fig. 6h and i). Taken together, these results showed that genes involved in FAAO could be used to stratify CRC patients as similar as peroxisome pathway genes; moreover, tumors in the PERLowTIM3High group and the FAAOLowTIM3High groups displayed some similar characteristics.
Although there is a continuously increasing interest in peroxisomes and their functions in health and disease, the prognostic value of the peroxisome pathway has not been investigated. Our work first revealed a prognostic role of the peroxisome pathway in CRC, based on resources of genomic data from public clinical datasets. Then, we identified an at-risk CRC subpopulation based on the peroxisome pathway and TIM3 expression, and patients in this group might exhibit an improved response to immunotherapy. This new stratification was further refined by looking at the specific peroxisome functions, such as beta oxidation of very long-chain fatty acids, FAAO, biosynthesis of bile acids and ether lipids, and ROS turnover. Unexpectedly, we found that only genes involved in FAAO can stratify patients similarly with genes in the peroxisome pathway.
Peroxisomes are critical metabolic organelles associated with lipid metabolism and cellular ROS turnover . Cablé and colleagues found reduced numbers of peroxisomes in colon carcinoma by electron microscopy . However, the role of the peroxisome pathway in colorectal cancer has not been determined. As multiple peroxisomal biogenesis proteins and various oxidases are involved in peroxisome activity, we defined a peroxisome score to summarize the activity of the peroxisome pathway. We found that patients with a low peroxisome score had shorter OS than those with a high peroxisome score. Moreover, the score was associated with peroxisome-related pathways (e.g., the bile acid metabolism pathway and fatty acid metabolism pathway). Within the past few years, the idea that peroxisomes are vital organelles in regulating inflammation and the anti-microbial response has emerged [11,12,13, 39]. Vijayan and colleagues observed that deletion of the peroxisomal biogenesis factors PEX14 or MFP2 leads to a pronounced hyperexpression of COX2 and TNF-α proteins . Deletion of PEX13, another peroxisomal biogenesis factor, causes the activation of Smad-dependent TGF-β signaling and the release of inflammatory cytokines, such as TGF-β and IL-6 . Our results also showed that patients with a low peroxisome score had high immune infiltration and enriched immune-related pathways, despite the poor prognosis.
IC pathways are crucial in preventing autoimmunity and maintaining immune homeostasis; however, their upregulation in the TME of various malignancies can lead to T cell inactivation and immune evasion . By blocking the interactions between ICs and their ligands, ICIs can efficiently initiate anti-tumor immune responses and have shown promise in eliminating tumor cells in some cancers. TIM3 (encoded by Havcr2), as one IC, is an inhibitory receptor expressed by various immune cells and was initially identified as a cell surface marker specific to CD4+ T helper 1 and CD8+ T cells [42, 43]. In cancer, the upregulation of TIM3 expression marks the most terminally exhausted CD8+ T cell subsets [44, 45]. Previous studies have shown that TIM3 expression in the tissues of CRC patients is associated with tumor progression and poor clinical outcomes [36, 37]. Our results showed that for tumors with a low peroxisome score, high TIM3 expression indicated poor OS. Intriguingly, even if it was not statistically significant, high TIM3 expression seemed to be correlated with a better prognosis for tumors with a high peroxisome score, suggesting the possible role of the peroxisome pathway in regulating biological behaviors of TIM3. Further investigations are needed to validate this hypothesis. Our genomic analysis revealed that the prognostic performance of TIM3 was dependent on the peroxisome score, and in agreement with this, we identified an at-risk CRC subpopulation (PERLowTIM3High) with a negative prognostic value.
Beyond the negative prognostic value of group III (PERLowTIM3High), an activated immune phenotype characterized by high levels of pro-tumorous and anti-tumorous immune cells was observed for these tumors. Tumor-infiltrating CD8+ T cells are the most potent tumor-suppressing immune cells. CD8+ T cell abundance was an independent predictor of better disease-free and OS outcomes in early CRC patients, leading to the establishment of the Immunoscore, a reliable independent prognostic marker for CRC [46, 47]. Our data show that tumors in group III had the highest CD8+ T cell abundance. Given that high TIM3 expression is associated with dysfunctional CD8+ T cells, it is not surprising that group III had a poor prognosis in the presence of high CD8+ T cell infiltration. Compared with tumors in other subgroups, group III tumors had the highest level of macrophages, as one of the important types of tumor-supportive immune cells. Tumor-associated macrophages, mostly in the type M2 form, are important for the inhibition of anti-tumor immune responses mediated by T cells within the TME [48, 49]. The risk subgroup largely overlapped with CMS1 and CMS4 subtypes, which are associated with high immune infiltration, but the enrichment of CD8 T cells and macrophages was not a simple reflection of the high proportion of these two CMS subtypes. Moreover, our multivariate cox regression model determined the independent prognostic value of group III for OS. Hence, the peroxisome/TIM3 stratification is distinct from the CMS classification of CRC . We also found an association between group III tumors and MSI status. Our observation might promote refinement of the current clinical dogma that MSI-high tumors are associated with a good prognosis and encourage the further stratification of MSI tumors by peroxisome score and TIM3 expression [51, 52].
Consistent with the concept that MSI-high CRC is associated with high tumor mutations, group III tumors had the highest TMB. Moreover, high TMB was determined to be a predictive biomarker for immunotherapy in several tumor types [53, 54]. Immunotherapy is increasingly being recognized as a major treatment strategy for multiple cancers, including melanoma and a subset of CRC patients [55, 56]. PD1 blockade, as the only immunotherapeutic strategy approved by the FDA, has shown efficacy in metastatic CRC patients with MSI-high status [6, 57, 58]. Further, based on preclinical data, co-blockade of TIM3 and PD1 has shown efficacy for solid tumors, leading to the clinical investigation of anti-TIM3 combined with anti-PD1 in the treatment of various malignancies, including CRC [44, 59,60,61]. The sensitivity to immunotherapy is likely determined by many factors, such as the presence of CD8+ T cells, high TMB, and the expression of ICs [54, 62]. It is well established that metastatic CRC patients with an MSI-high status benefit from anti-PD1/PD-L1 therapy . Our results suggest that patients in group III (enriched with MSI tumors) had poor clinical outcomes but are likely to benefit from treatment targeting PD1. As such, the peroxisome pathway might be an overlooked factor involved in predicting the response to immunotherapy. In addition to anti-PD1 therapy, the co-blockade of PD1 and TIM3 could prove to be a more suitable treatment strategy for group III patients. Although immunotherapy approaches have been unsuccessful for MSS tumors, highly immune-infiltrated MSS tumors in group III could be more likely to respond to ICIs than MSS tumors in other subgroups. Further clinical trials investigating the co-blockade of PD1 and TIM3 in selected group III patients might be warranted.
Since peroxisomes are pleiotropic organelles involved in multiple metabolic processes, the meaning of high peroxisome scores seems obscure. Thus, we also attempted to refine the stratification by genes involved in peroxisomal metabolic processes. We found that CRC patients can be stratified by genes only involved in FAAO, which was similarly with genes in the peroxisome pathway. FAAO is a specific peroxisomal metabolic process in which some unusual fatty acids, such as 2-hydroxy fatty acids and 3-methyl-branched fatty acids (i.e., phytanic acid), are shortened by one carbon atom, which is known to occur in peroxisomes [64, 65]. Peroxisomal FAAO is deficient in patients with peroxisome biogenesis disorders, which is reflected by the accumulation of phytanic acid (PA) in the plasma . However, the role of FAAO or PA in CRC is not clear yet. Our results showed that only for CRC patients with a relatively low peroxisome score or FAAO score was high TIM3 expression indicative of poor OS, suggesting the possible role for normal peroxisomal FAAO in regulating biological functions of TIM3 in CRC, which is a new direction for our future work.
The retrospective nature of this analysis, conducted using public datasets, is the main limitation of our study. Furthermore, although we recapitulated the peroxisome pathway and identified the prognostic role of the peroxisome for CRC, assessing the correlation between the peroxisome score and the abundance of peroxisomes was outside the scope of our study. However, our findings suggested the value of the peroxisome pathway and the peroxisomal FAAO for CRC stratification, as well as the possible role of peroxisomes in TME and immunotherapy.
Collectively, we evaluated the prognostic role of the peroxisome pathway in CRC and identified an at-risk CRC subpopulation (PERLowTIM3High), exhibiting high immune infiltration. Patients with PERLowTIM3High tumors might more likely respond to anti-PD1 based immunotherapy. When looking at individual peroxisomal functions to refine the stratification, only genes involved in FAAO can stratify patients as similar as genes in the peroxisome pathway. Combined evaluation of the expression of TIM3 and genes involved in the peroxisome pathway or FAAO might be used for CRC diagnostics and could be helpful for personalized treatment.
Gene set variant analysis
T-cell immunoglobulin and mucin domain 3
Consortium molecular subtype
Tumor mutational burden
Programmed cell death protein-1
Fatty acid alpha oxidation
Immune checkpoint inhibitor
Reactive oxygen species
The Cancer Genome Atlas
Transcripts per kilobase million
Gene Expression Omnibus
The Molecular Signatures Database
Gene set enrichment analysis
Kyoto Encyclopedia of Genes and Genomes
Estimation of STromal and Immune cells in MAlignant Tumours using Expression data
programmed cell death protein-1
cytotoxic T lymphocyte antigen–4
CpG-island methylator phenotype
beta oxidation of very long chain fatty acids
Primary bile acid biosynthesis
Ether lipid metabolism
Cellular response to reactive oxygen species
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–49. https://doi.org/10.3322/caac.21660. Epub 2021 Feb 4.
Biller LH, Schrag D. Diagnosis and treatment of metastatic colorectal Cancer. Jama. 2021;325:669–85.
Khalil DN, Smith EL, Brentjens RJ, Wolchok JD. The future of cancer treatment: immunomodulation, CARs and combination immunotherapy. Nat Rev Clin Oncol. 2016;13:273–90.
Hodi FS, O'Day SJ, McDermott DF, Weber RW, Sosman JA, Haanen JB, et al. Improved survival with ipilimumab in patients with metastatic melanoma. N Engl J Med. 2010;363(8):711–23. https://doi.org/10.1056/NEJMoa1003466. Epub 2010 Jun 5. Erratum in: N Engl J Med. 2010;363(13):1290.
Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, et al. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science. 2017;357:409–13.
Le DT, Uram JN, Wang H, Bartlett B, Kemberling H, Eyring A, et al. Programmed death-1 blockade in mismatch repair deficient colorectal cancer. J Clin Oncol. 2016;34:15_suppl:103–103.
Kaymak I, Williams KS, Cantor JR, Jones RG. Immunometabolic interplay in the tumor microenvironment. Cancer Cell. 2020;39:28–37.
Peddareddigari VG, Wang D, DuBois RN. The tumor microenvironment in colorectal carcinogenesis. Cancer Microenviron. 2010;3:149–66.
Kubli SP, Berger T, Araujo DV, Siu LL, Mak TW. Beyond immune checkpoint blockade: emerging immunological strategies. Nat Rev Drug Discov. 2021;20(12):899–919. https://doi.org/10.1038/s41573-021-00155-y. Epub 2021 Mar 8.
Islinger M, Voelkl A, Fahimi HD, Schrader M. The peroxisome: an update on mysteries 2.0. Histochem Cell Biol. 2018;150:443–71.
Cara FD, Sheshachalam A, Braverman NE, Rachubinski RA, Simmonds AJ. Peroxisome-mediated metabolism is required for immune response to microbial infection. Immunity. 2017;47:93–106.e7.
Terlecky SR, Terlecky LJ, Giordano CR. Peroxisomes, oxidative stress, and inflammation. World J Biological Chem. 2012;3:93–7.
Vijayan V, Srinu T, Karnati S, Garikapati V, Linke M, Kamalyan L, et al. A new immunomodulatory role for peroxisomes in macrophages activated by the TLR4 ligand lipopolysaccharide. J Immunol. 2017;198:2414–25.
Mock A, Zschäbitz S, Kirsten R, Scheffler M, Wolf B, Herold-Mende C, et al. Serum very long-chain fatty acid-containing lipids predict response to immune checkpoint inhibitors in urological cancers. Cancer Immunol Immunother. 2019;68:2005–14.
Baur G, Wendel A. The activity of the peroxide-metabolizing system in human colon carcinoma. J Cancer Res Clin. 1980;97:267–73.
Cablé S, Keller JM, Colin S, Haffen K, Kédinger M, Parache RM, et al. Peroxisomes in human colon carcinomas. Virchows Archiv B. 1992;62:221.
Keller J, Cablé S, Bouhtoury F, Heusser S, Scotto C, Armbruster L, et al. Peroxisome through cell differentiation and neoplasia. Biol Cell. 1993;77:77–88.
Song J, Park S, Oh J, Kim D, Ryu JH, Park WC, et al. NUDT7 loss promotes KrasG12D CRC development. Cancers. 2020;12:576.
Mounir M, Lucchetta M, Silva TC, Olsen C, Bontempi G, Chen X, et al. New functionalities in the TCGAbiolinks package for the study and integration of cancer data from GDC and GTEx. PLoS Comput Biol. 2019;15(3):e1006701. https://doi.org/10.1371/journal.pcbi.1006701.
Aran D, Sirota M, Butte AJ. Systematic pan-cancer analysis of tumour purity. Nat Commun. 2015;6:8971.
Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An integrated TCGA Pan-Cancer clinical data resource to drive high-quality survival outcome analytics. Cell. 2018;173:400–416.e11.
Mayakonda A, Lin D-C, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747–56.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-Seq data. Bmc Bioinformatics. 2013;14:7.
Tokumaru Y, Oshi M, Katsuta E, Yan L, Satyananda V, Matsuhashi N, et al. KRAS signaling enriched triple negative breast cancer is associated with favorable tumor immune microenvironment and better survival. Am J Cancer Res. 2020;10:897–907.
Oshi M, Takahashi H, Tokumaru Y, Yan L, Rashid OM, Nagahashi M, et al. The E2F pathway score as a predictive biomarker of response to neoadjuvant therapy in ER+/HER2− breast Cancer. Cells. 2020;9:1643.
Khatri P, Sirota M, Butte AJ. Ten years of pathway analysis: current approaches and outstanding challenges. PLoS Comput Biol. 2012;8:e1002375.
Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS One. 2010;5:e13984.
Jassal B, Matthews L, Viteri G, Gong C, Lorente P, Fabregat A, et al. The reactome pathway knowledgebase. Nucleic Acids Res. 2019;48:D498–503.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.
Consortium TGO, Carbon S, Douglass E, Good BM, Unni DR, Harris NL, et al. The gene ontology resource: enriching a GOld mine. Nucleic Acids Res. 2020;49:D325–34.
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44:D457–62.
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.
Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37:773–82.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Roh W, Chen P-L, Reuben A, Spencer CN, Prieto PA, Miller JP, et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci Transl Med. 2017;9:eaah3560.
Yu M, Lu B, Liu Y, Me Y, Wang L, Zhang P. Tim-3 is upregulated in human colorectal carcinoma and associated with tumor progression. Mol Med Rep. 2017;15:689–95.
Kuai W, Xu X, Yan J, Zhao W, Li Y, Wang B, et al. Prognostic impact of PD-1 and Tim-3 expression in tumor tissue in stage I-III colorectal Cancer. Biomed Res Int. 2020;2020:1–13.
Ghiringhelli F, Fumet J-D. Is there a place for immunotherapy for metastatic microsatellite stable colorectal Cancer? Front Immunol. 2019;10:1816.
Facciotti F, Ramanjaneyulu GS, Lepore M, Sansano S, Cavallari M, Kistowska M, et al. Peroxisome-derived lipids are self antigens that stimulate invariant natural killer T cells in the thymus. Nat Immunol. 2012;13:474–80.
Oruqaj G, Karnati S, Vijayan V, Kotarkonda LK, Boateng E, Zhang W, et al. Compromised peroxisomes in idiopathic pulmonary fibrosis, a vicious cycle inducing a higher fibrotic response via TGF-β signaling. Proc National Acad Sci. 2015;112:E2048–57.
Toor SM, Nair VS, Decock J, Elkord E. Immune checkpoints in the tumor microenvironment. Semin Cancer Biol. 2020;65:1–12.
Monney L, Sabatos CA, Gaglia JL, Ryu A, Waldner H, Chernova T, et al. Th1-specific cell surface protein Tim-3 regulates macrophage activation and severity of an autoimmune disease. Nature. 2002;415:536–41.
Gayden T, Sepulveda FE, Khuong-Quang D-A, Pratt J, Valera ET, Garrigue A, et al. Germline HAVCR2 mutations altering TIM-3 characterize subcutaneous panniculitis-like T cell lymphomas with hemophagocytic lymphohistiocytic syndrome. Nat Genet. 2018;50:1650–7.
Sakuishi K, Apetoh L, Sullivan JM, Blazar BR, Kuchroo VK, Anderson AC. Targeting Tim-3 and PD-1 pathways to reverse T cell exhaustion and restore anti-tumor immunity. J Exp Medicine. 2010;207:2187–94.
Fourcade J, Sun Z, Benallaoua M, Guillaume P, Luescher IF, Sander C, et al. Upregulation of Tim-3 and PD-1 expression is associated with tumor antigen–specific CD8+ T cell dysfunction in melanoma patientsNY-ESO-1–specific T cells co-express Tim-3 and PD-1. J Exp Medicine. 2010;207:2175–86.
Galon J, Costes A, Sanchez-Cabo F, Kirilovsky A, Mlecnik B, Lagorce-Pagès C, et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science. 2006;313:1960–4.
Pagès F, Mlecnik B, Marliot F, Bindea G, Ou F-S, Bifulco C, et al. International validation of the consensus Immunoscore for the classification of colon cancer: a prognostic and accuracy study. Lancet. 2018;391:2128–39.
Liu Y, Cao X. The origin and function of tumor-associated macrophages. Cell Mol Immunol. 2015;12:1–4.
Qian B-Z, Pollard JW. Macrophage diversity enhances tumor progression and metastasis. Cell. 2010;141:39–51.
Guinney J, Dienstmann R, Wang X, de Reyniès A, Schlicker A, Soneson C, et al. The consensus molecular subtypes of colorectal cancer. Nat Med. 2015;21:1350–6.
Popat S, Hubner R, Houlston RS. Systematic review of microsatellite instability and colorectal Cancer prognosis. J Clin Oncol. 2004;23:609–18.
Ribic CM, Sargent DJ, Moore MJ, Thibodeau SN, French AJ, Goldberg RM, et al. Tumor microsatellite-instability status as a predictor of benefit from fluorouracil-based adjuvant chemotherapy for Colon Cancer. New Engl J Med. 2003;349:247–57.
Muzny DM, Bainbridge MN, Chang K, Dinh HH, Drummond JA, Fowler G, et al. Comprehensive molecular characterization of human colon and rectal cancer. Nature. 2012;487:330–7.
Samstein RM, Lee C-H, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019;51:202–6.
Marin-Acevedo JA, Dholaria B, Soyano AE, Knutson KL, Chumsri S, Lou Y. Next generation of immune checkpoint therapy in cancer: new developments and challenges. J Hematol Oncol. 2018;11:39.
Zappasodi R, Merghoub T, Wolchok JD. Emerging concepts for immune checkpoint blockade-based combination therapies. Cancer Cell. 2018;33:581–98.
Le DT, Uram JN, Wang H, Bartlett BR, Kemberling H, Eyring AD, et al. PD-1 blockade in tumors with mismatch-repair deficiency. New Engl J Med. 2015;372:2509–20.
Overman MJ, McDermott R, Leach JL, Lonardi S, Lenz H-J, Morse MA, et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (CheckMate 142): an open-label, multicentre, phase 2 study. Lancet Oncol. 2017;18:1182–91.
Zhou Q, Munger ME, Veenstra RG, Weigel BJ, Hirashima M, Munn DH, et al. Coexpression of Tim-3 and PD-1 identifies a CD8+ T-cell exhaustion phenotype in mice with disseminated acute myelogenous leukemia. Blood. 2011;117:4501–10.
Borate U, Esteve J, Porkka K, Knapper S, Vey N, Scholl S, et al. Phase Ib study of the anti-TIM-3 antibody MBG453 in combination with Decitabine in patients with high-risk myelodysplastic syndrome (MDS) and acute myeloid leukemia (AML). Blood. 2019;134(Supplement_1):570.
Ciardiello D, Vitiello PP, Cardone C, Martini G, Troiani T, Martinelli E, Ciardiello F. Immunotherapy of colorectal cancer: Challenges for therapeutic efficacy. Cancer Treat Rev. 2019;76:22–32. https://doi.org/10.1016/j.ctrv.2019.04.003. Epub 2019 May 4.
Gibney GT, Weiner LM, Atkins MB. Predictive biomarkers for checkpoint inhibitor-based immunotherapy. Lancet Oncol. 2016;17:e542–51.
Overman MJ, Lonardi S, Leone F, McDermott RS, Morse MA, Wong KYM, et al. Nivolumab in patients with DNA mismatch repair deficient/microsatellite instability high metastatic colorectal cancer: Update from CheckMate 142. J Clin Oncol. 2017;35(4_suppl):519.
Casteels M, Sniekers M, Fraccascia P, Mannaerts GP, Van Veldhoven PP. The role of 2-hydroxyacyl-CoA lyase, a thiamin pyrophosphate-dependent enzyme, in the peroxisomal metabolism of 3-methyl-branched fatty acids and 2-hydroxy straight-chain fatty acids. Biochem Soc T. 2007;35:876–80.
Jansen GA, Wanders RJA. Alpha-oxidation. Biochimica Et Biophysica Acta Bba - Mol Cell Res. 2006;1763:1403–12.
Wierzbicki AS. Peroxisomal disorders affecting phytanic acid α-oxidation: a review. Biochem Soc T. 2007;35:881–6.
We would like to thank the members of our laboratory for their suggestions and Yujian Yang for the critical review of the manuscript.
This work was supported by the National Natural Science Foundation of China (No: 81870390) and the Discipline and Platform Construction Project of Zhongnan hospital of Wuhan University (No: ZLYNXM202017).
Ethics approval and consent to participate
This research is based on the collection and analysis of publicly available data, and do not involve any human specimens and animal experiments.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
: Pathway enrichment analysis of peroxisome-high and peroxisome-low groups in TCGA dataset. Table S2: Pathway enrichment analysis of peroxisome-high and peroxisome-low groups in GEO-GSE39582 dataset. Table S3. Univariate and multivariate analyses of the relationship between subgroups stratified by FAAO/TIM3 and clinical characteristics. Table S4. Univariate and multivariate analyses of the relationship between subgroups stratified by BOVLCFA/TIM3 and clinical characteristics.
. a The volcano plot shows the enrichment analysis using Hallmark gene sets in the GSE39582 cohort. The red gene sets were enriched in the Per-High group, whereas the blue gene sets were enriched in the Per-Low group. b Correlation between the peroxisome score versus immune score, stromal score, and tumor purity in the GSE39582 cohort. c Boxplots of the immune score, stromal score, and tumor purity from ESTIMATE of Per-High and Per-Low groups in the GSE39582 cohort. d Boxplot of total immune infiltrates (sum of absolute scores across 22 immune cell types) and correlation between peroxisome score and total immune infiltrates of patients in TCGA cohort. For Boxplots, p values in group comparison with Mann-Whitney U-test are shown. For panels B, Pearson’s rho (r) and statistical difference (p) are indicated. **P < 0.01; ****P < 0.0001. Figure S2. Boxplots depict the expression of immune checkpoint genes in the TCGA colorectal dataset. Statistical P values between groups were determined by Mann-Whitney U-test. HM_PEROXISOME: Hallmark Peroxisome gene set. **P < 0.01, ***P < 0.001, ****P < 0.0001. Figure S3. a The log-rank test score at the candidate cut-off across the log-transformed TIM3 gene expression values is plotted. b Kaplan-Meier curves are plotted for the TIM3-High group and TIM3-Low group by the optimal cut-off shown in panel a. Figure S4. GSE39582 group III tumors were highly infiltrated with CD8 T cells and macrophages. a Violin plot showing the total immune infiltrates of CIBERSORTx for each subgroup in the GSE39582 CRC dataset. b Violin plot showing the ESTIMATE tumor purity for each subgroup in the GSE39582 CRC dataset. c Boxplots showing enrichment levels of CD8 T cells and macrophages for each subgroup in the GSE39582 dataset. d, g Boxplots of enrichment level of CD8 T cells (d) and macrophages (g) MSS tumors. e, h Boxplots of enrichment level of CD8 T cells (e) and macrophages (h) for GSE39582 CMS1 and CMS4 tumors. f, i Boxplots of enrichment level of CD8 T cells (f) and macrophages (i) for GSE39582 CMS2, CMS3, and indeterminate tumors. *P < 0.05; **P < 0.01; ***P < 0.001. Figure S5. a Scatter plot of BOVLCFA score and log2-transformed TIM3 gene expression values are shown for TCGA cohort. The Pearson’s rho (r) and statistical difference (p) are shown in the scatter plot. MSI (blue diamond) and MSS (black circle) status are labeled for CRC patients. The median value of TIM3 expression is indicated in a gray dashed line. The candidate cut-offs (black lines) for high BOVLCFA and low BOVLCFA groups are shown. Using the median value of BOVLCFA score and these two candidate cut-offs, we separated patients into four different groups labeled. Kaplan-Meier curves show the OS for the optimal cut-off of TIM3 expression in TCGA CRC subgroups with high BOVLCFA score and low BOVLCFA score, respectively. P-value was calculated by the log-rank test and shown for each plot. b The same as panel a, but for FAAO score. c The same as panel a, but for PBAB score. d The same as panel a, but for ELM score. e The same as panel a, but for CRTROS score. BOVLCFA: beta oxidation of very long chain fatty acids; PBAB: primary bile acid biosynthesis; ELM: ether lipid metabolism; CRTROS: cellular response to reactive oxygen species.
About this article
Cite this article
Yin, J., Wang, H., Hong, Y. et al. Identification of an at-risk subpopulation with high immune infiltration based on the peroxisome pathway and TIM3 in colorectal cancer. BMC Cancer 22, 44 (2022). https://doi.org/10.1186/s12885-021-09085-9
- Colorectal cancer (CRC)
- T-cell immunoglobulin and mucin domain 3 (TIM-3)
- Gene set variant analysis (GSVA)
- Fatty acid alpha oxidation (FAAO)