Multi-omics analysis reveals immunosuppression in oesophageal squamous cell carcinoma induced by creatine accumulation and HK3 deficiency

Abstract

Background: Deep insights into the metabolic remodelling effects on the immune microenvironment of oesophageal squamous cell carcinoma (ESCC) are crucial for advancing precision immunotherapies and targeted therapies. This study aimed to provide novel insights into the molecular landscape of ESCC and identify clinically actionable targets associated with immunosuppression driven by metabolic changes.

Methods: We performed metabolomic and proteomic analyses combined with previous genomic and transcriptomic data, identified multi-omics-linked molecular features, and constructed metabolic-immune interaction-based ESCC classifiers in a discovery cohort and an independent validation cohort. We further verified the molecular characteristics and related mechanisms of ESCC subtypes.

Results: Our integrated multi-omics analysis revealed dysregulated proteins and metabolic imbalances characterizing ESCC, with significant alterations in metabolites and proteins linked to genetic traits. Importantly, ESCC patients were stratified into three subtypes (S1, S2, and S3) on the basis of integrated metabolomic and proteomic data. A robust subtype prediction model was developed and validated across two independent cohorts. Notably, patients classified under the poorest prognosis subtype (S3 subtype) exhibited a significant immunosuppressive microenvironment. We identified key metabolism-related biomarkers for the S3 subtype, specifically creatine and hexokinase 3 (HK3). Creatine accumulation and HK3 protein deficiency synergistically reprogrammed macrophage metabolism, driving M2-like TAM polarization. This metabolic shift fostered an immunosuppressive microenvironment that accelerated tumour progression. These results highlight the potential of targeting creatine metabolism to improve the efficacy of immunotherapy and targeted therapy for ESCC.

Conclusions: Our analysis reveals molecular variation in multi-omics linkages and identifies targets that reverse the immunosuppressive microenvironment through metabolic remodelling improving immunotherapy and targeted therapy for ESCC.

Keywords: Creatine; Hexokinase 3; Immunosuppressive microenvironment; Metabolomics; Multi-omics; Oesophageal squamous cell carcinoma; Proteomics; Subtypes.

Fig. 1

Overview of the experimental strategy, and proteomic landscape of oesophageal squamous cell carcinoma (ESCC) compared with NATs. A Overview of the experimental strategy, image created with BioRender.com, with permission. B Principal component analysis (PCA) of the proteome. Red, tumours; blue, NATs. C Number of proteins identified in tumours (red dots) and NATs (blue dots). The 95% confidence intervals are shown by the hidden region. D Heatmap showing the differentially expressed proteins in tumours and NATs. Significantly up-regulated and down-regulated pathways are shown on the right. E Heatmap showing the proteins associated with patient survival in these significantly changed pathways, and log2-based hazard ratio are shown on the right. HR: hazard ratio for OS. F Representative proteins of four biological pathways and their relationship with prognosis (FDR-corrected log-rank P values). The cut-off values of each protein are as follows: SNRPB (18.1421), SF3A1 (16.6450), EIF2S1 (18.6551), DNAJC10 (13.1860), DLD (17.6574), BCKDHA (15.9849), COX6B1 (17.2751), UQCRC1 (18.0475). See also Fig. S1, S2 and Table S1, S3, S4

Fig. 2

The metabolomic landscape of oesophageal squamous cell carcinoma. A Proportions of annotated metabolites identified in our study. B Volcano plots of the annotated metabolites. C A pathway-based analysis of metabolomic changes between tumour and NATs. The DA score captures the average, gross changes for all metabolites in a pathway. A score of 1 indicates that all measured metabolites in the pathway increase in the tumour compared to normal tissues, and a score of −1 indicates that all measured metabolites in a pathway decrease. Pathways with no fewer than three measured metabolites were used for the DA score calculation. D Pathway abundance (PA) scores between tumour and NATs. The PA score was calculated as the mean log2 fold change in the abundances of the measured metabolites in this pathway. E A metabolic map profiling the synthesis and degradation of several amino acids, based on metabolomics data. See also Fig. S2 and Table S5

Fig. 3

Effects of mutations and copy number alterations (CNAs) on mRNA and protein abundance. A Functional effect of mutations on mRNA and proteins. The y-axis shows significant mutant genes in ESCC, and the x-axis is cis- or trans- effected genes and related pathways. B Correlation of CNA to mRNAs and protein abundance. Positive and negative correlations are indicated in red and blue, respectively. Genes are ordered by chromosomal location on the x and y axes. The diagonal lines indicate the cis-effects of CNA on mRNAs or proteins. C Overlap of cis-effects observed at mRNA and proteins (FDR < 0.05). D KEGG pathway enrichment analysis of overlapped RNA and proteins in C. E CNA frequency diagram. Red for amplification, blue for deletion. F Volcano plot showing log2-based hazard ratio for each significant CNA peak regions. G Kaplan–Meier curves for overall survival analysis of patients with 1p34 gain or 13q22 gain (P value from log-rank test). H Heatmap showing the normalized expression of cis- and trans-effecting proteins and mRNAs significantly associated with copy number amplification in the 1p34 region (left) and log2-based hazard ratio (right). HR: hazard ratio for OS. The P-values were adjusted using Benjamini–Hochberg false discovery rate (FDR) correction. I Left: Heatmap showing the score of the proteasome pathway and the normalized expression of associated proteins. Centre: Log2-based hazard ratios and FDR of these proteins. Right: Spearman’s correlation coefficients of these proteins with PPCS expression (blue) and proteasome pathway score (red). HR: hazard ratio for OS. See also Fig. S3 and Table S7

Fig. 4

ESCC molecular subtyping based on integrated proteomics and metabolomics analysis. A Heatmap with clinical characteristics of ESCC samples into three SNF-derived subtypes: S1 (n = 47), S2 (n = 32), and S3 (n = 47) based on integrated analysis of proteomics and metabolomics. Pathways are significantly enriched in each subtype. B Kaplan–Meier curves for overall survival and progression-free survival of different subtypes (P value from log-rank test). C Number of patients with lymphatic metastasis between three subtypes (P value from chi-square test). D Heatmap showing the median number of immune cell infiltration in NATs and three subtypes, and log2-based hazard ratio for each immune cell infiltration score. E Multiple immunofluorescence results and intensity statistical analysis. Box plots show the percentage of positive cells in subtypes S1 (n = 7) and S3 (n = 7). Data presented as mean ± s.d.; P values by two-tailed Student’s t test

Fig. 5

Molecular subtype validation in another independent cohort. A The contribution of six signatures to the subtype diagnostic model. B Receiver operating characteristic (ROC) curves with reported areas under the curve (AUCs) demonstrated the efficacy of the subtype diagnostic model in identifying subtypes of ESCC. C A abundance of the six signatures in three subtypes. Data presented as mean ± s.d.; P values by Wilcoxon rank-sum test. D IHC and intensity statistics of four characteristic proteins in the predicted S3 subtype (n = 12) and S1/2 subtype (n = 40) in the independent ESCC cohort 2 (P value from Wilcoxon rank-sum test). Data presented as mean ± s.d.; P values by Wilcoxon rank-sum test. E, F The intensity of creatine (E) and 2’DG (F) in the predicted S3 subtype (n = 12) and S1/2 (n = 40) subtype in cohort 2. Data presented as mean ± s.d.; P values by Wilcoxon rank-sum test. G Kaplan–Meier curves for overall survival analysis of predicted S3 subtype and S1/2 subtype (P value from the log-rank test). H, I Spearman correlation analysis of the protein expression of HK3 protein and the intensity of creatine in our study (H) or in cohort 2 (I). J ROC curves with reported AUCs demonstrated the efficacy of the model containing only creatine and HK3 protein. K, L Kaplan–Meier curves for overall survival analysis of the expression of HK3 protein (K) and the intensity of creatine in our study (L) (P value from log-rank test). See also Table S10

Fig. 6

Analysis of immune infiltration associated with the abundance of creatine and HK3. A Changes in the abundance of metabolites, proteins and mRNAs of the creatine synthesis. The log2-fold changes of metabolites, proteins and mRNAs in S3 and S1 subtypes are colour-coded (*P < 0.05; ** P < 0.01; *** P < 0.001, see also Fig. S6 A, Supporting Information). B Spearman’s correlation coefficients of the abundance of creatine synthesis-related molecules and each immune cell infiltration score (only shown FDR < 0.05). C Spearman’s correlation coefficients of immune score and the expression of HK3 protein and mRNA. D Spearman’s correlation coefficients of each immune cell infiltration score and the expression of HK3 protein and mRNA (only shown P < 0.05). E Functional enrichment analysis of the proteins with a significant correlation (|Spearman’s correlation coefficient|> 0.3, FDR < 0.05) with the expression of HK3 in proteomics. Red: positive correlation. Blue: negative correlation. F Left: Heatmap showing the expression of immune-related proteins significantly positive correlation with the expression of HK3 in proteomics. Right: Log10-based hazard ratio and FDR-corrected P value of these immune-related proteins. HR: hazard ratio for OS. See also Fig. S6 and S7

Fig. 7

Hk3 deficiency and creatine treatment reprogrammed macrophage polarization in vivo and in vitro. A Experimental design. Hk3-ko or WT mice were fed creatine-containing water or not. B, C Hk3-ko and creatine-containing diets promoted tumour growth in subcutaneously injected model mice (B). Boxplot showing the difference of tumour weight among different groups (C) (Data presented as mean ± s.d.; P values were calculated via two-tailed Student’s t test; n = 6). Dots refer to independent samples. D–F Representative flow cytometry analysis of tumour fractions from different groups. Barplot showing the discrepancy in pan TAMs (D, Cd45⁺F4/80⁺Cd11b⁺ cells), M2-like TAMs (E, Cd45⁺F4/80⁺Cd11b⁺Cd206⁺ cells) and M1-like TAMs (F, Cd45⁺F4/80⁺Cd11b⁺Cd86⁺Mhc-II⁺cells) among different groups. Data presented as mean ± s.d.; P values were calculated via two-tailed Student’s t test; n = 3. G, H Effect of Hk3 knockdown (Hk3 sh) and 1 mM creatine pretreatment (+ Cr.) on the chemotactic ability of RAW264.7 cells, with IL-4 for 24 h (H) or not (G), co-cultured with mEC25 cells. Data presented as mean ± s.d.; P values were calculated via two-tailed Student’s t test; n = 5. I Hk3 knockdown (Hk3 sh) and 1 mM creatine pretreatment (+ Cr.) promoted M2-like TAM polarization in RAW264.7 cells. Data presented as mean ± s.d.; P values were calculated via two-tailed Student’s t test; n = 3. See also Fig. S8 and Table S2