- Research
- Open access
- Published:
Novel hypoxia- and lactate metabolism-related molecular subtyping and prognostic signature for colorectal cancer
Journal of Translational Medicine volume 22, Article number: 587 (2024)
Abstract
Background
Colorectal cancer (CRC) is a serious global health burden because of its high morbidity and mortality rates. Hypoxia and massive lactate production are hallmarks of the CRC microenvironment. However, the effects of hypoxia and lactate metabolism on CRC have not been fully elucidated. This study aimed to develop a novel molecular subtyping based on hypoxia-related genes (HRGs) and lactate metabolism-related genes (LMRGs) and construct a signature to predict the prognosis of patients with CRC and treatment efficacy.
Methods
Bulk and single-cell RNA-sequencing and clinical data of CRC were downloaded from the TCGA and GEO databases. HRGs and LMRGs were obtained from the Molecular Signatures Database. The R software package DESeq2 was used to perform differential expression analysis. Molecular subtyping was performed using unsupervised clustering. A predictive signature was developed using univariate Cox regression, random forest model, LASSO, and multivariate Cox regression analyses. Finally, the sensitivity of tumor cells to chemotherapeutic agents before and after hypoxia was verified using in vitro experiments.
Results
We classified 575 patients with CRC into three molecular subtypes and were able to distinguish their prognoses clearly. The C1 subtype, which exhibits high levels of hypoxia, has a low proportion of CD8 + T cells and a high proportion of macrophages. The expression of immune checkpoint genes is generally elevated in C1 patients with severe immune dysfunction. Subsequently, we constructed a predictive model, the HLM score, which effectively predicts the prognosis of patients with CRC and the efficacy of immunotherapy. The HLM score was validated in GSE39582, GSE106584, GSE17536, and IMvigor210 datasets. Patients with high HLM scores exhibit high infiltration of CD8 + exhausted T cells (Tex), especially terminal Tex, and oxidative phosphorylation (OXPHOS)−Tex in the immune microenvironment. Finally, in vitro experiments confirmed that CRC cell lines were less sensitive to 5-fluorouracil, oxaliplatin, and irinotecan under hypoxic conditions.
Conclusion
We constructed novel hypoxia- and lactate metabolism-related molecular subtypes and revealed their immunological and genetic characteristics. We also developed an HLM scoring system that could be used to predict the prognosis and efficacy of immunotherapy in patients with CRC.
Background
Colorectal cancer (CRC) is currently the third most common cancer and the second leading cause of cancer-related deaths worldwide [1]. Growing evidence suggests that the tumor microenvironment (TME) plays an important role in the development, progression, and drug resistance of CRC.
Hypoxia is a prominent feature of the microenvironments of CRC and other solid tumors. The generation of hypoxic conditions may be related to the massive consumption of oxygen due to uncontrolled tumor proliferation and impaired oxygen supply due to irregular and disorganized tumor-associated neovascularization. Hypoxic conditions increase the levels of hypoxia-inducible factors (HIFs) in cells. As transcription factors, HIFs further increase the transcription of downstream target genes, thus playing a regulatory role in various biological processes, such as cell metabolism, proliferation, metastasis, epithelial-mesenchymal transition (EMT), and angiogenesis [2,3,4,5,6,7,8,9,10]. Additionally, hypoxia and HIF accumulation have been found to be associated with resistance to chemotherapy as well as worse prognosis in patients with CRC [11,12,13,14]. Recent studies have also found that hypoxia can be considered a biomarker for predicting the outcome of immunotherapy and that the efficacy of immunotherapy can be improved by ameliorating hypoxia and targeting HIF-1 [15,16,17,18,19]. Several studies have demonstrated that hypoxia regulates the function and differentiation of immunosuppressive cells, such as myeloid-derived suppressor cells, tumor-associated macrophages (TAMs), and regulatory T cells (Tregs), thereby promoting immunosuppression and tumor immune escape [20,21,22,23,24,25].
The hypoxia-mediated shift of tumor metabolism towards glycolysis, as well as the native characteristics of aerobic glycolysis (the so-called Warburg effect), leads to massive glucose consumption by tumor cells, consequently increasing lactate production and decreasing the pH of the TME [26]. Lactate can directly modulate endothelial cell phenotype and drive tumor angiogenesis through multiple pathways [27,28,29]. Furthermore, considerable lactate accumulation stimulates macrophage polarization to the M2 phenotype in tumors [30]. Lactate has been found to reduce interferon-γ (IFN-γ) production by CD8 + T cells and NK cells, inhibiting their cytotoxic function [31]. Lactate also decreases the motility of CD4 + and CD8 + T cells, which might reduce their infiltration and movement into the TME [32]. Additionally, lactate in the highly glycolytic TME may increase the expression of programmed cell death-1 (PD-1) in Tregs; PD-1 blockade therapy could activate PD-1-expressing Tregs, leading to immunotherapy failure [33]. Lactate has also recently been found to provide metabolic support to tumor-infiltrating Tregs [34]. Moreover, patients with metastatic CRC have been reported to have higher serum lactate concentrations than those with non-metastatic CRC [35].
However, considering the heterogeneity of CRC and the complex interaction between hypoxia and lactate metabolism, the effects of hypoxia and lactate metabolism on CRC have not been fully elucidated. Therefore, it is necessary to conduct a landscape assessment of the fundamental combination of hypoxia and lactate metabolism on CRC prognosis, TME, and immunotherapy. In this study, we performed a new subtype classification of CRC by combining hypoxia-related genes (HRGs) and lactate metabolism-related genes (LMRGs). We then identified the characteristics of the corresponding subtypes from multiple perspectives. Overall, this study aimed to determine a precise treatment strategy for CRC using this new subtype classification method, thereby improving the treatment efficacy and survival of patients with CRC.
Methods
Data acquisition
We obtained RNA sequencing data of CRC, as well as clinical information from the TCGA database (https://portal.gdc.cancer.gov/). After removing patients with less than 1Â month of follow-up, we included count and transcripts per million (TPM) data from 575 patients with CRC and 51 normal tissues. The GSE39582, GSE106584, and GSE17536 expression matrices and clinical information were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/). The IMvigor210 data were obtained from IMvigor210CoreBiologies (http://research-pub.gene.com/IMvigor210CoreBiologies). HRGs and LMRGs were obtained from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/downloads.jsp). Specifically, 326 HRGs were identified from Harris hypoxia, hallmark hypoxia, GOBP regulation of cellular responses to hypoxia, and reactome cellular responses to hypoxia. Overall, 288 LMRGs were obtained from the GOBP glucose catabolic process to lactate via pyruvate, GOBP lactate metabolic process, GOBP lactate transmembrane transport, GOMF L-lactate dehydrogenase activity, GOMF lactate dehydrogenase activity, GOMF lactate transmembrane transporter activity, HP abnormal brain lactate level by mrs, HP abnormal lactate dehydrogenase level, HP-elevated lactate pyruvate ratio, HP-increased circulating lactate dehydrogenase concentration, HP-increased CSF lactate level, and HP-increased serum lactate level.
Differential expression analysis
The R software package DESeq2 (version 1.32.0) was used for differential expression analysis. Genes with a false discovery rate (FDR) < 0.05 and |fold change|≥ 2 were screened as differentially expressed genes (DEGs). We analyzed the correlations among the 35 DEGs using Spearman correlation analysis. GeneMANIA was used to assess gene interactions and predict gene functions [36]. Additionally, we evaluated the prognostic significance of each gene using the Cox method with the survival package (version 3.5-5), which integrates survival time, survival status, and gene expression data.
Molecular subtyping
Cluster analysis was performed using ConsensusClusterPlus [37] with agglomerative PAM clustering with Euclidean distance; 80% of the samples were resampled for 10 repetitions. The optimal number of clusters was determined using an empirical cumulative distribution function plot and the average consistency within the group.
The prognostic differences between the groups were analyzed using the Survfit function of the R package survival (version 3.5-5). The log-rank test was used to assess the statistical significance, and Kaplan–Meier curves were plotted.
We classified the CRC samples in the TCGA database into four consensus molecular subtypes (CMS), following the report of Guinney et al. [38]. Gene Set Enrichment Analysis (GSEA) software (version 3.0) was utilized to perform GSEA. The h.all.v2023.1.Hs.symbols.gmt subset of the MSigDB was downloaded to assess the related pathways and molecular mechanisms.
Immune landscape analysis
We calculated 22 immune cell infiltration profiles per sample based on the CIBERSORT method in the R package IOBR [39]. Next, we analyzed and visualized the anticancer immune status of each sample and the proportion of tumor-infiltrating immune cells in the seven-step cancer immune cycle based on RNA-seq data in the Tracking Tumor Immunophenotype (TIP, http://biocc.hrbmu.edu.cn/TIP/) [40]. The immunophenoscore (IPS) and Tumor Immune Dysfunction and Exclusion (TIDE) score were used to assess the potential clinical efficacy of immunotherapies in the different subtypes and indicate the potential for tumor immune evasion. The IPS and TIDE scores for the samples were derived from The Cancer Immunome Atlas (TCIA, https://tcia.at/home/) and an online tool called TIDE (http://tide.dfci.harvard.edu/), respectively [41]. Finally, we compared the expression levels of well-known immune checkpoint genes across the three subtypes.
Prediction of patients’ clinical drug response
The R package oncoPredict (version 0.2) [42] was used to predict each patient’s response to multiple clinical medications based on cell line drug response and gene expression data from the Broad Institute’s Cancer Therapeutics Response Portal (CTRP) and Sanger Genomics of Drug Sensitivity in Cancer (GDSC).
Characterization of somatic mutations in patients
To characterize the somatic mutations in patients with different subtypes, mutation annotation format (MAF) files of patients were downloaded from the TCGA database. These files were analyzed and visualized using the R package maftools (version 2.10.05) [43].
Construction and validation of the prognostic model
Differential expression analysis was performed as previously described. Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses were performed using the R package clusterProfiler (version 3.14.3) to obtain the results of the gene set enrichment. The KEGG Rest API (https://www.kegg.jp/kegg/rest/keggapi.html) was also used. Genes in the R package org.Hs.eg.db (version 3.1.0) were used for GO annotation.
The R package randomForest (version 4.7-1.1) was used to develop optimal models. The importance of each explanatory variable was assessed by evaluating the mean decrease in accuracy and the Gini coefficient. The top 20 genes in terms of importance were selected for the least absolute shrinkage and selection operator (LASSO) regressions. The R package glmnet (version 4.1-8) was used to perform the LASSO analysis. A stepwise multivariate Cox proportional hazards model was used to optimize the model. Finally, we constructed an HLM scoring system based on six genes for prognostic prediction. HLM score = 0.43*(AGXT expression) + 0.03*(TRIB2 expression) + 0.17*(ELFN2 expression) + 0.16*(PCDHB10 expression) + 0.02*(CALCA expression) + 0.01*(CD79A expression). To validate the efficacy of the model, receiver operating characteristic (ROC) analysis was conducted at one, three, and five years using the R package Proc (version 1.17.0.1). The area under the curve (AUC) and confidence intervals were evaluated using the ci function of Proc.
To validate the stability of our model, we acquired GSE39582, GSE106584, GSE17536, and IMvigor210 expression matrices and clinical information. These data were divided into the high- and low-risk groups based on the HLM score using optimal truncation, and survival analysis was performed.
Single-sample gene set enrichment analysis (ssGSEA)
We obtained and defined the CD8 + exhausted T cell (Tex), GZMK + Tex, terminal Tex, oxidative phosphorylation (OXPHOS)−Tex, and TCF7 + Tex gene sets from previous studies (Supplementary Table 1) [44, 45], set the minimum gene set to 5 and the maximum gene set to 5000, and calculated the enrichment scores for each sample in each gene set using the R package GSVA (version 1.40.1).
Single-cell RNA sequencing (scRNA-seq) data processing and analysis
The scRNA-seq data used in this study were downloaded from the GEO database under the accession code GSE132465. This dataset contains single-cell 3'-RNA sequencing data from 23 patients with primary CRC [46]. The Seurat R package (version 5.0) was used to analyze scRNA-seq data according to standard analysis procedures [47]. Cells were selected for analysis based on the following criteria: cells with unique molecular identifier (UMI) counts greater than 1000, cells with more than 200 but less than 6000 unique genes, and cells with less than 20% mitochondrial gene expression in their UMI counts. The ElbowPlot function was used to determine the dimensionality of each dataset. The t-distributed stochastic neighbor embedding (t-SNE) projection was used to visualize the cell clusters. Major cell types and subtypes were annotated by comparing the typical marker genes and differentially expressed genes in each cluster.
Pseudo-bulk analysis of scRNA-seq data
Pseudo-bulk analysis of scRNA-seq data was performed to characterize the gene expression profile of each sample. Specifically, scRNA-seq data were converted to bulk-like data by aggregating the gene counts of individual cells belonging to each sample.
Construction of a predictive nomogram
To visualize the prognosis prediction in patients with CRC, we constructed a nomogram based on survival time, survival status, age, sex, T stage, N stage, M stage, pathological type, and HLM score using the R package rms (version 6.7-1). Harrell’s concordance index (C-index), calibration curves, and decision curve analysis (DCA) were used to assess the nomogram performance.
Protein interactions
Protein interactions mediate a range of physiological functions and pathological developments in organisms. We obtained the interaction networks of proteins that interacted with the characterized genes from the BioPlex Interactome database (https://bioplex.hms.harvard.edu/) [48]. The purple circles represent the quered protein, the green circles represent the bait protein, the gray diamonds represent the prey protein, and the arrows represent the directed edge (bait-to-prey).
Cell viability assay
The HCT116 and LS174T cell lines were purchased from the American Type Culture Collection (ATCC) website. The cells were cultured in 1640 medium (Invitrogen, C11875500BT) containing 10% fetal bovine serum (Vivacell, C04001-500) and 1% penicillin and streptomycin (Gibco, 15140122) at 37 °C in a 5% CO2 incubator. A hypoxic environment was created using a CO2 three-gas incubator (Thermo Fisher Scientific, 51901137); the O2 concentration was adjusted to 1%.
We seeded HCT116 and LS174T cells in 96-well plates (Costar 3599) at densities of 8,000 and 15,000 cells/well, respectively, and placed them in normoxic or hypoxic environments. The following day, different concentrations of 1640 complete medium, 5-fluorouracil (5-FU, Solarbio, F8300), oxaliplatin (Solarbio, O8390), or irinotecan (Solarbio, II0140), in combination with or without BAY87-2243 (a potent and selective HIF-1 inhibitor; Beyotime, SC1193), were added to the 96-well plates. The cells were then incubated under normoxic or hypoxic environments for another 48Â h.
For Trypan blue staining, the spent medium in the wells was discarded, 100 µl of 0.4% Trypan blue (LABLEAD, T6146) solution was added to each well, and then the Trypan blue stain was washed away with PBS (Invitrogen, C20012500BT) after 5 min. Finally, the cells were observed and photographed under the microscope.
For Cell Counting Kit-8 (CCK8) experiments, the spent medium in each well of the 96-well plate was discarded, 100 µl of 1640 medium containing 10% CCK8 (Beyotime, C0038) was added to each well, and the plates were incubated at 37 °C for 2 h. Finally, the absorbance at 450 nm was measured using an enzyme labeling instrument (Thermo Fisher Scientific, VLBL00D0).
Statistical analysis
SPSS 26.0 (SPSS, Inc., Chicago, IL, USA) and GraphPad Prism 8 (GraphPad, Inc., CA, USA) software were used for the analyses. Differences between two groups were calculated using the paired two-tailed Student’s t-test or the Mann–Whitney–Wilcoxon test. Comparisons among three groups were performed using ANOVA or the Kruskal–Wallis rank-sum test. The chi-square test was used to compare the clinical characteristics. Statistical significance was set at p < 0.05.
Results
Identification of hypoxia-associated DEGs and lactate metabolism-associated DEGs in CRC
Figure 1 illustrates the analytical framework used in this study. RNA-seq data were obtained from the TCGA database for 575 CRC samples with more than 1 month of follow-up, and 51 normal colorectal tissues were used for comparison. Differential expression analysis was performed, and 8168 DEGs were identified in CRC when the threshold was set at |fold change|> 2 and FDR < 0.05. Of these genes, 4496 were upregulated and 3672 were downregulated (Fig. 2A). Additionally, 326 HRGs and 288 LMRGs were retrieved from the MSigDB database. We then conducted univariate Cox regression analysis on all protein-coding genes in the RNA-seq data of patients and identified 2207 genes associated with overall survival (OS). After taking the intersection of DEGs, genes associated with OS, HRGs, and LMRGs, we identified 26 hypoxia-related DEGs (HRDEGs) and 9 lactate metabolism-related DEGs (LMRDEGs, Fig. 2B, C). Analysis of these 35 genes revealed a wide range of correlations and intergene interactions (Fig. 2D, E).
Constructing molecular subtypes based on the HRDEGs and LMRDEGs
Unsupervised clustering was performed on 575 CRC samples using the 35 identified HRDEGs and LMRDEGs. The highest average within-group agreement was observed when classified into three subtypes (Fig. 2F, G). Therefore, we classified the 575 samples into three subtypes: C1 (n = 92), C2 (n = 324), and C3 (n = 159). A comprehensive evaluation revealed that different molecular subtypes exhibited distinct hypoxia and lactate metabolism-related microenvironments (Figure S1A, B). The Kaplan–Meier survival curves indicated that patients with the C1 subtype had the worst OS than patients with subtypes C2 (HR = 1.67, 95% CI 1.07–2.62, p = 0.02) and C3 (HR = 2.48, 95% CI 1.41–4.36, p < 0.01). Additionally, there was a trend towards better OS in the C3 subtype compared with C2, although the difference was not statistically significant (C2 vs. C3, HR = 1.48, 95% CI 0.92–2.40, p = 0.11, Fig. 3A). The heat map in Fig. 3B shows the expression of the 35 genes in different subtypes and their association with clinical features. In particular, subtype C1 was associated with a higher number of patients with mucinous adenocarcinoma (MAC, p < 0.01) and advanced TNM stages (p < 0.01, Table 1). Additionally, we assessed the correlation between the molecular subtypes related to hypoxia and lactate metabolism and recognized CMS. Our findings revealed that the C1 subtype had a higher proportion of CMS4 (68.48%), whereas CMS2 dominated the C3 subtype (61.01%, Fig. 3C). GSEA of subtypes C1 and C3 revealed that the hypoxia pathway was significantly upregulated in C1 (FDR = 0.02). In contrast, the oxidative phosphorylation pathway was significantly downregulated (FDR = 0.01, Fig. 3D), consistent with our expectations. In addition, C1 exhibited a significant upregulation of KRAS signaling, EMT, and angiogenesis (Fig. 3E). In contrast, C3 showed a significant upregulation of MYC targets v2, MYC targets v1, and E2F targets (Fig. 3F).
Immune landscape of the different molecular subtypes
As hypoxia is closely related to the immune microenvironment, we evaluated immune cell infiltration in the three subtypes. We analyzed the infiltration of 22 immune cells using CIBERSORT and found that the C1 subtype predominantly contained macrophages, including M1 (p = 0.02) and M2 (p < 0.01). Meanwhile, an increase in the number of CD8 + T cells was observed in subtypes C2 and C3 (p = 0.01). Resting CD4 + memory T cells (p < 0.01), activated CD4 + memory T cells (p < 0.01), and activated dendritic cells (p < 0.01) were significantly more abundant in subtypes C2 and C3 than in the C1 subtype (Fig. 4A).
The cancer immune cycle (CIC) comprises seven major steps necessary for the immune-mediated control of tumor growth, beginning with the release of antigens from cancer cells and ending with the killing of cancer cells. The three subtypes were also evaluated for the seven major steps of CIC. The C1 subtype was found to play a significant role in the release of cancer antigens, priming and activation, trafficking of immune cells to tumors, and infiltration of immune cells into tumors (Fig. 4B). Additionally, Fig. 4C suggests that immune cells such as T cells, macrophages, and NK cells were widely recruited in the C1 subtype. However, the performance of the C1 subtype was slightly lower than that of the C2 and C3 subtypes in the recognition of cancer cells by T cells and killing of cancer cells, which are key steps in antitumor immunity, despite its ability to enrich more immune cells.
To evaluate the immunocompetence status of the three subtypes, we assessed their IPS. The IPS of the C1 subtype was significantly lower than those of subtypes C2 and C3 (p < 0.01, Fig. 4D), indicating that the C1 subtype was less responsive to immunotherapy. Specifically, the C1 subtype exhibited higher scores for MHC molecules (p < 0.01) and effector cells (p < 0.01) and lower scores for suppressor cells (p < 0.01) and checkpoints (p < 0.01, Fig. 4E). These findings appear to be inconsistent with the poor prognosis associated with the C1 subtype. To further investigate this phenomenon, we analyzed the immune functions of the samples using the TIDE algorithm. The results showed that the C1 subtype had a significantly higher overall TIDE score (p < 0.01), tumor immune dysfunction score (p < 0.01), and tumor immune exclusion score (p < 0.01) than the C2 and C3 subtypes (Fig. 4F), suggesting that immune escape was more common in the C1 subtype. Additionally, the PD-L1 score (p < 0.01, Fig. 4G) and cancer-associated fibroblast (CAF) score (p < 0.01, Fig. 4H) were significantly higher in the C1 subtype than in subtypes C2 and C3. We also evaluated the expression of 24 well-known immune checkpoint genes across the different subtypes. The results showed that these molecules, including CD274 (PD-L1), LAG3, CTLA4, and CD276, were upregulated in the C1 subtype (Fig. 4I).
Mutation mapping of the different molecular subtypes
To better understand the distinctions between the various subtypes, we examined them at the gene mutation level. Subtypes C1 and C2 exhibited a high total mutation burden (TMB), whereas subtype C3 exhibited a low TMB (p < 0.01, Fig. 5A–C). The top 20 genes in terms of mutation frequency were analyzed for each subtype. APC showed increasing mutation frequency in subtypes C1, C2, and C3 (67% vs. 72% vs. 82%), whereas MUC16 (37% vs. 30% vs. 13%) and PIK3CA (32% vs. 27% vs. 20%) showed decreasing mutation frequency in subtypes C1, C2, and C3, respectively. Additionally, TP53 mutation frequency was the highest in the C3 subtype (78%), followed by C1 (63%) and C2 (49%). The C2 subtype had the highest KRAS mutation frequency (48%), followed by subtypes C1 (40%) and C3 (33%).
Drug sensitivity
In clinical practice, stage II CRC patients with high-risk factors and patients with stage III and IV CRC require postoperative chemotherapy. Sensitivity to commonly used chemotherapeutic agents greatly affects the outcomes and prognosis of patients with CRC. Therefore, we predicted the sensitivity of patients with the three subtypes to the most commonly used chemotherapeutic agents for CRC, namely 5-FU, oxaliplatin, and irinotecan, based on the GDSC2 and CTRP2 databases. The sensitivities of the C1, C2, and C3 subtypes to 5-FU did not differ significantly in the GDSC2 database (p = 0.51, Fig. 5D). However, patients in the C1 subtype had a significantly higher half-maximal inhibitory concentration (IC50) for oxaliplatin than those in subtypes C2 (p < 0.01) and C3 (p < 0.01); no significant difference in oxaliplatin sensitivity was found between patients in subtypes C2 and C3 (p = 0.45, Fig. 5E). There was no significant difference in the IC50 for irinotecan between subtypes C1 and C2 (p = 0.39) or C3 (p = 0.22). Patients with the C3 subtype had slightly higher IC50 values than those with the C2 subtype (p < 0.01, Fig. 5F). When we switched to the CTRP2 database, we found that the AUC for 5-FU was higher in the C1 subtype than in the C2 (p < 0.01) and C3 (p < 0.01, Fig. 5G) subtypes, suggesting that the C1 subtype was the least sensitive to 5-FU. In addition, ssGSEA found higher stemness enrichment scores in C1 (Figure S2A) and stemness-associated genes such as LGR5 and CD34 were significantly overexpressed in the C1 subtype than in the other two subtypes (Figure S2B, C). The differences in characteristics of stemness may explain the differences in drug sensitivity between the subgroups.
To gain a better understanding of the sensitivity of tumor cells to chemotherapeutic drugs under hypoxic conditions, we performed in vitro experiments under hypoxic (1% O2) and normoxic (21% O2) conditions. The experimental workflow is shown in Fig. 6A. The sensitivity of the CRC cell lines HCT116 and LS174T to the conventional chemotherapeutic agents 5-FU, oxaliplatin, and irinotecan was characterized. Moreover, the cells were treated with or without the HIF-1 inhibitor BAY87-2243 under hypoxic and normoxic conditions. Trypan blue staining revealed that the survival rate of tumor cells was higher after 48 h of hypoxia compared with that under normoxic conditions at the same concentrations of 5-FU (7.5 µM), oxaliplatin (7.5 µM), or irinotecan (10 µM in HCT116 cells and 50 µM in LS174T cells) (Fig. 6B). These results were confirmed by CCK8 experiments (Fig. 6C–H). Notably, under hypoxic conditions, BAY87-2243 increased the sensitivity of HCT116 and LS174T cells to 5-FU to normoxic levels (Fig. 6C, F). However, under both normoxic and hypoxic conditions, BAY87-2243 decreased the sensitivity of HCT116 cells to oxaliplatin (Fig. 6D); this was also observed in LS174T cells under normoxic conditions (Fig. 6G). Under hypoxic conditions, BAY87-2243 increased the sensitivity of HCT116 and LS174T cells to irinotecan (Fig. 6E, H).
Construction of a prognostic prediction model
To clarify the molecular mechanisms underlying hypoxic tumors and their impact on the OS of patients with CRC, we compared the C1 subtype, which had the worst prognosis, with the C3 subtype, which had the best prognosis. Differential expression analysis using DESeq2 revealed that 1953 genes were upregulated in the C1 subtype, while 362 genes were downregulated (Fig. 7A). The heatmap in Fig. 7B displays the top 20 genes that were upregulated and downregulated. KEGG analysis of these 2315 DEGs revealed that they were predominantly enriched in cell adhesion molecules, ECM-receptor interactions, and the PI3K-Akt signaling pathway (Fig. 7C). GO functional enrichment analysis revealed that these DEGs were significantly enriched in extracellular structure organization, extracellular matrix organization, collagen-containing extracellular matrix, and extracellular matrix structural components, suggesting a close association between DEGs and extracellular matrix (Fig. 7D).
Univariate Cox regression analysis was performed to identify 2910 OS-related genes. After determining their intersection with the obtained 2315 DEGs, 830 OS-related DEGs were identified (Fig. 7E). Subsequently, the 830 DEGs were screened using a random forest model, and the top 20 genes were selected (Fig. 7F). To avoid problems with covariance and overfitting, we further screened these 20 DEGs using the LASSO method, which identified eight potential genes (Fig. 7G, H). Finally, we conducted a multivariate Cox regression analysis based on these eight genes and obtained six DEGs to construct a prognostic prediction signature, which we named the HLM score (Fig. 8A).
Validation of the predictive efficacy of the HLM score
We divided the 251 patients in the C1 and C3 subtypes into high-risk and low-risk groups based on their HLM scores. A marked decrease in the OS of CRC patients was observed as the HLM score increased (Fig. 8B). The Kaplan–Meier curves indicated that patients in the high-risk group had significantly worse OS (HR = 3.14, 95% CI 1.66–5.94, p < 0.01, Fig. 8C). The AUCs of the ROC for predicting 1-, 3-, and 4-year OS were 0.78 (95% CI 0.90–0.651), 0.80 (95% CI 0.90–0.71), and 0.70 (95% CI 0.83–0.57), respectively, suggesting superior predictive efficacy (Fig. 8D). Additionally, compared with clinical factors alone, we found that when we combined the HLM score with clinical factors such as age and T, N, and M stage, it significantly improved the AUCs (Figure S3). When the HLM score was applied to the entire TCGA CRC cohort (n = 575), the C1 subtype had the highest HLM score, the C3 subtype had the lowest HLM score, and the C2 subtype had an intermediate HLM score, consistent with our previous molecular subtyping results (Fig. 8E). The TCGA CRC cohort was also divided into high-risk and low-risk groups based on the HLM scores. GSEA of the high-risk and low-risk groups revealed that the hypoxia and angiogenesis pathway were significantly upregulated in the high-risk group. In contrast, the oxidative phosphorylation pathway was significantly downregulated (Figure S4A). The expression of HIF-1α and lactate metabolism-related genes such as IGFBP3, TGFB3, and CAV1 were significantly higher in the high-risk group than in the low-risk group (Figure S4B, D). Patients in the high-risk group had a significantly worse OS (HR = 1.64, 95% CI 1.16–2.38, p < 0.01, Fig. 8F). In addition, the high-risk group, similar to the C1 subtype, exhibited a higher proportion of CMS4 while the low-risk group exhibited a higher proportion of CMS2 (Figure S5A). The effect of CMS alone on CRC prognostic stratification was insignificant. However, the HLM score in concert with CMS demonstrated the capacity to refine CRC prognostic stratification (Figure S5B, C).
To validate the efficacy of the HLM score more extensively, we obtained three datasets from the GEO database. These datasets were categorized into high- and low-risk groups based on the HLM scores. Results showed that in the GSE106584 (HR = 2.10, 95% CI 1.25–3.53, p < 0.01, Fig. 8G), GSE17536 (HR = 1.94, 95% CI 1.14–3.31, p = 0.01, Fig. 8H), and GSE39582 (HR = 1.48, 95% CI 1.09–2.00, p = 0.02, Fig. 8I) datasets, patients in the high-risk group displayed inferior OS than patients in the low-risk group. Additionally, we analyzed the IMvigor210 immunotherapy cohort, in which patients in the high-risk group had a higher proportion of progressive disease (54.42% vs. 43.48%) and fewer patients with complete response/partial response/stable disease (30.61% vs. 40.79%, p = 0.02; Fig. 8J) than the low-risk group. The Kaplan–Meier curve also suggested worse OS for patients in the high-risk group in the IMvigor210 cohort than in the low-risk group (HR = 1.36, 95% CI 1.05–1.76, p = 0.02, Fig. 8K). Besides, we assessed the efficacy of the HLM score in these four validation sets through ROC curves, risk heatmaps, and calibration curves, all of which suggest that the HLM score is robust (Figure S6).
ssGSEA and scRNA-seq analysis to assess immune cell infiltration
The ssGSEA results indicated that the enrichment scores for CD8 + Tex, GZMK + Tex, terminal Tex, OXPHOS- Tex, and TCF7 + Tex were significantly higher in the C1 subtype than in subtypes C2 and C3 (p < 0.01, Fig. 9A). Similar results were obtained when TCGA CRC samples were categorized into high- and low-risk groups based on the HLM score (p < 0.01, Fig. 9B). To further validate these findings, publicly available scRNA-seq data were analyzed. After selection and filtering, 47,285 cells were included in the analysis. These cells clustered into five major types: B cells, T cells, myeloid cells, stromal cells, and epithelial cells (Fig. 9C, D). Based on the HLM score of each sample after pseudo-bulk analysis, the 23 CRC samples from the scRNA-seq data source were categorized into the high-risk (n = 11) and low-risk (n = 12) groups. The proportion of myeloid cells was slightly higher in the low-risk group than in the high-risk group (p = 0.05, Fig. 9E, F). We then performed subcluster annotation on 16,065T cells from these CRC samples (Fig. 9G). The proportion of exhausted CD8 + T cells (p = 0.05, Fig. 9H) and CD8 + intraepithelial lymphocytes (p = 0.04, Figure S7A) was higher in the high-risk group than in the low-risk group. To determine specific differences in the proportion of CD8 + T cells between the two groups, we performed further subpopulation annotation of CD8 + T cells (Fig. 9I). The results showed that the proportions of terminal Tex (p < 0.05, Fig. 9J), OXPHOS-Tex (p < 0.05, Fig. 9K) and GZMK + early Tem (effective memory T cells) (p < 0.05, Figure S7B) were significantly higher in the high-risk group than in the low-risk group. These findings suggest that HRGs and LMRGs may play a role in CD8 + T cell exhaustion in CRC.
In addition, we classified the 23 samples of scRNA-seq into C1/C2/C3 subtypes. It was found that, similar to the HLM score classification results, for the T cell subpopulation, exhausted CD8 + T cells especially OXPHOS- Tex were significantly higher in C1 than in C2 and C3 subtypes (Figure S7C, D). Furthermore, we found that naïve T (Tn) cells were lower in C1 than in C2, C3 subtypes (Figure S7D). Tn cells are precursors of effector and memory T cell subsets. Differences in Tn proportion between the three subtypes may be attributed to the fact that the proportion of exhausted T cells varies between the subtypes.
Establishment of a predictive nomogram
A predictive nomogram based on seven factors, namely age, sex, T stage, N stage, M stage, pathological type, and HLM score, was constructed to visualize the prognosis of patients in the total TCGA CRC cohort. The nomogram could predict the OS of patients with CRC with a C-index of 0.78 (Fig. 10A). Calibration curves for the nomogram also showed ideal prediction accuracy (Fig. 10B). The DCA curve demonstrated that the nomogram provided more net clinical benefits than the clinical characteristics alone (Fig. 10C). The nomogram based on the HLM score could predict both short- and long-term OS in patients with CRC and assist their clinical management.
Single gene analysis
To understand how the HLM score predicts the prognosis of patients with CRC, we performed a single gene analysis of six characterized genes in C1 and C3 subtypes, and the Kaplan–Meier curve indicated that the overexpression of TRIB2 or ELFN2 predicted poorer OS in patients with CRC, with HRs of 2.39 (95% CI 1.30–4.40, p < 0.01, Fig. 10D) and 1.88 (95% CI 1.05–3.38, p = 0.03, Fig. 10G), respectively. The expression matrix revealed that TRIB2 expression was significantly higher in CRC tissues than in normal tissues. In CRC, the expression of TRIB2 increased sequentially in subtypes C3, C2, and C1 (p < 0.01, Fig. 10E). Additionally, ELFN2 was highly expressed in CRC, particularly in the C1 subtype (p < 0.01, Fig. 10H). However, TRIB2 and ELFN2 are not included in the existing set of HRGs and LMRGs. To obtain information on the proteins interacting with TRIB2 and ELFN2, we consulted the BioPlex Interactome database. Our findings indicated that ISCA1, which interacts with TRIB2, is closely related to lactate metabolism (Fig. 10F). Similarly, KLHL24, which interacts with ELFN2, is closely associated with hypoxia (Fig. 10I). Therefore, we infer that TRIB2 and ELFN2 are associated with hypoxia and lactate metabolism.
Discussion
Hypoxia and acidosis within the TME, caused by the increased secretion of lactic acid and H+ as end products of glycolysis [49], are markers of solid tumors, including CRC. These conditions determine the selection of invasive and aggressive malignant clones that display resistance to radiotherapy, conventional chemotherapy, targeted therapy, or immunotherapy [50]. Due to the extensive intratumor heterogeneity and complex genetic and biological characteristics of CRC [51], predicting the outcome of patients with CRC remains challenging. To address this, we performed molecular subtyping of CRC based on HRGs and LMRGs and developed an HLM scoring system to predict prognosis and treatment response in patients with CRC.
Currently, molecular subtyping of CRC has shifted from a mutation-centered to a transcriptome-centered approach and from supervised to unsupervised clustering [52]. For instance, the CMS approach proposed by Guinney et al. in 2015, which is based on transcriptome analysis, has gained widespread acceptance [38]. To the best of our knowledge, this is the first study to propose an unsupervised clustering method for the molecular subtyping of CRC, based on HRGs and LMRGs. Our molecular subtyping method classified CRC into three subtypes and distinguished between their prognostic differences. In addition, the molecular subtypes generated correlated with the CMS and can complement it to guide the individualized treatment of patients with CRC.
Subsets of immune cells form a complicated network and cross-talk in different ways within the tumor [53]. Rather than their mere presence, the type of immune cells is more likely to be crucial. The C1 subtype predominantly contains macrophages, including M1 and M2 macrophages. Through specific differentiation, macrophages can evolve into two different polarization states in CRC: classically activated M1 (pro-inflammatory) and M2 (anti-inflammatory) macrophages. CRC originates from the epithelium; hence, TAMs are predominantly pro-inflammatory M1 macrophages in the early stages. However, they tend to convert to the anti-inflammatory and cancer-promoting M2 phenotype with tumorigenic activity as tumor cells utilize them to support the growth and progression of advanced CRC [54].
The results indicate that the hypoxic C1 subtype has abundant immune cell infiltration. However, the ability of these immune cells to kill tumor cells is inferior to that of the other two subtypes, suggesting that the number of infiltrating immune cells and the distribution of immune cell subsets play an important role in influencing antitumor immunity in the TME. TIDE analysis similarly demonstrated this, with the C1 subtype exhibiting higher tumor immune dysfunction and tumor immune exclusion scores than the other subtypes. This suggests that patients with the hypoxic subtype have a higher prevalence of immune escape from the TME, as evidenced by the upregulation of 24 well-known immune checkpoint genes. Additionally, the C1 subtype had a significantly higher CAF score than the C2 and C3 subtypes. CAFs play a crucial role in the reactive stroma of the TME and significantly influence tumor biology, including angiogenesis, invasion, immune evasion, metastasis, and drug resistance [55, 56]. Hypoxia activates resident fibroblasts by activating the reactive oxygen species (ROS) and HIF-1α pathways, transforming them into CAFs. CAFs undergo metabolic reprogramming to adapt to the TME and support glycolysis under hypoxic conditions. Increased glycolysis contributes to tumor cell proliferation [57, 58]. Furthermore, hypoxia affects chemokine and cytokine production by CAFs, leading to an increase in the proportion of immune cells associated with an immunosuppressive microenvironment [59].
Our study developed the HLM scoring system to predict the OS and efficacy of immunotherapy in patients with CRC using univariate Cox regression, the random forest model, LASSO, and multivariate Cox regression analyses. Notably, the group with a high HLM score exhibited a high proportion of CD8 + Tex infiltration, especially terminal Tex and OXPHOS CD8 + T cells. CD8 + T-cells are critical for eliminating cancer cells. We observed a reduction in both the number and function of CD8 + T cells in patients with high levels of hypoxia. In cancer, the continuous stimulation of antigen receptors can lead to an alternative differentiation trajectory known as T cell exhaustion, when antigens cannot be completely eliminated [60]. In CD8 + T cells, sustained exposure to hypoxia rapidly accelerates their differentiation to terminal exhaustion and represses antitumor immunity [61]. T cell exhaustion can be classified into four stages: T cell exhaustion progenitors 1, T cell exhaustion progenitors 2, T cell exhaustion intermediate and T cell exhaustion terminally [62]. Hypoxia results in abnormal OXPHOS in the mitochondria, which hinders T cell proliferation and promotes the transcriptional and metabolic reprogramming of the progenitor Tex (Tpex) and their differentiation into terminal Tex [63, 64]. Monoclonal antibody-based immunotherapies, such as anti-PD-L1/PD-1 agents, can become ineffective when T cells reach a stage of terminal exhaustion [62]. This explains why we observed that the hypoxic C1 subtype highly expresses numerous immune checkpoint molecules but still suffers from poor immunotherapeutic outcomes.
The HLM score was constructed based on six genes, AGXT, TRIB2, ELFN2, PCDHB10, CALCA, and CD79A, which have been previously reported to be associated with the tumorigenesis and progression of various cancer types, including CRC, as well as chemotherapy resistance [65,66,67,68,69]. When considering individual genes, we discovered that the high expression of TRIB2 and ELFN2 was associated with poor prognosis in patients with CRC. TRIB2 is a specialized member of the tribbles pseudokinase family that interacts with MAPK, AP4, CDC25, OCT3/4, C/EBP-alpha, ubiquitin E3 ligases, PCBP2, and AKT to regulate cellular processes, such as the cell cycle, senescence, stem cell pluripotency, protein degradation, and cell survival [70]. One study discovered that TRIB2 impedes tumor cell senescence, stimulates proliferation, and triggers cell cycle arrest via AP4/p21 signaling in CRC [66]. ELFN2 is a newly discovered hypomethylated gene that can inhibit the formation of protein phosphatase complexes and suppress the activity of protein phosphatase 1 by regulating it [71]. ELFN2 is also involved in glioblastoma cell autophagy [71], pancreatic cancer radiotherapy resistance [72], gastric cancer invasion [67], and endometrial cancer progression [73].
However, this study has certain limitations. First, most of our analyses and conclusions are based on data from public databases. Although our findings were validated across multiple datasets, it is important to acknowledge the possibility of bias. Moreover, we did not delve into the mechanisms by which hypoxia causes insensitivity to chemotherapeutic agents in patients with CRC. More in-depth in vitro and in vivo experiments are required to further explore the mechanisms underlying hypoxia-induced resistance to chemotherapy.
Conclusions
In conclusion, we molecularly typed patients with CRC based on their HRGs and LMRGs and revealed the immunologic and genetic characteristics of the different molecular subtypes. In addition, we constructed an HLM score that can be used to predict the prognosis and efficacy of immunotherapy of patients with CRC. This score has been validated in multiple datasets and has the potential to guide individualized and precise diagnosis and treatment of patients with CRC.
Availability of data and materials
Publicly available datasets were analyzed in this study. These data can be found in the TCGA (https://portal.gdc.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/), and IMvigor210CoreBiologies (http://research-pub.gene.com/IMvigor210CoreBiologies) databases.
Abbreviations
- AUC:
-
Area under the curve
- CAF:
-
Cancer-associated fibroblast
- CIC:
-
Cancer immune cycle
- CMS:
-
Consensus molecular subtypes
- CTRP:
-
Cancer Therapeutics Response Portal
- DCA:
-
Decision curve analysis
- DEG:
-
Differentially expressed gene
- EC:
-
Effector cell
- EMT:
-
Epithelial-mesenchymal transition
- GDSC:
-
Genomics of drug sensitivity in cancer
- GO:
-
Gene ontology
- GSEA:
-
Gene set enrichment analysis
- HIF:
-
Hypoxia-inducible factor
- HRG:
-
Hypoxia-related gene
- KEGG:
-
Kyoto Encyclopedia of Genes and Genomes
- LASSO:
-
Least absolute shrinkage and selection operator
- LMRG:
-
Lactate metabolism-related genes
- MAF:
-
Mutation annotation format
- MHC:
-
Major histocompatibility complex
- OS:
-
Overall survival
- ROC:
-
Receiver operating characteristic
- ROS:
-
Reactive oxygen species
- SC:
-
Suppressor cells
- scRNA-seq:
-
Single-cell RNA sequencing
- TAM:
-
Tumor-associated macrophage
- TCIA:
-
The cancer immunome atlas
- Tem:
-
Effective memory T cells
- Tex:
-
Exhausted T cells
- TIDE:
-
Tumor immune dysfunction and exclusion
- TMB:
-
Tumor mutation burden
- TPM:
-
Transcripts per million
- UMI:
-
Unique molecular identifier
References
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:209–49. https://doi.org/10.3322/caac.21660.
Ni T, He Z, Dai Y, Yao J, Guo Q, Wei L. Oroxylin A suppresses the development and growth of colorectal cancer through reprogram of HIF1α-modulated fatty acid metabolism. Cell Death Dis. 2017;8: e2865. https://doi.org/10.1038/cddis.2017.261.
Di Conza G, Trusso Cafarello S, Loroch S, Mennerich D, Deschoemaeker S, Di Matteo M, et al. The mTOR and PP2A pathways regulate PHD2 phosphorylation to fine-tune HIF1α levels and colorectal cancer cell survival under hypoxia. Cell Rep. 2017;18:1699–712. https://doi.org/10.1016/j.celrep.2017.01.051.
Simiantonaki N, Taxeidis M, Jayasinghe C, Kurzik-Dumke U, Kirkpatrick CJ. Hypoxia-inducible factor 1 alpha expression increases during colorectal carcinogenesis and tumor progression. BMC Cancer. 2008;8:320. https://doi.org/10.1186/1471-2407-8-320.
Li Q, Li Y, Li J, Ma Y, Dai W, Mo S, et al. FBW7 suppresses metastasis of colorectal cancer by inhibiting HIF1α/CEACAM5 functional axis. Int J Biol Sci. 2018;14:726–35. https://doi.org/10.7150/ijbs.24505.
Zhang W, Shi X, Peng Y, Wu M, Zhang P, Xie R, et al. HIF-1α promotes epithelial-mesenchymal transition and metastasis through direct regulation of ZEB1 in colorectal cancer. PLoS ONE. 2015;10: e0129603. https://doi.org/10.1371/journal.pone.0129603.
Ju S, Wang F, Wang Y, Ju S. CSN8 is a key regulator in hypoxia-induced epithelial-mesenchymal transition and dormancy of colorectal cancer cells. Mol Cancer. 2020;19:168. https://doi.org/10.1186/s12943-020-01285-4.
Li H, Rokavec M, Jiang L, Horst D, Hermeking H. Antagonistic effects of p53 and HIF1A on microRNA-34a regulation of PPP1R11 and STAT3 and hypoxia-induced epithelial to mesenchymal transition in colorectal cancer cells. Gastroenterology. 2017;153:505–20. https://doi.org/10.1053/j.gastro.2017.04.017.
Zhang J, Zhu L, Fang J, Ge Z, Li X. LRG1 modulates epithelial-mesenchymal transition and angiogenesis in colorectal cancer via HIF-1α activation. J Exp Clin Cancer Res. 2016;35:29. https://doi.org/10.1186/s13046-016-0306-2.
Sui H, Zhao J, Zhou L, Wen H, Deng W, Li C, et al. Tanshinone IIA inhibits β-catenin/VEGF-mediated angiogenesis by targeting TGF-β1 in normoxic and HIF-1α in hypoxic microenvironments in human colorectal cancer. Cancer Lett. 2017;403:86–97. https://doi.org/10.1016/j.canlet.2017.05.013.
Tang YA, Chen YF, Bao Y, Mahara S, Yatim SMJM, Oguz G, et al. Hypoxic tumor microenvironment activates GLI2 via HIF-1α and TGF-β2 to promote chemoresistance in colorectal cancer. Proc Natl Acad Sci USA. 2018;115:E5990–9. https://doi.org/10.1073/pnas.1801348115.
Dong S, Liang S, Cheng Z, Zhang X, Luo L, Li L, et al. ROS/PI3K/Akt and Wnt/β-catenin signalings activate HIF-1α-induced metabolic reprogramming to impart 5-fluorouracil resistance in colorectal cancer. J Exp Clin Cancer Res. 2022;41:15. https://doi.org/10.1186/s13046-021-02229-6.
Malier M, Gharzeddine K, Laverriere MH, Marsili S, Thomas F, Decaens T, et al. Hypoxia drives dihydropyrimidine dehydrogenase expression in macrophages and confers chemoresistance in colorectal cancer. Cancer Res. 2021;81:5963–76. https://doi.org/10.1158/0008-5472.CAN-21-1572.
Wei TT, Lin YT, Tang SP, Luo CK, Tsai CT, Shun CT, et al. Metabolic targeting of HIF-1α potentiates the therapeutic efficacy of oxaliplatin in colorectal cancer. Oncogene. 2020;39:414–27. https://doi.org/10.1038/s41388-019-0999-8.
Scharping NE, Menk AV, Whetstone RD, Zeng X, Delgoffe GM. Efficacy of PD-1 blockade is potentiated by metformin-induced reduction of tumor hypoxia. Cancer Immunol Res. 2017;5:9–16. https://doi.org/10.1158/2326-6066.CIR-16-0103.
Jayaprakash P, Ai M, Liu A, Budhani P, Bartkowiak T, Sheng J, et al. Targeted hypoxia reduction restores T cell infiltration and sensitizes prostate cancer to immunotherapy. J Clin Invest. 2018;128:5137–49. https://doi.org/10.1172/JCI96268.
Chamoto K, Chowdhury PS, Kumar A, Sonomura K, Matsuda F, Fagarasan S, et al. Mitochondrial activation chemicals synergize with surface receptor PD-1 blockade for T cell-dependent antitumor activity. Proc Natl Acad Sci U S A. 2017;114:E761–70. https://doi.org/10.1073/pnas.1620433114.
Najjar YG, Menk AV, Sander C, Rao U, Karunamurthy A, Bhatia R, et al. Tumor cell oxidative metabolism as a barrier to PD-1 blockade immunotherapy in melanoma. JCI Insight. 2019. https://doi.org/10.1172/jci.insight.124989.
Lequeux A, Noman MZ, Xiao M, Van Moer K, Hasmim M, Benoit A, et al. Targeting HIF-1 alpha transcriptional activity drives cytotoxic immune effector cells into melanoma and improves combination immunotherapy. Oncogene. 2021;40:4725–35. https://doi.org/10.1038/s41388-021-01846-x.
Wang Y, Yin K, Tian J, Xia X, Ma J, Tang X, et al. Granulocytic myeloid-derived suppressor cells promote the stemness of colorectal cancer cells through exosomal S100A9. Adv Sci. 2019;6:1901278. https://doi.org/10.1002/advs.201901278.
Corzo CA, Condamine T, Lu L, Cotter MJ, Youn JI, Cheng P, et al. HIF-1α regulates function and differentiation of myeloid-derived suppressor cells in the tumor microenvironment. J Exp Med. 2010;207:2439–53. https://doi.org/10.1084/jem.20100587.
Chiu DK, Tse AP, Xu IM, Di Cui J, Lai RK, Li LL, et al. Hypoxia inducible factor HIF-1 promotes myeloid-derived suppressor cells accumulation through ENTPD2/CD39L1 in hepatocellular carcinoma. Nat Commun. 2017;8:517. https://doi.org/10.1038/s41467-017-00530-7.
Kumar V, Cheng P, Condamine T, Mony S, Languino LR, McCaffrey JC, et al. CD45 phosphatase inhibits STAT3 transcription factor activity in myeloid cells and promotes tumor-associated macrophage differentiation. Immunity. 2016;44:303–15. https://doi.org/10.1016/j.immuni.2016.01.014.
Noman MZ, Janji B, Hu S, Wu JC, Martelli F, Bronte V, et al. Tumor-promoting effects of myeloid-derived suppressor cells are potentiated by hypoxia-induced expression of miR-210. Cancer Res. 2015;75:3771–87. https://doi.org/10.1158/0008-5472.CAN-15-0405.
Clambey ET, McNamee EN, Westrich JA, Glover LE, Campbell EL, Jedlicka P, et al. Hypoxia-inducible factor-1 alpha-dependent induction of FoxP3 drives regulatory T-cell abundance and function during inflammatory hypoxia of the mucosa. Proc Natl Acad Sci USA. 2012;109:E2784–93. https://doi.org/10.1073/pnas.1202366109.
Lunt SY, Vander Heiden MG. Aerobic glycolysis: meeting the metabolic requirements of cell proliferation. Annu Rev Cell Dev Biol. 2011;27:441–64. https://doi.org/10.1146/annurev-cellbio-092910-154237.
Végran F, Boidot R, Michiels C, Sonveaux P, Feron O. Lactate influx through the endothelial cell monocarboxylate transporter MCT1 supports an NF-κB/IL-8 pathway that drives tumor angiogenesis. Cancer Res. 2011;71:2550–60. https://doi.org/10.1158/0008-5472.CAN-10-2828.
Sonveaux P, Copetti T, De Saedeleer CJ, Végran F, Verrax J, Kennedy KM, et al. Targeting the lactate transporter MCT1 in endothelial cells inhibits lactate-induced HIF-1 activation and tumor angiogenesis. PLoS ONE. 2012;7: e33418. https://doi.org/10.1371/journal.pone.0033418.
Lee DC, Sohn HA, Park ZY, Oh S, Kang YK, Lee KM, et al. A lactate-induced response to hypoxia. Cell. 2015;161:595–609. https://doi.org/10.1016/j.cell.2015.03.011.
Colegio OR, Chu NQ, Szabo AL, Chu T, Rhebergen AM, Jairam V, et al. Functional polarization of tumour-associated macrophages by tumour-derived lactic acid. Nature. 2014;513:559–63. https://doi.org/10.1038/nature13490.
Brand A, Singer K, Koehl GE, Kolitzus M, Schoenhammer G, Thiel A, et al. LDHA-associated lactic acid production blunts tumor immunosurveillance by T and NK cells. Cell Metab. 2016;24:657–71. https://doi.org/10.1016/j.cmet.2016.08.011.
Haas R, Smith J, Rocher-Ros V, Nadkarni S, Montero-Melendez T, D’Acquisto F, et al. Lactate regulates metabolic and pro-inflammatory circuits in control of T cell migration and effector functions. PLOS Biol. 2015;13: e1002202. https://doi.org/10.1371/journal.pbio.1002202.
Kumagai S, Koyama S, Itahashi K, Tanegashima T, Lin YT, Togashi Y, et al. Lactic acid promotes PD-1 expression in regulatory T cells in highly glycolytic tumor microenvironments. Cancer Cell. 2022;40:201-218.e9. https://doi.org/10.1016/j.ccell.2022.01.001.
Watson MJ, Vignali PDA, Mullett SJ, Overacre-Delgoffe AE, Peralta RM, Grebinoski S, et al. Metabolic support of tumour-infiltrating regulatory T cells by lactic acid. Nature. 2021;591:645–51. https://doi.org/10.1038/s41586-020-03045-2.
Wei Y, Xu H, Dai J, Peng J, Wang W, Xia L, et al. Prognostic significance of serum lactic acid, lactate dehydrogenase, and albumin levels in patients with metastatic colorectal cancer. BioMed Res Int. 2018;2018:1804086. https://doi.org/10.1155/2018/1804086.
Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P, et al. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 2010;38:W214–20. https://doi.org/10.1093/nar/gkq537.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554:544–8. https://doi.org/10.1038/nature25501.
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. https://doi.org/10.1038/nm.3967.
Zeng D, Ye Z, Shen R, Yu G, Wu J, Xiong Y, et al. IOBR: multi-omics immuno-oncology biological research to decode tumor microenvironment and signatures. Front Immunol. 2021;12:687975. https://doi.org/10.3389/fimmu.2021.687975.
Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, et al. TIP: a web server for resolving tumor immunophenotype profiling. Cancer Res. 2018;78:6575–80. https://doi.org/10.1158/0008-5472.CAN-18-0689.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24:1550–8. https://doi.org/10.1038/s41591-018-0136-1.
Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021. https://doi.org/10.1093/bib/bbab260.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747–56. https://doi.org/10.1101/gr.239244.118.
Zhang L, Yu X, Zheng L, Zhang Y, Li Y, Fang Q, et al. Lineage tracking reveals dynamic relationships of T cells in colorectal cancer. Nature. 2018;564:268–72. https://doi.org/10.1038/s41586-018-0694-x.
Zheng L, Qin S, Si W, Wang A, Xing B, Gao R, et al. Pan-cancer single-cell landscape of tumor-infiltrating T cells. Science. 2021;374: abe6474. https://doi.org/10.1126/science.abe6474.
Lee HO, Hong Y, Etlioglu HE, Cho YB, Pomella V, Van den Bosch B, et al. Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nat Genet. 2020;52:594–603. https://doi.org/10.1038/s41588-020-0636-z.
Hao Y, Stuart T, Kowalski MH, Choudhary S, Hoffman P, Hartman A, et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol. 2024;42:293–304. https://doi.org/10.1038/s41587-023-01767-y.
Huttlin EL, Bruckner RJ, Navarrete-Perea J, Cannon JR, Baltier K, Gebreab F, et al. Dual proteome-scale networks reveal cell-specific remodeling of the human interactome. Cell. 2021;184:3022–40. https://doi.org/10.1016/j.cell.2021.04.011.
Kung-Chun Chiu D, Pui-Wah Tse A, Law CT, Ming-Jing XuI, Lee D, Chen M, et al. Hypoxia regulates the mitochondrial activity of hepatocellular carcinoma cells through HIF/HEY1/PINK1 pathway. Cell Death Dis. 2019;10:934. https://doi.org/10.1038/s41419-019-2155-3.
Kopecka J, Salaroglio IC, Perez-Ruiz E, Sarmento-Ribeiro AB, Saponara S, De Las RJ, et al. Hypoxia as a driver of resistance to immunotherapy. Drug Resist Updat. 2021;59:100787. https://doi.org/10.1016/j.drup.2021.100787.
Li J, Ma X, Chakravarti D, Shalapour S, DePinho RA. Genetic and biological hallmarks of colorectal cancer. Genes Dev. 2021;35:787–820. https://doi.org/10.1101/gad.348226.120.
Wang W, Kandimalla R, Huang H, Zhu L, Li Y, Gao F, et al. Molecular subtyping of colorectal cancer: recent progress, new challenges and emerging opportunities. Semin Cancer Biol. 2019;55:37–52. https://doi.org/10.1016/j.semcancer.2018.05.002.
Picard E, Verschoor CP, Ma GW, Pawelec G. Relationships between immune landscapes, genetic subtypes and responses to immunotherapy in colorectal cancer. Front Immunol. 2020;11:369. https://doi.org/10.3389/fimmu.2020.00369.
Wang H, Tian T, Zhang J. Tumor-associated macrophages (TAMs) in colorectal cancer (CRC): from mechanism to therapy and prognosis. Int J Mol Sci. 2021. https://doi.org/10.3390/ijms22168470.
Yang D, Liu J, Qian H, Zhuang Q. Cancer-associated fibroblasts: From basic science to anticancer therapy. Exp Mol Med. 2023;55:1322–32. https://doi.org/10.1038/s12276-023-01013-0.
Bigos KJ, Quiles CG, Lunj S, Smith DJ, Krause M, Troost EG, et al. Tumour response to hypoxia: understanding the hypoxic tumour microenvironment to improve treatment outcome in solid tumours. Front Oncol. 2024;14:1331355. https://doi.org/10.3389/fonc.2024.1331355.
Zhang D, Wang Y, Shi Z, Liu J, Sun P, Hou X, et al. Metabolic reprogramming of cancer-associated fibroblasts by IDH3α downregulation. Cell Rep. 2015;10:1335–48. https://doi.org/10.1016/j.celrep.2015.02.006.
Becker LM, O’Connell JT, Vo AP, Cain MP, Tampe D, Bizarro L, et al. Epigenetic reprogramming of cancer-associated fibroblasts deregulates glucose metabolism and facilitates progression of breast cancer. Cell Rep. 2020;31:107701. https://doi.org/10.1016/j.celrep.2020.107701.
Pei L, Liu Y, Liu L, Gao S, Gao X, Feng Y, et al. Roles of cancer-associated fibroblasts (CAFs) in anti- PD-1/PD-L1 immunotherapy for solid cancers. Mol Cancer. 2023;22:29. https://doi.org/10.1186/s12943-023-01731-z.
McLane LM, Abdel-Hakeem MS, Wherry EJ. CD8 T cell exhaustion during chronic viral infection and cancer. Annu Rev Immunol. 2019;37:457–95. https://doi.org/10.1146/annurev-immunol-041015-055318.
Scharping NE, Rivadeneira DB, Menk AV, Vignali PDA, Ford BR, Rittenhouse NL, et al. Mitochondrial stress induced by continuous stimulation under hypoxia rapidly drives T cell exhaustion. Nat Immunol. 2021;22:205–15. https://doi.org/10.1038/s41590-020-00834-9.
Beltra JC, Manne S, Abdel-Hakeem MS, Kurachi M, Giles JR, Chen Z, et al. Developmental relationships of four exhausted CD8+ T cell subsets reveals underlying transcriptional and epigenetic landscape control mechanisms. Immunity. 2020;52:825-41.e8. https://doi.org/10.1016/j.immuni.2020.04.014.
Wu H, Zhao X, Hochrein SM, Eckstein M, Gubert GF, Knöpper K, et al. Mitochondrial dysfunction promotes the transition of precursor to terminally exhausted T cells through HIF-1α-mediated glycolytic reprogramming. Nat Commun. 2023;14:6858. https://doi.org/10.1038/s41467-023-42634-3.
Mao Y, Zhang J, Zhou Q, He X, Zheng Z, Wei Y, et al. Hypoxia induces mitochondrial protein lactylation to limit oxidative phosphorylation. Cell Res. 2024;34:13–30. https://doi.org/10.1038/s41422-023-00864-6.
Kjersem JB, Thomsen M, Guren T, Hamfjord J, Carlsson G, Gustavsson B, et al. AGXT and ERCC2 polymorphisms are associated with clinical outcome in metastatic colorectal cancer patients treated with 5-FU/oxaliplatin. Pharmacogenomics J. 2016;16:272–9. https://doi.org/10.1038/tpj.2015.54.
Hou Z, Guo K, Sun X, Hu F, Chen Q, Luo X, et al. TRIB2 functions as novel oncogene in colorectal cancer by blocking cellular senescence through AP4/p21 signaling. Mol Cancer. 2018;17:172. https://doi.org/10.1186/s12943-018-0922-x.
Zhang YQ, Zhang JJ, Song HJ, Li DW. Overexpression of CST4 promotes gastric cancer aggressiveness by activating the ELFN2 signaling pathway. Am J Cancer Res. 2017;7:2290–304.
Nabbi A, Danesh A, Espin-Garcia O, Pedersen S, Wellum J, Fu LH, et al. Multimodal immunogenomic biomarker analysis of tumors from pediatric patients enrolled to a phase 1–2 study of single-agent atezolizumab. Nat Cancer. 2023;4:502–15. https://doi.org/10.1038/s43018-023-00534-x.
Sorrentino C, D’Antonio L, Fieni C, Ciummo SL, Di Carlo E. Colorectal cancer-associated immune exhaustion involves T and B lymphocytes and conventional NK cells and correlates with a shorter overall survival. Front Immunol. 2021;12:778329. https://doi.org/10.3389/fimmu.2021.778329.
Mayoral-Varo V, Jiménez L, Link W. The critical role of TRIB2 in cancer and therapy resistance. Cancers. 2021. https://doi.org/10.3390/cancers13112701.
Liu C, Fu H, Liu X, Lei Q, Zhang Y, She X, et al. LINC00470 coordinates the epigenetic regulation of ELFN2 to distract GBM cell autophagy. Mol Ther. 2018;26:2267–81. https://doi.org/10.1016/j.ymthe.2018.06.019.
Souchek JJ, Baine MJ, Lin C, Rachagani S, Gupta S, Kaur S, et al. Unbiased analysis of pancreatic cancer radiation resistance reveals cholesterol biosynthesis as a novel target for radiosensitisation. Br J Cancer. 2014;111:1139–49. https://doi.org/10.1038/bjc.2014.385.
Dong Y, Zhang T, Li X, Yu F, Yu H, Shao S. Identification of key prognostic-related miRNA-mRNA pairs in the progression of endometrial carcinoma. Gynecol Obstet Investig. 2022;87:12–21. https://doi.org/10.1159/000520339.
Acknowledgements
The authors thank Yaqi Wang of the School of Life Sciences, Peking University, for providing revisions to this study.
Funding
This study was funded by the National Natural Science Foundation of China (82073223).
Author information
Authors and Affiliations
Contributions
Conceptualization: JG. Data curation: AH. Formal analysis: AH and ZS. Funding acquisition: JG. Methodology: AH. Software: HH. Supervision: JG. Validation: ZS. Writing the original draft: AH and ZS. Writing, review and editing: YY, JC, and ZG. All the authors contributed to the manuscript and approved the submitted version.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
12967_2024_5391_MOESM1_ESM.docx
Supplementary Material 1. Figure S1. Hypoxia- and lactate metabolism-related microenvironments of the three subtypes. A Hypoxia microenvironment. B Lactate metabolism-related microenvironment. Figure S2. Stemness enrichment scores and expression of stemness-associated genes in different subtypes. A Stemness enrichment scores. B Exprission of LGR5. C Exprission of CD34. Figure S3. Time-dependent ROC curves for predicting 1-, 3-, and 5-year OS. A Clinical factors alone. B The HLM score combined with clinical factors. Figure S4. Hypoxia- and lactate metabolism-related microenvironments of the high-risk and low-risk groups. A GSEA of the high-risk and low-risk groups. B Hypoxia- and lactate metabolism-related microenvironment of the high-risk and low-risk groups. Figure S5. Comparison of HLM Score and CMS. A Distribution of the CMS in the high-risk and low-risk group. B CMS for CRC prognostic stratification. C HLM score combined with CMS for CRC prognostic stratification. Figure S6. ROC curves, risk heatmaps, and calibration curves in the validation sets. A GSE106584, B GSE17536, C GSE39582 and D IMvigor210. Figure S7. scRNA-seq analysis to assess immune cell infiltration. A Percentage of T cell subpopulation infiltration in high-risk vs. low-risk groups. B Percentage of CD8+ T cell subpopulation infiltration in high-risk vs. low-risk groups. C Percentage of T cell subpopulation infiltration in different molecular subtypes. D Percentage of CD8+ T cell subpopulation infiltration in different molecular subtypes. Tem, effective memory T cells; IEL, intraepithelial lymphocyte; Tc17, IL-17-producing CD8+ T cells; Tm, memory T cells; Tn, naïve T cells; Trm, tissue-resident memory T cells.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Huang, A., Sun, Z., Hong, H. et al. Novel hypoxia- and lactate metabolism-related molecular subtyping and prognostic signature for colorectal cancer. J Transl Med 22, 587 (2024). https://doi.org/10.1186/s12967-024-05391-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12967-024-05391-5