Deciphering lactate metabolism in colorectal cancer: Prognostic modeling, immune infiltration, and gene mutation insights

Abstract

BACKGROUND

Colorectal cancer (CRC) remains a major global health burden due to its high incidence and mortality, with treatment efficacy often hindered by tumor heterogeneity, drug resistance, and a complex tumor microenvironment (TME). Lactate metabolism plays a pivotal role in reshaping the TME, promoting immune evasion and epithelial-mesenchymal transition, making it a promising target for novel therapeutic strategies and prognostic modeling in CRC.

AIM

To offer an in-depth analysis of the role of lactate metabolism in CRC, highlighting its significance in the TME and therapeutic response.

METHODS

Utilizing single-cell and transcriptomic data from the Gene Expression Omnibus and The Cancer Genome Atlas, we identified key lactate metabolic activities, particularly in the monocyte/macrophage subpopulation.

RESULTS

Seven lactate metabolism-associated genes were significantly linked to CRC prognosis and used to construct a predictive model. This model accurately forecasts patient outcomes and reveals notable distinct patterns of immune infiltration and transcriptomic profiles mutation profiles between high- and low-risk groups. High-risk patients demonstrated elevated immune cell infiltration, increased mutation frequencies, and heightened sensitivity to specific drugs (AZD6482, tozasertib, and SB216763), providing a foundation for personalized treatment approaches. Additionally, a nomogram integrating clinical and metabolic data effectively predicted 1-, 3-, and 5-year survival rates.

CONCLUSION

This report underscored the pivotal mechanism of lactate metabolism in CRC prognosis and suggest novel avenues for therapeutic intervention.

Key Words: Colorectal cancer; Lactate metabolism; Prognostic model; Immune infiltration; Gene mutation analysis

Core Tip: This study presents a comprehensive analysis of lactate metabolism in colorectal cancer, identifying its pivotal role in shaping the tumor microenvironment and influencing immune infiltration. A predictive model utilizing a panel of seven lactate metabolism-related genes was developed, accurately predicting patient outcomes and drug sensitivity. Integration of transcriptomic and single-cell data highlights monocyte/macrophage metabolic activity, offering novel insights into personalized colorectal cancer therapy and the potential of targeting lactate metabolism to enhance treatment efficacy.

INTRODUCTION

Globally, colorectal cancer (CRC) stands out as a major contributor to cancer-related disease burden, affecting millions annually. Its incidence is notably higher in developed countries, with lifestyle and dietary factors contributing to the rising case numbers. Despite continuous advances in treatment, CRC remains a leading cause of cancer-related mortality, highlighting the need for early detection and improved therapeutic strategies[1,2]. The treatment of CRC involves multimodal approaches, including surgery, chemotherapy, radiotherapy, and targeted therapy. Treatment options depend on the stage of the disease, with surgery being the primary option for localized tumors. Chemotherapy and radiotherapy are used in advanced stages and as adjuvant therapies to reduce the risk of recurrence[3]. Targeted therapies and immunotherapies have shown potential, especially for metastatic CRC, to attack cancer cells by targeting specific molecular pathways and enhancing the immune response[1]. However, these treatments are usually only beneficial for patients with specific genetic mutations.

One of the key challenges in CRC management is drug resistance, which often lead to treatment failure and disease progression[4]. Developing effective therapies is further hindered by the highly variable tumor microenvironment (TME) and the diverse pathological features of CRC[5,6]. In addition, metastatic disease, where cancer spreads beyond the primary site, presents significant therapeutic challenges, as exemplified by the aggressive nature of advanced CRC[7]. Recent studies have focused on overcoming these challenges by exploring new therapeutic targets, improving the effectiveness of existing treatments, and developing more precise prognostic models to guide therapy[6,7]. The integration of genomic data with other data types and advances in precision medicine are expected to improve CRC prognosis by enabling more personalized and effective treatment options[8].

Tumor cells preferentially undergo glycolysis to produce energy through lactate metabolism, aerobic glycolysis persists in tumor cells, a phenomenon referred to as the Warburg effect[9,10]. Lactate serves not only as an energy source but also plays a critical role in promoting tumor invasiveness by acidifying the microenvironment and suppressing immune cell function, thereby facilitating tumor evasion of immune surveillance[11]. Understanding the role of lactate metabolism in tumor development could lead to new therapeutic strategies; such as inhibiting glycolytic pathways or modulating the acidity of the TME, ultimately improving treatment efficacy for CRC patients[12].

Lactate metabolism is closely associated with epithelial-mesenchymal transition activity. Epithelial-mesenchymal transition promotes more aggressive cellular activity, including metastasis[13]. Thus, an in-depth understanding of the mechanisms underlying lactate metabolism in CRC could uncover novel approaches to improve patient outcomes. We systematically investigate the role of lactate metabolism in CRC to gain mechanistic insights by examining its role in the TME, characterizing metabolism-related gene signatures implicated in lactate activity, followed by developing a prognostic framework for CRC outcome estimation. These findings could offer valuable insights for developing targeted therapeutic strategies and improving CRC prognosis.

MATERIALS AND METHODS

Single-cell data sources and data analysis

The dataset GSE200997, containing single-cell RNA sequencing data for CRC (sequencing platform: GPL21697), which included information on 16 tumor samples and 7 normal samples, was acquired via the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) database. The Seurat R package[14] was used for single-cell data analysis. First, to ensure data quality, cells with the following features were excluded from the analysis: Cells with fewer than 300 total genes detected; cells with 5000 or more total genes detected; cells with 100000 or more RNA molecules; cells with a mitochondrial gene percentage of 20% or more; and cells with an erythrocyte gene percentage of 50% or more. Genes detected in fewer than 3 cells were excluded from the analysis. The data were normalized using the NormalizeData function to ensure comparable data across samples. The top 2000 most variable genes were selected by the FindVariableFeatures function, and these highly variable genes were linearly transformed using ScaleData to eliminate the possible effects of different levels of data. Principal component analysis (PCA) was performed on these genes via the RunPCA method to reduce the dimensionality of the data. The harmony package[15] was used for batch correction between samples to ensure the consistency of the data from different samples. Based on the PCA fragmentation map, we identified the first 28 principal components for cluster analysis to reveal the intrinsic structure among cell populations. Leiden’s algorithm was used to cluster single cells and further define subpopulations of cells. The uniform manifold approximation and projection method was used for dimensionality reduction and visualization to make the differences between cell populations more intuitive. Finally, the cells were annotated with the SingleR package[16] according to the expression of specific marker genes.

CRC transcriptome data sources

Sequencing data for colon and rectal adenocarcinomas (COADREAD) were sourced from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/); this dataset contains information on 647 tumor samples and 51 normal samples. Sample data in TCGA-COADREAD were preprocessed, where normal samples and tumor patients with survival times of less than 30 days were removed, 618 samples were retained. Moreover, the sequencing data of CRC patients with the registration number GSE103479 (sequencing platform: GPL23985) were acquired from the Gene Expression Omnibus database. Tumor data derived from GSE103479 were also screened using the same data preprocessing standards, 156 tumor samples were retained for subsequent analysis. When analyzing the transcriptome data, model development was based on the TCGA dataset, with model validation was carried out externally using the GSE103479 sample set.

Analysis of lactate metabolism

Twenty-five genes related to lactate metabolism were obtained from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb). Eleven cellular subpopulations were identified in the single-cell analysis. By analyzing transcriptional patterns of 25 genes involved in lactate metabolism, the AUCell package[17] was selected to calculate the lactate metabolism activity in 11 cell subpopulations. The scMetabolism package[18] was used to assess lactate metabolism activity and glycolysis/gluconeogenesis pathway activity in 11 cell subpopulations at single-cell resolution. In addition, to analyze the strength of cellular interactions between cellular subpopulations, cell communication analysis was performed based on the CellChat package[19].

Derivation of a prognostic gene signature based on lactate metabolism characteristics

The edgeR package[20] was used to identify genes related to lactate metabolism that were differentially expressed in the monocyte/macrophage subpopulation; these genes were stratified into high- and low-activity subgroups according to median values of lactate metabolism scores [P < 0.05, |log2 fold change| > log2(1.5)]. Subsequently, based on univariate Cox analysis, LASSO Cox analysis and multivariate Cox analysis, lactate metabolism genes significantly associated with CRC prognosis were identified. A lactate metabolism-driven prognostic signature for CRC was formulated through integration of transcriptomic data and regression-based weighting. Patients were stratified into distinct prognostic categories based on the median model-derived risk score, and heatmaps illustrating expression profiles of selected genes were generated using the pheatmap R package[21], visualizations of risk score patterns and corresponding survival outcomes were generated. The genes used to construct the risk model were subjected to PCA via the FactoMineR[22] and factoextra packages[23] to assess the effects of each feature gene on the classification of patients into high- and low-risk score groups. The survival package in R was applied to conduct Kaplan-Meier survival analysis[24], and survival curves were plotted using the survminer package. Using the timeROC algorithm, we evaluated survival prediction performance at multiple time points, and validated the area under the curve (AUC) metrics in the external GSE103479 cohort[25].

Analysis of immune infiltration

Immune gene sets defined in the Thorsson et al[26] reference were used to characterize immune cell types. To assess TME heterogeneity, single-sample gene set enrichment analysis was applied to compare immune pathway enrichment between stratified risk groups. We also assessed TME activity and immune cell infiltration abundance in CRC patients according to the estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE)[27] and cell-type identification by estimating relative subsets of RNA transcript algorithms[28], respectively. Immunopheno score (IPS), reflecting immune response potential, were collected from The Cancer Imaging Archive (https://tcia.at/home) and applied to assess therapeutic variation between stratified risk groups.

Mutation characterization

The copy number variation profiles from the TCGA-COADREAD sequencing project were accessed via the TCGA repository for analysis. Graphical representation of the variation positions of high-frequency mutated genes based on the functionality of the maftools package[29] was generated to facilitate the exploration of the phenomena of comutation and mutual exclusivity between pairs of these genes.

Analysis of drug sensitivity

Information on drug sensitivity and molecular markers of drug response in cancer cells was downloaded from the Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/) database. Based on the functionality of the OncoPredict package[30], we downloaded Genomics of Drug Sensitivity in Cancer 2 gene expression profiles and corresponding drug response information to generate ridge regression models that could be applied to CRC transcriptome data and then generated sensitivity scores to predict half-maximal inhibitory concentrations for all drugs in the database in CRC patients.

Nomogram for predicting CRC patient survival

Independent prognostic indicators for CRC were determined through univariate and multivariate Cox regression based on variables including age, sex, tumor-node-metastasis stage, and model-derived risk score, using data from the TCGA-COADREAD cohort. Then, based on the identified independent clinical prognostic factors and risk score, nomogram models were developed to estimate overall survival probabilities at 1, 3, and 5 years for CRC patients. In addition, we plotted 1-, 3-, and 5-year calibration curves and decision curves to assess the accuracy of the nomogram models.

Western blot

Five pairs of CRC samples and paracancerous tissue samples were collected for the trial, as well as the procurement of normal colon cells lines (FHC, Sands Circle, VA, United States), and CRC cell lines SW620 (Sunn Bio, Shanghai, China), HCT116 (Sunn Bio, Shanghai, China), and WS480 (Shanghai Cell Bank, Shanghai, China). Total protein was extracted with radioimmunoprecipitation assay buffer. After protein quantification, aliquots of protein samples were loaded and separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and then transferred to 0.45-μm polyvinylidene fluoride membranes. 5% skimmed milk was blocked for 2 hours. The corresponding primary antibodies were incubated at 4 °C overnight. Primary antibodies were: Anti-leptin receptor overlapping transcript-like 1 (anti-LEPROTL1) antibody (TA331958, Origene, MD, United States), anti-inducible T cell costimulator (anti-ICOS) antibody (AP13108PU-N, Origene, MD, United States), anti-RNA-binding protein 3 antibody (EPR6061, Abcam, United Kingdom), anti-basic leucine zipper transcription factor ATF-like (BATF) antibody (EPR21911-55, Abcam, United Kingdom), anti-solute carrier family 2 member 3 (anti-SLC2A3) antibody (LS-B7024-50, Lifespan Biosciences, WA, United States), anti-acyl-CoA dehydrogenase very-long-chain (anti-ACADVL) antibody (AP32801PU-N, Origene, MD, United States), anti-mitochondrial-dihydronicotinamide adenine dinucleotide dehydrogenase subunit 2 (anti-MT-ND2) antibody (LS-B15632-50, Lifespan Biosciences, WA, United States), anti-glyceraldehyde-3-phosphate dehydrogenase antibody (LS B4075, Lifespan Biosciences, WA, United States). Subsequently, the membranes were incubated with secondary antibodies for 1 hour at room temperature. Immunoblots were displayed with an electrochemiluminescence fluorescence detection kit (Beyotime, Shanghai, China) and visualized with a Tanon 4600 system (Tanon Science and Technology Co., Ltd., Shanghai, China).

Immunofluorescence staining and quantification of fluorescence intensity

The collected CRC tissue and paracancerous tissue cells were fixed in phosphate buffered saline (PBS) containing 4% formaldehyde for 20 minutes at room temperature. After fixation, the samples were washed three times with PBS containing 1% bovine serum albumin (BSA). After fixation, samples were permeabilized with 1% BSA/PBS solution containing 0.5% Triton X-100 for 20 minutes. After permeabilization, the samples were again washed three times with 1% BSA/PBS. The samples were closed for 30 minutes using 1% BSA/PBS solution containing 0.1% Triton X-100. Primary antibodies were incubated overnight at 4 °C. After primary antibody incubation, samples were washed and incubated with secondary antibody. After secondary antibody incubation, samples were washed three times with 1% BSA/PBS. Nuclear staining was then performed using 4,6-diamidino-2-phenylindole for 10 minutes to help visualize nuclei. After secondary antibody incubation, samples were washed three times with 1% BSA/PBS. Nuclear staining was then performed using 4,6-diamidino-2-phenylindole for 10 minutes to help visualize the nuclei. The stained samples were made into slides and encapsulated using appropriate anti-fade encapsulation media to maintain fluorescence. Immunofluorescence images were captured using a LSM710 laser scanning confocal microscope (Carl Zeiss, Germany).

Antibodies

Primary antibodies and dilutions: Anti-Pan-Kla (PTM-1401, PTM Bio, Hangzhou, China) = 1:500, anti-histone H3 lysine 18 lactylation (H3K18 La) (PTM-1406RM, PTM Bio, Hangzhou, China) = 1:200. Secondary antibodies and dilutions: Anti-rabbit immunoglobulin G fragment antigen-binding 2 Alexa Fluor 488 molecular probes (Cell Signaling Technology, MA, United States) = 1:200. DyLight 594 conjugated to AffiniPure goat anti-rabbit immunoglobulin G (Jackson IR, PA, United States) = 1:100.

Quantitative real-time fluorescence polymerase chain reaction

TRIzol (Thermo Fisher, MA, United States) reagent was used to extract total RNA from normal colonic epithelial cell and CRC cell lines (SW620, HCT116, WS480), CRC tissues, and paracancerous tissues. Complementary DNA was extracted from 500 ng of RNA using HiScript II SuperMix (Vazyme, Nanjing, China). Polymerase chain reaction amplification conditions were 94 °C for 10 minutes, 94 °C for 10 seconds, and 60 °C for 45 seconds. Glyceraldehyde-3-phosphate dehydrogenase served as an internal reference, and the primer pair sequences of all genes used in this study are shown in Table 1.

Table 1 Primer sequences for all genes.

Gene symbol	Forward sequence	Reverse sequence
MT-ND2	CTTCTGAGTCCCAGAGGTTACC	GAGAGTGAGGAGAAGGCTTACG
ACADVL	TAGGAGAGGCAGGCAAACAGCT	CACAGTGGCAAACTGCTCCAGA
SLC2A3	TGCCTTTGGCACTCTCAACCAG	GCCATAGCTCTTCAGACCCAAG
BATF	GATGTGAGAAGAGTTCAGAGGAG	GTTTCTCCAGGTCTTCGCTCTC
RBM3	GACCACTTCAGCAGTTTCGGAC	TGGCTCTCATGGCAACTGAAGC
ICOS	CCCATAGGATGTGCAGCCTTTG	GGCTGTGTTCACTGCTCTCATG
LEPROTL1	GGACTCCCTATTGTATTTGCCAG	GTCGTCCTTGCTTCCAAAGACC
GAPDH	GTCTCCTCTGACTTCAACAGCG	ACCACCCTGTTGCTGTAGCCAA

MT-ND2: Mitochondrial-dihydronicotinamide adenine dinucleotide dehydrogenase subunit 2; ACADVL: Acyl-CoA dehydrogenase very-long-chain; SLC2A3: Solute carrier family 2 member 3; BATF: Basic leucine zipper transcription factor ATF-like; RBM3: RNA binding motif protein 3; ICOS: Inducible T cell costimulatory; LEPROTL1: Leptin receptor overlapping transcript-like 1; GAPDH: Glyceraldehyde-3-phosphate dehydrogenase.

Statistical analysis

R software (version 4.3.1) served as the primary platform for all analytical procedures. P < 0.05 was considered to indicate statistical significance. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001, ^dP < 0.0001, and ns represents non-significant difference.

RESULTS

Cell mapping in CRC

The sequencing data revealed a positive correlation between RNA feature counts (nFeature_RNA) were positively associated with total RNA abundance (nCount_RNA) per cell (r = 0.69), and only a minimal or negligible association was found between nFeature_RNA/nCount_RNA and mitochondrial gene content (Figure 1A). After quality control, sequencing data for 23476 genes in 49269 cells were used for subsequent analysis (Figure 1B). A selection of 2000 genes with prominent transcriptional heterogeneity. were identified after ANOVA with feature selection for the next step of PCA (Figure 1C). We selected the top 28 principal components for cluster analysis. Uniform manifold approximation and projection downscaling analysis revealed 16 cell clusters (Figure 1D). Based on annotations according to the expression levels of marker genes, 16 cell clusters were classified into 11 cell subgroups: Plasmablasts, B cells, monocytes/macrophages, mast cells, endothelial cells, fibroblasts, natural killer cells, CD8+ T cells, epithelial cells, CD4+ T cells, and regulatory T (Treg) cells (Figure 1E and F). Finally, we compared the proportions of 11 cell subpopulations in tumor samples and normal samples. The proportion of monocytes/macrophages was significantly greater in the tumor samples than in the control samples, indicating that monocytes/macrophages play an important role in tumors (Figure 1G).

Figure 1 Cell mapping in colorectal cancer. A: Plot of RNA feature counts (nFeature_RNA) vs total RNA abundance (nCount_RNA) correlation after quality control (QC); plot of nCount_RNA vs mitochondrial percentage after QC; plot of nFeature_RNA vs mitochondrial percentage after QC; B: Violin plots of nFeature_RNA, nCount_RNA, mitochondrial gene percentage, and erythrocyte gene percentage for each sample after QC; C: Volcano plot of highly variable genes. The graph shows 23476 genes in all cells; red dots indicate highly variable genes, and the top 10 most variable genes are labeled in the graph; D: Uniform manifold approximation and projection distribution map of 16 cell clusters identified after cluster analysis and dimensionality reduction analysis; E: Uniform manifold approximation and projection distribution map of 11 cell subclusters identified after annotation; F: Expression bubble plots of marker genes in 11 cell subgroups; G: Ratio of the 11 cell subpopulations in tumor tissues and normal tissues. nFeature_RNA: RNA feature counts; nCount_RNA: Total RNA abundance; UMAP: Uniform manifold approximation and projection; NK: Natural killer.

Monocytes/macrophages exhibit high lactate metabolic activity

Lactate metabolism pathway activity of the 11 cell subpopulations was assessed via the AUCell method, and the results showed that monocytes/macrophages exhibited the highest lactate metabolism activity (Figure 2A and B). Additional confirmation was achieved through AUCell-based evaluation of glycolytic and gluconeogenic activity within defined cellular compartments (Figure 2C), with monocytes/macrophages exhibiting high lactate metabolic activity. Cell communication analysis demonstrated that monocytes/macrophages interact closely with B cells, CD4+ T cells, monocytes/macrophages, and Tregs (Figure 2D). Specifically, monocytes/macrophages interact with B cells and CD8+ T cells through the migration inhibitory factor (CD74+ C-X-C motif chemokine receptor-4) interaction pair and with CD4+ T cells and Tregs through the migration inhibitory factor (CD74+CD44) interaction pair (Figure 2E).

Figure 2 Monocytes/macrophages exhibit high lactate metabolic activity. A: Uniform manifold approximation and projection heatmap of lactate metabolic activity in 11 cell subpopulations; B: Boxplot of lactate metabolic activity in 11 cell subpopulations; C: Uniform manifold approximation and projection heatmap of glycolysis/gluconeogenesis pathway activity in 11 cellular subpopulations; D: Cellular communication network of monocytes/macrophages with the remaining 10 cell subpopulations; E: Ligand-receptor dot plot of monocytes/macrophages with B cells, CD4+ T cells, monocytes/macrophages, and regulatory T cells. UMAP: Uniform manifold approximation and projection; AUC: Area under the curve; NK: Natural killer; TNFSF: Tumor necrosis factor superfamily member; TNFRSF: Tumor necrosis factor receptor superfamily; MIF: Migration inhibitory factor; CXCR4: C-X-C motif chemokine receptor-4; ANXA: Annexin A.

A prognostic signature was established using CRC-associated genes involved in lactate metabolism

Monocyte/macrophage populations were stratified into high- and low-activity subgroups based on the median AUC of lactate metabolic scores. Subsequently, 251 genes exhibiting relevance to lactate metabolism were identified. Through univariate Cox proportional hazards analysis followed by LASSO regression, 14 candidate genes with dual relevance to lactate metabolism and patient prognosis were screened (Figure 3A and B). 7 independent prognostic factors among these genes were identified, which were used to construct a prognostic model for CRC (Figure 3C). The risk score of the prognostic model for CRC was calculated as follows: Risk score = 1.180 × MT-ND2 - 0.214 × RNA binding motif protein 3 (RBM3) + 0.364 × ACADVL + 0.370 × SLC2A3 + 0.130 × BATF -0.205 × ICOS - 0.346 × LEPROTL1. For the TCGA-COADREAD cohort, the cohort was divided into two prognostic groups using the median value as the cutoff. MT-ND2, ACADVL, SLC2A3, and BATF exhibited significantly higher transcriptional activity in the high-risk subgroup, and RBM3, ICOS, and LEPROTL1 displayed downregulated expression patterns in high-risk patients (Figure 3D). The distribution of the risk score across patients was assessed, and an elevated risk score was associated with increased mortality and shortened overall survival (Figure 3E and F). The PCA results showed that the risk score model could better distinguish the high-risk score samples from the low-risk score samples (Figure 3G). In the TCGA-COADREAD cohort, a notable survival advantage was observed in patients classified as low risk (Figure 3H). Time-dependent AUCs for 1-, 3-, and 5-year survival prediction were calculated as 0.67, 0.67, and 0.64, respectively (Figure 3I). The predictive value of the risk score was further validated in the GSE17536 validation set; the model had AUC values of 0.64, 0.69, and 0.76 for predicting patient survival at 1, 3, and 5 years, respectively, in this dataset (Figure 3J). These results indicate that the risk score model, which consists of seven genes related to lactate metabolism, has good predictive performance.

Figure 3 Prognostic modeling of lactate metabolism characteristics in colorectal cancer. A: Screening for the optimal penalty parameter (λ) in LASSO Cox analysis; B: Coefficients of genes in LASSO Cox analysis; C: Forest plot of seven independent prognostic factors; D: Heatmap of the expression levels of 7 lactate metabolism genes in the high- and low-risk score groups; E: Risk score distribution of samples in The Cancer Genome Atlas database-colon and rectal adenocarcinomas (TCGA-COADREAD) cohort; F: Scatterplot of the survival status of COADREAD patients in the TCGA-COADREAD cohort according to the risk score distribution; G: Principal component analysis distribution plot of patients with high and low risk scores based on the expression profiles of 7 genes related to lactate metabolism; H: Kaplan-Meier survival curves of patients in the TCGA-COADREAD cohort with high and low risk scores; I: Receiver operating characteristic curves of patients in the TCGA-COADREAD cohort with the risk score for predicting 1-, 3-, and 5-year survival; J: Receiver operating characteristic curves for the risk score for predicting patient survival at 1, 3, and 5 years in the GSE103479 cohort. ^aP < 0.05, ^cP < 0.001. MT-ND2: Mitochondrial-dihydronicotinamide adenine dinucleotide dehydrogenase subunit 2; RBM3: RNA binding motif protein 3; ACADVL: Acyl-CoA dehydrogenase very-long-chain; SLC2A3: Solute carrier family 2 member 3; BATF: Basic leucine zipper transcription factor ATF-like; ICOS: Inducible T cell costimulator; LEPROTL1: Leptin receptor overlapping transcript-like 1; PCA: Principal component analysis; ROC: Receiver operating characteristic; AUC: Area under the curve; TCGA: The Cancer Genome Atlas database.

Immune infiltration and immunotherapy sensitivity in the risk score subgroups

According to the single-sample gene set enrichment analysis results, the infiltration of 26 immune cell types was greater in the high-risk score group (Figure 4A). ESTIMATE analysis was also performed in this study, which revealed that high-risk subgroup demonstrated elevated stromal and ESTIMATE scores relative to the low-risk subgroup. Higher tumor purity was observed in patients classified as low-risk (Figure 4A-E). Cell-type identification by estimating relative subsets of RNA transcript analysis revealed that the immune landscape of low-risk patients was characterized by elevated levels of plasma cells, activated CD4 memory T cells, and antigen-presenting dendritic cells, while enhanced accumulation of immunosuppressive M0 macrophages and mast cell subsets was noted in the high-risk cohort (Figure 4F). In addition, associations were assessed between risk scores, immune checkpoint gene expression, and IPS values. An inverse relationship was detected between the risk score and the expression of cytotoxic T lymphocyte antigen 4 and programmed death-ligand 1 (CD274) (Figure 4G and H). Finally, stratified IPS analysis was conducted under different cytotoxic T lymphocyte antigen 4/programmed death-1 expression scenarios in both risk groups. Across all immune checkpoint subgroups, IPSs were higher in low-risk individuals (Figure 4I-L) (P < 0.05). High IPS may indicate strong immune cell infiltration, suggesting that a heightened immunotherapeutic responsiveness in patients from the low-risk group.

Figure 4 Immune infiltration and immunotherapy sensitivity in the risk score subgroups. A: Heatmap of immune cell infiltration levels according to single-sample gene set enrichment analysis results and estimation of stromal and immune cells in malignant tumor tissues using expression data results; B: Violin plot demonstrating the immune score in the high- and low-risk score groups; C: Violin plot demonstrating the stromal score in the high- and low-risk score groups; D: Violin plot showing the estimation of stromal and immune cells in malignant tumor tissues using expression data score in the high- and low-risk score groups; E: Violin plot demonstrating tumor purity in the high- and low-risk score groups; F: Cell type identification by estimating relative subsets of RNA transcripts results demonstrating the infiltration levels of 10 immune cells in the high- and low-risk score groups; G: Scatter plot of the correlation between the risk score and cytotoxic T lymphocyte associated protein 4 expression level; H: Scatterplot of the correlation between the risk score and CD274 expression level; I: Immunopheno score (IPS) of cytotoxic T lymphocyte antigen 4 (CTLA-4)-negative/programmed death-1(PD-1)-negative colorectal cancer (CRC) patients in the high- and low-risk score subgroups; J: IPS of CTLA-4-negative/PD-1-positive CRC patients in the high- and low-risk score subgroups; K: IPS of CTLA-4-positive/PD-1-negative CRC patients in the high- and low-risk score subgroups; L: IPS of CTLA-4-positive and PD-1-positive CRC patients in the high- and low-risk score subgroups. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001, ^dP < 0.0001, denote significant difference from the low-risk score subgroups. ESTIMATE: Estimation of stromal and immune cells in malignant tumor tissues using expression data; CTLA: Cytotoxic T lymphocyte antigen.

Gene mutation landscape of patients in the risk score subgroups

By comparing mutation patterns, tumor mutation burden (TMB), and interactions between genes, we aimed to reveal differences in molecular features between different risk score groups. The classification of missense mutations, frameshift mutations, and nonsense mutations and single nucleotide variant category analysis were used to further refine mutations into T>G, T>A, T>C, C>T, C>G, and C>A variant types. The median number of variants in the low-risk score group was 86 (Figure 5A), and the median number of variants in the high-risk score group was 92 (Figure 5B), suggesting increased TMB in the high-risk score group. Mutation detection was achieved in 92.23% of low-risk group samples (Figure 5C) and mutation-positive rates reached 96.89% in high-risk group samples, with higher TP53 mutation frequency was noted in the high-risk subgroup (Figure 5D). We also found a more dispersed relationship between gene comutation and mutually exclusivity in the low-risk group, while the high-risk group demonstrated a more pronounced pattern of comutation (Figure 5E and F). Overall, a strong correlation was observed between patient risk classification and mutation-based molecular characteristics. An increased TMB was identified in the high-risk cohort and greater prevalence of specific gene mutations.

Figure 5 Landscape of mutations in patients in the risk score subgroups. A: Distribution of single nucleotide variant types in tumor samples in the low-risk score group; B: Distribution of single nucleotide variant types in tumor samples in the high-risk score group; C: Most common mutated genes in the low-risk score group and their mutation frequencies in the samples; D: The most common mutated genes in the high-risk score group and their mutation frequencies in the samples; E: Analysis of the relationships between cooccurring genes and mutually exclusive genes in the low-risk score group; F: Analysis of gene cooccurrences and mutually exclusive relationships in the high-risk score group. ^aP < 0.05. SNV: Single nucleotide variant; SNP: Single nucleotide polymorphisms; ONP: Oligonucleotide polymorphism; INS: Insertion; DEL: Deletion; TMB: Tumor mutation burden.

Small molecule drug therapy sensitivity

Tumors in the high-risk score group showed greater sensitivity to the three drugs (AZD6482, tozasertib, and SB216763) than did those in the low-risk score group (Figure 6). The significant difference in the half-maximal inhibitory concentration values suggested that a tumor’s risk score may be related to its responsiveness to specific drugs. This finding provides support for personalized treatment strategies based on the risk score, suggesting that patients with a high risk score may be better suited for treatment with these drugs (P < 0.05).

Figure 6 Sensitivity to small molecule drugs. A: Half-maximal inhibitory concentration (IC₅₀) of AZD6482 in the high- and low-risk score groups; B: IC₅₀ of tozasertib in the high- and low-risk score groups; C: IC₅₀ of SB216763 in the high- and low-risk score groups. ^aP < 0.05, ^dP < 0.0001, denote significant difference from the low-risk score groups. IC₅₀: Half-maximal inhibitory concentration.

Nomogram model for predicting survival in CRC patients

We performed univariate Cox and multivariate Cox analyses. The results showed that the risk score was an independent prognostic factor for CRC (Figure 7A and B). Based on the patients’ clinical information, we plotted a nomogram model for predicting the 1-, 3-, and 5-year survival rates of CRC patients (Figure 7C). The accuracy of the prediction model was validated by comparing the actual survival rate with the survival rate predicted by the model. For 1-, 3-, and 5-year survival predictions, the model showed good predictive ability, especially for predicting 1-year survival (Figure 7D-F). In addition, we compared the net clinical benefit of clinical parameter-based models and the nomogram model with taking no action or taking all possible actions via decision curves. The results showed that using the nomogram prediction model provided a greater net gain in decision-making at most threshold probabilities (Figure 7G-I).

Figure 7 Nomogram model for predicting survival in colorectal cancer patients. A: Forest plot of univariate Cox analysis combining clinical information and the risk score in The Cancer Genome Atlas database-colon and rectal adenocarcinomas cohort; B: Forest plot of multivariate Cox analysis combining clinical information and the risk score in The Cancer Genome Atlas database-colon and rectal adenocarcinomas cohort; C: Nomogram model constructed by integrating clinical information and the risk score; D: Calibration curve for 1-year survival; E: Calibration curve for 3-year survival; F: Calibration curve for 5-year survival; G: Decision curve for 1-year survival; H: Decision curve for 3-year survival; I: Decision curve for 5-year survival. T: Tumor; N: Node; M: Metastasis.

Expression levels of biomarkers in CRC

H3K18 La expression was assessed in paired CRC and adjacent non-tumor tissues, revealing notably elevated levels in tumor samples (Figure 8A). Red fluorescence (representing the level of H3K18 La lactonization) showed intensified expression in malignant tissues compared with adjacent normal tissues. This intense fluorescence demonstrating increased lactonization of H3K18 La in cancerous vs non-tumorous samples (Figure 8B). The gene expression levels and protein expression levels of MT-ND2, ACADVL, SLC2A3, BATF, RBM3, ICOS, and LEPROTL1 were detected in FHC and CRC cell lines (SW620, HCT116, WS480), CRC tissues, and paracancerous tissues. Immunoblotting results confirmed the upregulation of MT-ND2, ACADVL, SLC2A3, and BATF were upregulated in SW620, HCT116, WS480, and CRC tissues, and these proteins (RBM3, ICOS, LEPROTL1) showed downregulated expression patterns (Figure 8C). The quantitative real-time fluorescence polymerase chain reaction analysis revealed altered expression levels of MT-ND2, ACADVL, SLC2A3, and BATF were significantly elevated in SW620, HCT116, WS480, and CRC tissues, and conversely, the expression levels of RBM3, ICOS, and LEPROTL1 were decreased (P < 0.05) (Figure 8D and E).

Figure 8 Expression levels of biomarkers in colorectal cancer. A: Immunohistochemistry results demonstrating the expression level of histone H3 lysine 18 lactylation in colorectal cancer (CRC) tissues and paracancerous tissues; B: Immunofluorescence results of histone H3 lysine 18 lactylation levels in tumor tissues and paracancerous tissues; C: Western blot results demonstrating the protein expression levels of mitochondrial-dihydronicotinamide adenine dinucleotide dehydrogenase subunit 2 (MT-ND2), acyl-CoA dehydrogenase very-long-chain (ACADVL), solute carrier family 2 member 3 (SLC2A3), basic leucine zipper transcription factor ATF-like (BATF), RNA binding motif protein 3 (RBM3), inducible T cell costimulatory (ICOS), leptin receptor overlapping transcript-like 1 (LEPROTL1) in CRC tissues/paracancerous tissues, and normal colon cell/CRC cell lines (SW620, HCT116, WS480); D: The quantitative real-time fluorescence polymerase chain reaction results demonstrating the expression levels of MT-ND2, ACADVL, SLC2A3, BATF, RBM3, ICOS, and LEPROTL1 in CRC tissues and paracancerous tissues; E: The quantitative real-time fluorescence polymerase chain reaction results demonstrate the expression levels of MT-ND2, ACADVL, SLC2A3, BATF, RBM3, ICOS, and LEPROTL1 in normal colon cell and CRC cell lines (SW620, HCT116, WS480). ^cP < 0.001, compared with normal group. H3K18 La: Histone H3 lysine 18 lactylation; MT-ND2: Mitochondrial-dihydronicotinamide adenine dinucleotide dehydrogenase subunit 2; ACADVL: Acyl-CoA dehydrogenase very-long-chain; SLC2A3: Solute carrier family 2 member 3; BATF: Basic leucine zipper transcription factor ATF-like; RBM3: RNA binding motif protein 3; ICOS: Inducible T cell costimulatory; LEPROTL1: Leptin receptor overlapping transcript-like 1; FHC: Normal colon cells.

DISCUSSION

This study provides a comprehensive analysis of lactate metabolism in CRC and its association with immune infiltration, gene mutations, and patient prognosis. Our findings reveal that lactate metabolism plays a critical role in regulating monocyte/macrophage function, including their functions in inflammation, the immune response, and the TME. Lactate significantly suppresses the secretion of tumor necrosis factor and interleukin-1beta by human macrophages in response to Mycobacterium tuberculosis, and the glycolytic product lactate exerts a negative feedback effect on macrophages, leading to attenuated glycolytic shifts and reduced proinflammatory cytokine production upon subsequent stimulation[31]. Lactate and related H⁺ ions are regarded as negative regulators of immunosuppression, and they potentiate the point-of-function-specific effects of M2-type macrophages[31]. M1-type macrophages initiate and maintain inflammatory responses, secrete proinflammatory cytokines, activate endothelial cells, and induce other immune cells to enter inflamed tissues. On the other hand, M2-type macrophages promote the resolution of inflammation, phagocytose apoptotic cells, promote collagen deposition, orchestrate tissue integrity, and release anti-inflammatory mediators. These phenotypic and functional changes are accompanied by macrophage lactate metabolism or other types of metabolic activity. Macrophage function and polarization are closely linked to metabolic changes, with M1-type macrophages relying primarily on aerobic glycolysis and M2-type macrophages relying primarily on oxidative metabolism[32]. Lactate metabolism in cancer is a complex metabolic reprogramming process whereby cancer cells produce lactate via aerobic glycolysis, even in the presence of an adequate oxygen supply. Moreover, lactate is utilized not only as an energy source by oxidized cancer cells in the TME but also as a signaling molecule and a potential target for cancer therapy[10,33]. Taken together, our results indicate that CRC tissues have a greater monocyte/macrophage ratio and exhibit high lactate metabolic activity than normal tissues. We speculate that enhanced lactate metabolism in macrophages may increase the lactate concentration in the TME, thereby enhancing the immunosuppressive microenvironment and supporting tumor cell survival and proliferation. In a high-lactate environment, macrophages may be more inclined to polarize toward the M2 phenotype, which is associated with the promotion of tissue repair, anti-inflammatory effects, and tumor growth.

The prognostic model we developed, based on seven lactate metabolism-related genes (MT-ND2, ACADVL, SLC2A3, BATF, RBM3, ICOS, and LEPROTL1), demonstrated robust predictive power for patient outcomes. MT-ND2 is a gene that encodes a subunit of the dihydronicotinamide adenine dinucleotide dehydrogenase complex, which is involved in the oxidative phosphorylation of cells. It may play an important role in cancer metabolism, especially in environments associated with mitochondrial function and metabolic reprogramming[34]. ACADVL is an enzyme involved in the beta-oxidation of long-chain fatty acids, which regulates fatty acid metabolism[35]. SLC2A3 is a glucose-transport protein involved in the cellular uptake of glucose. In the context of cancer, changes in SLC2A3 expression alter glycolytic processes in tumor cells[36]. Another study revealed that SLC2A3 promotes macrophage infiltration in gastric cancer through glycolytic reprogramming[36]. BATF is a transcription factor participates in immune cell lineage specification and immunological function and positively regulates immune responses in the TME[37]. RBM3, an RNA-binding protein, is involved in cancer progression by regulating the cell cycle and cellular stress responses[38]. ICOS is an immune costimulatory molecule that plays an important role in T-cell activation and function. In CRC, the interaction between ICOS and its ligands may affect the antitumor immune response[39]. The function of LEPROTL1 in cancer is unknown, and further studies are needed to explore its function. Overall, we provided a useful tool for assessing the prognosis of CRC. This prognostic model reflects the central role of lactate metabolism in CRC development and provides a potential biomarker for future therapeutic strategies. The immune infiltration analysis further underscores the role of lactate metabolism in shaping the TME. Significant variations in immune cell infiltration patterns were identified between the high- and low-risk groups, underscoring the value of integrating risk stratification with immune profiling to inform individualized immunotherapy strategies. Mutational landscape analysis indicated an increased prevalence of specific gene mutations and elevated TMB in the high-risk group, potentially accounting for the observed immune microenvironment disparities. Furthermore, drug sensitivity assessment suggested that tumors classified as high-risk were more responsive to certain small-molecule agents, reinforcing the rationale for tailoring treatment strategies based on distinct molecular features.

The nomogram model, which integrates clinical and molecular data, further enhances the utility of our prognostic model by serving as a practical instrument for short- and long-term survival prediction in CRC cohorts. The model demonstrated good accuracy in predicting patient outcomes, particularly for short-term survival. This reinforces the importance of combining molecular profiling with clinical characteristics to improve prognostic accuracy and guide therapeutic decision-making in CRC. Notably, we constructed a model with an AUC around 0.7. We speculate that it is due to the following reasons; CRC is a highly heterogeneous cancer, and the possible reasons include gene mutations, microsatellite instability, epigenetic alterations, and heterogeneity of the TME[40,41]. Secondly, patients may have differences in the expression of prognostic markers in the patient’s body due to different treatment regimens. Due to the complexity of biological data, this means that even with good methods and models, the AUC may still be at a low level, as differences between individuals may not be perfectly captured by simple statistical models. This study also analyzed the expression level of H3K18 La in CRC. H3K18 La is a direct epigenetic marker of lactate metabolism. Lactate is not only a product of energy metabolism, but also modifies histones through lactylation and directly regulates gene expression. Elevated H3K18 La indicates active lactate metabolism, which may promote facilitating malignant evolution via epigenetic modulation[42,43]. H3K18 La levels were significantly higher in CRC tissues than in paracancerous tissues, confirming the biological effects of lactate accumulation and its metabolic reprogramming in the TME. The high expression of H3K18 La may reflect tumor aggressiveness, which is consistent with the expression pattern of lactate metabolism genes in the prognostic model, further supporting the central role of lactate metabolism in CRC. The analysis of H3K18 La revealed a significant influence of lactate metabolism on tumor progression at an epigenetic level, and a significant role of lactate metabolism in CRC. H3K18 La level revealed the regulatory mechanism of lactate metabolism in CRC, providing key evidence for understanding the tumor metabolism-immunity-epigenetic interaction network.

This study still has limitations. In this study, we mainly screened prognostic factors by univariate Cox, LASSO Cox, and multivariate Cox, although it showed good prognostic prediction performance. However, we consider the subsequent introduction of multiple machine learning algorithms to screen the best prognostic factors and prognostic models for CRC by comparing the Matthews correlation coefficients or through receiver operating characteristic curves. On the other hand, most of the 7 prognostic genes identified in this study were found to be associated with cancer prognosis for the first time, and the potential molecular mechanisms deserve to be subsequently investigated in depth, which is the focus of our subsequent work plan.

CONCLUSION

This work emphasizes the importance of lactate metabolic activity in driving CRC progression, modulating immune responses, and influencing patient prognosis. By integrating lactate metabolism-related genes into a prognostic model, we provide a valuable tool for predicting patient outcomes and guiding personalized therapy. Our findings suggest that targeting lactate metabolism, in combination with immunotherapy or small molecule inhibitors, holds promise for improving treatment outcomes in CRC. Additional investigations are needed to assess the therapeutic value of targeting lactate metabolism and its potential relevance to the TME in various malignancies.