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.

Background

Oesophageal cancer (EC) is the 11th most commonly diagnosed cancer and the seventh leading cause of cancer death worldwide, with oesophageal squamous cell carcinoma (ESCC) accounting for approximately 90% and a 5-year survival rate of 20% [1]. Compared with other common cancer types, EC, especially ESCC, is deeply understudied, and there has been limited progress in therapeutics in recent decades. Although immune checkpoint inhibitors (ICIs) show good efficacy, only 30–40% of patients benefit from ICIs and resistance to ICI therapy is very common [2,3,4]. In-depth understanding of microenvironmental heterogeneity and immunosuppression mechanism has led to the identification of new targets for activating tumour immunity and developed a series of immunomodulatory agents in cancers. The monotherapy of these agents and combination with other anti-tumour therapies have been shown to overcome the tolerance and resistance of immunotherapy. For example, supplementation of arginine or arginase inhibitors treatment improves cancer immunotherapy outcomes by promoting the proliferation of T cells and NK cells [5]. A2aR antagonists have been shown to potentiate anti-tumour immunity in renal cell carcinoma [6], prostate cancer [7], and colorectal cancer [8] by inhibiting the secretion of immunosuppressive factors by dendritic cells (DCs).

In the last decade, hundreds of ESCCs have been analysed at a multi-omics level including genomics, transcriptomics, epigenomics, single-cell RNA sequencing and spatial transcriptomics. These studies identified molecular drivers of carcinogenesis, cell lineage trajectories, and prognostic subtypes [9,10,11,12,13], but lacked integration with metabolomics to link metabolic remodelling to the immune microenvironment. No agents have been identified to reverse the immunosuppressive microenvironment, specifically through metabolic remodelling, improving immunotherapy and targeted therapy for ESCC.

To address the gap in understanding metabolic-immune interactions in ESCC, we integrated metabolomic and proteomic profiling with existing genomic and transcriptomic data from treatment-naive tumours and paired normal tissues. This multi-omics approach aimed to elucidate the metabolic abnormalities linked to multi-omics disruption and systematically map metabolic reprogramming modulating the tumour microenvironment. We further established molecular subtypes and constructed subtype diagnostic and prognostic models with clinical relevance, providing a foundation for targeted therapeutic strategies in ESCC.

Methods

Clinical sample acquisition

In this study, two cohorts of ESCC patients were recruited from Shanxi Cancer Hospital, with Cohort 1 (ESCC cases for metabolomics and proteomics study) as the discovery cohort and Cohort 2 (52 cases for immunohistochemistry and targeted metabolomics analysis) as the validation cohort. In cohort 1, 136 tumour tissues and 136 paired normal adjacent tissues (NATs) were used for metabolomics study, of which 127 tumour tissues and 128 paired NATs underwent proteomic analysis. Informed consent was obtained from all subjects, and the study was approved by the ethical committees of the Shanxi Medical University. NATs were collected with a minimum distance of 5 cm from the tumours. None of the patients received radiotherapy, chemotherapy or other therapy before the operation. Post-resection, the tissue specimens were rapidly frozen using liquid nitrogen within 30 min and preserved in a freezer at − 80 °C. Every patient was categorized based on the pTNM system's eighth edition as per the American Joint Committee on Cancer (AJCC). Haematoxylin and eosin (H&E)-stained sections were reviewed by at least three independent pathologists to confirm that the tumour specimen was consistent with ESCC and that the NATs contained no tumour cells. The detailed clinicopathological information of all patients is summarized in Table S1 (Additional file 2). Among the tumour-NAT pairs in cohort 1, 120 tumour-NAT pairs had previously been included in other multi-omics studies. Among them, 80 pairs with transcriptome and whole-genome sequencing (WGS) data came from our previous study of PMID:36584672 [13, 14], and 40 pairs with transcriptome and whole-exome sequencing (WES) data from our another previous study of PMID: 38803565 [15, 16].

Mice

Hk3 knockout C57BL/6 mice were purchased from Cyagen Biosciences, Inc. (Jiangsu, China). Three- to four-week-old female C57BL/6 mice and 6-week-old female NTG mice were used. These mice were purchased from SPF (Beijing) Biotechnology Co., Ltd. The animal experiments were approved by the Ethics Committee of Shanxi Medical University (SYDL2023040).

Cell lines

The ESCC cell lines (KYSE-150, TE-1 and mEC25) were kindly provided by Professor Qimin Zhan (Peking University, Beijing, China) and Professor Li Fu (Shenzhen University, Shenzhen, China). RAW264.7 cells were purchased from Procell Life Science & Technology Co., Ltd. (Wuhan, China). KYSE-150 and TE-1 were maintained in RPMI 1640 (VivaCell, #C3001 - 0500), 10% foetal bovine serum, and 1 × penicillin/streptomycin (VivaCell, #C3420 - 0100). mEC25 and RAW264.7 were maintained in DMEM (VivaCell, #C3103 - 0500), 10% foetal bovine serum, and 1 × penicillin/streptomycin. All cells were kept in a 5% CO₂, 37 ℃ incubators. All cell lines were routinely verified using short tandem repeat analysis and were tested for mycoplasma before their use in experiments.

Data-independent acquisition (DIA) proteomics analysis

Tissues were homogenized in a centrifuge tube containing SDS L3 and EDTA buffer. Subsequently, the specimens were cooled on ice and combined with DTT to achieve a concentration of 10 mM. The tissue was crushed via a grinder operating at 60 Hz for 2 min. Post-centrifugation, the supernatant was discarded, followed by the addition of DTT to achieve a concentration of 10 mM. Following a 1 h incubation at 56 °C in a water bath, IAM was introduced to achieve a concentration of 55 mM, after which the samples were kept in darkness for 45 min. Subsequently, chilled acetone was introduced into the protein mixture, underwent centrifugation, and the resulting supernatant was disposed of. After the samples were air-dried, lysis buffer was added, followed by centrifugation, and the resulting supernatant was carefully gathered. The quality of the extracted proteins underwent quality assurance via Bradford quantification and SDS‒PAGE. In the process of protein enzymatic hydrolysis, a mixture of 2.5 μg of trypsin enzyme and 100 μg of protein solution was prepared at a 40:1 ratio, followed by a 4 h digestion at 37 °C. The peptides were desalted via a Strata X column and vacuum-dried.

Peptides were extracted, pooled, and diluted in mobile phase A (5% ACN, pH 9.8) for high-pH reverse-phase separation. A Shimadzu LC-20 AD liquid phase system with a 5 μm 4.6 × 250 mm Gemini C18 column was used to separate all the samples. The sample was then eluted at a flow rate of 1 mL/min by a gradient: 5% mobile phase B (95% ACN, pH 9.8) for 10 min, 5% to 35% mobile phase B for 40 min, 35% to 95% mobile phase B for 1 min, and flow phase B for 3 min and equilibration for 10 min in 5% mobile phase B. Throughout the process, elution peaks were meticulously monitored at a wavelength of 214 nm, with components being systematically collected every minute. The amassed components were partitioned into a total of 10 fractions and subsequently subjected to freeze-drying procedures.

Metabolic profiling

The extraction of metabolites was carried out primarily according to methods previously reported [17]. By pooling the same volume of each sample, a quality control sample (QC) was prepared. Separation and detection of metabolites were conducted using a tandem Q Exactive HF high-resolution mass spectrometer (Thermo Fisher Scientific, USA) from Waters 2D UPLC. Samples were analysed with a Waters 2D UPLC coupled to a Q Exactive mass spectrometer controlled by Xcalibur [18] software (v 2.3, Thermo Fisher Scientific, USA). A Waters ACQUITY UPLC BEH C18 column (1.7 μm, 2.1 mm × 100 mm) was used for chromatographic separation. A positive mode mobile phase consisted of 0.1% formic acid and acetonitrile, and a negative mode mobile phase consisted of 10 mM ammonium formate and acetonitrile.

Compound Discoverer software [19] (v3.1) was used to process the raw LC–MS/MS data that was gathered in both positive and negative ion modes. The processing included extracting peaks, retention time adjustment, additive ion pooling, filling in missing values and labelling background peaks. The metabolites were identified using the mzCloud [20] database (https://www.mzcloud.org/), ChemSpider [21] (data source from HMDB [22], KEGG [23], and LipidMaps [24]) databases (https://www.chemspider.com/), and the standard library built in-house by BGI. MetaX [25] was used for subsequent processing. To determine the relative peak area, the data were normalized using probabilistic quotient normalization (PQN). In order to correct the batch effect, a robust LOESS signal correction (QC-RLSC) based on quality control is used. All QC samples were calculated for their coefficient of variation (CV), and compounds with CVs greater than 30% were excluded.

Differential abundance (DA) score

The differential abundance (DA) score indicates a pathway’s ability to have greater amounts of metabolites than a control group. To compute the DA score, a nonparametric differential abundance test was first applied to all of the metabolites contained in a particular pathway. Following the assessment of metabolites that markedly altered their levels, the DA score was established in the following manner [26]:

$$\text{DA score}= \frac{\text{Number of metabolites increased}-\text{ Number of metabolites decreased}}{\text{Number of measured metabolites in pathway}}$$

The maximum number of the DA score is 1, and the minimum score is − 1. An increase in abundance of all metabolites in a pathway indicates a DA score of 1, while a decrease in abundance indicates a DA score of − 1.

Differential expression analysis

The R package ‘DEP’ [27] was utilized to identify differentially expressed proteins (DEPs). DEPs were defined based on a minimum |fold change| of 1.5 and an adjusted P value below 0.05. For the identification of differentially abundant metabolites, the R package ‘MetaboAnalystR’ [28] was used. The criteria were a minimum |fold change| of 1.2 or more, and an adjusted P value below 0.05. Differentially expressed genes (DEGs) were identified by the R package ‘edgeR’ [29], based on criteria of a |fold change| of 2 or more and an adjusted P value below 0.05.

Survival analysis

Survival outcomes, including overall survival and progression-free survival, were analysed using Kaplan–Meier curves and log-rank tests. Specifically, all survival analyses were conducted using the ‘survival’ [30] and ‘survminer’ [31] packages in R software. This approach ensured robust statistical evaluation of the survival data, with Kaplan–Meier curves providing visual representation of time-to-event distributions and log-rank tests enabling comparisons between groups. ‘Surv_cutpoint’ function of ‘survminer’ package in R software was used to calculate the cut-off value between high/low groups [32, 33].

Expression quantitative trait loci (eQTL) analysis

One hundred twenty patients from cohort 1 who had undergone WGS or WES in our previous studies were included in the eQTL analysis. VCFtools [34] and the R package ‘vcfR’ [35] were used to extract and integrate genotype information from these patients. The eQTL analysis of single nucleotide polymorphisms (SNPs) was performed using the R package ‘MatrixEQTL’ [36]. In general, eQTLs can be categorized into two main types: (1) cis-eQTLs, which refer to eQTLs located in close proximity to the regulated gene, within a 1 Mb region upstream or downstream of the gene; and (2) trans-eQTLs, which are eQTLs located further away from the regulated gene.

Somatic copy number alterations (CNAs) analysis

One hundred twenty patients from cohort 1 who had undergone WGS or WES in our previous studies were included in the CNAs analysis. CNAs were identified using BIC-seq2 [37], and the amplified genomic regions and the deleted genomic regions were identified using GISTIC2.0 (GenePattern) [38]. A default q value threshold of 0.25 of GISTIC2.0 was applied to determine significantly amplified regions or deleted regions, both at the focal levels and arm levels. Regarding the focal-level CNA analysis, G-scores were computed for significant genomic regions and gene-coding regions based on the frequency and magnitude of amplification or deletion impacting all genes.

Integrated analysis of the genome and transcriptome with the metabolome

The R package ‘multiOmicsViz’ [39] was used to display CNAs that influence mRNA and protein levels in either ‘cis’ (within the same aberrant locus) or ‘trans’ (remote locus) modes. All mRNA-CNA pairs for 18,155 genes and for all protein-CNA pairs for 4423 genes were used to compute spearman’s correlation coefficients and their corresponding adjusted P values for multiple tests. Recon3D, an escalating and expanding network of human metabolites reconstruction, was utilized to identify gene–metabolite pairs [40]. This dataset provided extensive information on biochemical reactions, including substrates, products, reaction reversibility, and related catalysing genes. The R package ‘cosmosR’ [41] was used to determine the interconnections between genes and metabolites in the Recon3D database.

In exploring the associations between somatic mutations and the prevalence of metabolites in ESCC, linear regression models (controlling for confounding factors) were applied following a previous study [42]. We adjusted for gender, tumour size, number of lymph node metastasis, age, and TNM stage to account for confounding effects. Each regression model included only one metabolite as a covariate along with the confounding factors. The formula used for the model was as follows:

$$logit\left(\pi \left(Y=1\right)\right)= \beta 0+\beta 1*(\text{abundance of metabolite})+\beta 2*(\text{gender})+\beta 3*(\text{tumour size})+\beta 4*(\text{number of lymph node metastasis})+\beta 5*(\text{age})+\beta 6*(\text{stage})$$

logit(π(Y = 1)): the mutation status of one gene; the range is from zero to one.

Only genes known to be related to cancer and mutated in at least 6% of our cohort were included in the analysis. A binary indication (1/0) was established for every single mutation variable. Then use t-statistics and associated P-values to assess the correlation. The P-values were adjusted using Benjamini–Hochberg false discovery rate (FDR) correction. The mutation features were ranked based on their associations with every single metabolite, allowing for comparisons based on statistical significance.

Identification of molecular subtypes and immune cell infiltration analysis

One hundred twenty-six tumour samples overlapping between metabolomics and proteomics datasets were used for molecular subtyping. Molecular subtyping based on integrated proteomic and metabolomic analysis was performed using similarity network fusion [43] (SNF). Proteomic and metabolomic data were pre-processed using the R package ‘SNFtool’ [43]. Then, the Euclidean distance matrix for the sample level of the two omics datasets was calculated. The distance matrix was used to create the sample similarity matrix network, and then a similar network fusion of the two networks was carried out. Finally, the fusion network of the samples in the two groups was obtained. Complex heatmaps showing the proteomic and metabolomic clustering results generated by the R package ‘ComplexHeatmap’ [44]. The relative infiltration expression of 28 kinds of immune cells in the TME was quantified by an enrichment score in single-sample gene set enrichment analysis (ssGSEA, R package ‘GSVA’ [45]). The tumour purity, immune score, and stromal score were computed by the estimate method [46] (R package ‘Estimate’).

Subtype diagnostic model for ESCC with machine learning methods

A total of 126 samples in molecular subtypes were divided into training and test cohorts at a 7:3 ratio. The least absolute shrinkage and selection operator (Lasso) method was used to select the most useful predictive features. Tenfold cross-validation was performed to identify the lambda value that produced the lowest test mean squared error (MSE). The predictive ability of ESCC subtypes was evaluated by the area under the ROC curve (AUC). The LASSO regression was performed by the R package ‘glmnet’ [47], and the R package ‘pROC’ [48] was used to determine the ROC curve and AUC.

Each predictor’s contribution in the LASSO regression model was defined by:

$$\text{Contribution }= \frac{|\text{coefficient of this feature}|}{\Sigma |\text{coefficient of this feature}|}$$

Multiplex immunofluorescence (MIF)

Multiplex immunofluorescence staining was performed via a 5-colour manual according to the manufacturer’s protocol. ESCC sections were incubated with primary antibodies after deparaffinization, hydration, antigen unmasking, and quenching. The following primary antibodies were used: rabbit monoclonal anti-CD4 (1:200, Cell Signaling Technology, #48274), rabbit monoclonal anti-CD8α (1:200, Cell Signaling Technology, #85336), rabbit monoclonal anti-CD56 (1:100, Cell Signaling Technology, #99746), and rabbit monoclonal anti-CD163 (1:400, Cell Signaling Technology, #93498). Spectral DAPI and TSA Opal 570, 620, 690, and 780 fluorophores were incubated with moisture at room temperature for 10 min. Seven samples from each of the S1 and S3 subtypes were randomly chosen for multiplex immunofluorescence analysis. The number of positive cells in various fluorescence channels was quantitated via ImageJ [49] software.

Immunohistochemistry (IHC)

The sections were deparaffinized, rehydrated, and boiled for 3 min in citrate buffer (pH 6.0). The slides were blocked with 5% goat serum for 1 h and incubated with primary antibody overnight. Then, the slides were incubated with a secondary antibody conjugated with HRP for 1 h (PV6001, PV6002, Zsbio), followed by incubation with 3% hydrogen peroxide solution for 15 min and DAB (ZLI9018, Zsbio) for 1 min. Aperio Image Scope [50] software was used for analysis.

Animal experiments

All mice were fed under specific pathogen-free conditions with unlimited access to food and water. The housing conditions included a 12-h light/dark cycle, a temperature range of 22–26 °C, and a humidity range of 45–65%. The animal experiments were approved by the Ethics Committee of Shanxi Medical University (SYDL2023040). The mice were given creatine-containing water (100 mM, 400 ml/time, once per day) or control water by oral gavage two weeks before the tumour cells were injected, and the treatments were continued until the termination of the studies. At the termination of the study, the mice were euthanized by intraperitoneal injection of sodium pentobarbital at a dosage of 150 mg/kg, followed by confirmation of death through cessation of heartbeat and respiration.

For the mouse xenograft model, the experiments were randomized. Six mice were set up for each different treatment condition to establish subcutaneous xenograft models. Subcutaneous tumours were established by injecting 4 × 10⁶ mEC25 cells combined with Matrigel into wild-type and Hk3 knockout C57BL/6 mice. A total of 2 × 10⁶ KYSE-150, TE-1, or mEC25 cells were implanted into NTG mice. One week after the tumour cells were injected, the mice were peritumorally injected with 300 mL of 100 mM creatine injection solution 3 times a week, while the control group mice were injected with the same volume of stroke-physiological saline solution in the same manner.

Flow cytometry

Subcutaneous tumour tissues were minced into small fragments with scalpels and fractionated into a single-cell suspension using an enzyme mixture consisting of collagenase II, collagenase IV, hyaluronidase II, and DNase I (Solarbio). Red blood cells were lysed with red blood cell lysis buffer (R1010, Solarbio), and tumour-infiltrating leukocytes were isolated with a percoll density gradient. Specific surface or intracellular antibodies were used to stain the cells in a cell staining buffer (E-CK-A107, Elabscience). The antibodies used were anti-CD45 Elab Fluor Red 780, anti-F4/80 PE-Cy7, anti-CD11b FITC, anti-MHC-II PE-Cy5, anti-CD86 PE and anti-CD206 PE (Elabscience). The cells were captured on a CytoFLEX flow cytometer (Beckman Coulter), and the resulting data were analysed with FlowJo V10 software [51].

Chemotaxis assays

Chemotaxis assays were performed using transwell plates (8 µm, Corning, Inc.). RAW264.7 cells with or without polarization stimulus were pretreated for three days under different conditions, after which 5 × 10⁴ RAW264.7 cells were placed in the upper chamber, while 1 × 10⁵ mEC25 cells were placed in the lower chamber. After 48 h of coculture, the non-migrated cells were removed, and the migrated cells were fixed and stained. All chemotaxis assays were independently repeated three times, with consistent trends observed across biological replicates. For each group, the number of migrated cells was quantified by counting cells in five non-overlapping 100 × magnification fields randomly selected to avoid positional bias. Data are presented as bar graphs showing the mean values, with error bars indicating the standard deviation (s.d).

RNA extraction and real-time PCR

Total RNA was extracted from RAW264.7 cell lines via an RNAiso plus kit (Takara, Dalian, China). Complementary DNA (cDNA) synthesis was performed using the PrimeScript® Reverse-Transcription reagent Kit with gDNA Eraser (Takara, Dalian, China), and real-time-qPCR was carried out using the TB Green® Premix Ex Taq® II Kit (Takara, Dalian, China). The sequences of primers used in this study are listed in Table S2 (Additional file 2). Gapdh was used for normalization, and calculations were performed via the 2-ΔΔCt method.

Cell migration and invasion assays

Migration and invasion assays were performed using transwell plates with a size of 8 μm (Corning, Inc.). For migration assays, a total of 50,000–70,000 KYSE-150 cells or mEC25 cells per well were seeded into the upper compartment and cultured in a 200 μL serum-free medium. The lower compartment was filled with 600 μL complete medium (10% FBS). After 24 or 36 h, the cells in the upper chamber were discarded. The cells that migrated through the membrane were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet dye. Transmigrated cells in five random fields per membrane were counted under 100 × magnification. For invasion assays, the upper chamber was precoated with BD Matrigel Matrix diluted 1:8 in serum-free medium, and the number of cells was increased to 100,000 cells per well. All other conditions were maintained as in the migration assays.

Targeted metabolomics

The LC–MS segment of the platform was established via a Shimadzu Nexera X2 LC-30 AD system coupled with an ACQUITY UPLC BEH Amide column (1.7 μm, 2.1 mm × 100 mm, Waters) and a triple quadrupole mass spectrometer (QTRAP 5500, AB SCIEX). Two-microliter samples were successively injected using an LC autosampler. The ACQUITY UPLC BEH amide column (1.7 μm, 2.1 mm × 100 mm, Waters) was thermostatic at 45 °C and flowed at a rate of 300 μL/min. Compound separation was carried out using a gradient consisting of 20 mM ammonium acetate and 5% acetonitrile at pH 9.45 (solvent A) in addition to 100% acetonitrile (solvent B). The gradient started at 5% solvent A for 1 min, gradually increased to 45% solvent A over 12 min, then to 60% solvent A over 1 min, followed by a 2-min retention before returning to the initial mixture for 0.1 min and a subsequent 3-min re-equilibration. During data acquisition, quality control samples were interspersed every six or eight samples.

The MS parameters were configured as follows: source temperature was set to 550 °C, ion source gas 1 (GAS1) was set to 40, ion source gas 2 (GAS2) was set to 50, curtain gas (CUR) was set to 35, and ion spray voltage floating (ISVF) was set to − 4500 V. The mass spectrometer operated with a dwell time of 200 ms. To create the metabolite MRM library, each metabolite standard (50 mg/mL) was initially analysed to determine the optimal MRM transition parameters. Next, the retention times of 40 energy-related metabolites were determined by measuring their corresponding MRM (Q1/Q3) transitions individually. A reference standard containing the 40 energy-related metabolites was prepared for LC–MS analysis by diluting a series of samples.

The raw MRM data files were processed via MultiQuant [52] software for peak identification, alignment, extraction, and filtration.

Statistics

For categorical variables versus categorical variables, Fisher’s exact test was used in a 2 × 2 contingency table; otherwise, the chi-square test was used. The Wilcoxon rank-sum test was used to identify the DEPs between normal adjacent tissues and tumours or between patients with different mutation statuses and CNAs. The Kruskal‒Wallis test was used to test whether protein abundance, metabolite abundance and gene expression were differentially expressed among the three subtypes or among the other subgroups. To account for multiple testing, the P values were adjusted via Benjamini‒Hochberg FDR correction. Overall survival and progression-free survival outcomes were depicted using Kaplan–Meier plots and log-rank tests. Adjustment of P-values was performed via the Benjamini–Hochberg FDR correction. Univariate Cox proportional hazards regression models were used to identify variables associated with overall survival. Benjamini–Hochberg false discovery rate was used to correct the P values. In both in vivo and in vitro experimental settings, statistical analyses to determine group differences were conducted employing a T-test or, where appropriate, an ANOVA. For chemotaxis assays, the results were presented as the mean ± standard deviation (s.d.). Spearman's rank correlation test was used to evaluate all correlations. Statistical analysis was performed via R [53] software (version 4.0.4) unless otherwise indicated. P values less than 0.05, 0.01, and 0.001 are indicated with *, **, and ***, respectively.

Results

Proteomic characteristics of ESCC tissues compared with those of normal adjacent tissues

This study selected surgically resected primary ESCC tissues and paired NATs from treatment-naive Chinese patients, with 127 samples used for proteomics and 136 samples used for metabolomics. Among these, 120 samples had previously been included in other multi-omics studies, encompassing genomics (80 samples analysed by WGS and 40 samples by WES), and transcriptomics analysis [13, 16]. A schematic of the experimental design and sample distribution is shown in Fig. 1A and Fig. S1A–C (Additional file 1). In summary, 126 tumour samples were subjected to both metabolomics and proteomics analyses, and 120 tumour samples were shared across all four omics (metabolomics, proteomics, transcriptomics, and genomics). The detailed clinicopathological information of all patients is summarized in Table S1 (Additional file 2).

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

For proteomic analysis, Fusion Lumos was used to acquire mass spectrometry (MS) data for 254 samples in DIA mode. The quality of the DIA data was evaluated on the basis of the intragroup coefficient of variation (CV), principal component analysis (PCA) and quantitative correlation of the samples. The quality control (QC) group, a mixture of all the samples, was intermittently injected to monitor the stability and reproducibility of the proteomics data. The median CV% for the QC samples was less than 20%, and PCA revealed distinct profiles between the QC, tumour and nontumour samples (Additional file 1: Fig. S2A, B). The average Spearman’s correlation coefficient among the control samples was 0.93, demonstrating the consistent stability of the MS platform. The correlations of tumour samples were between 0.62 and 0.96 (mean = 0.85) (Additional file 1: Fig. S2 C). The results of the quality control process indicated that our proteomic analysis system was robust. Additionally, off-column fractionation with data-dependent acquisition (DDA) was employed to generate spectral libraries for data extraction. The unique peptide distribution, protein mass distribution and protein coverage distribution in the spectral library are shown in Fig. S2D–F (Additional file 1), respectively.

Proteomics identified 6415 proteins in tumours and 6201 in NATs, with tumour samples showing higher proteomic complexity (Additional file 2: Table S3). On average, 4210 protein groups were identified per sample, ranging from a minimum of 3074 in NATs to a maximum of 5107 in tumours. Principal component analysis revealed that the nonneoplastic oesophageal tissues were aggregated together and obviously separated from the ESCC tissues (Fig. 1B), highlighting profound proteomic remodelling during carcinogenesis. A case-by-case analysis revealed that the number of proteins identified in the ESCC tissues was significantly greater than that identified in the NATs (Fig. 1C). We identified 2318 differentially expressed proteins (1696 upregulated, 622 downregulated; FDR < 0.05, fold change > 1.5) between ESCC and NATs (Additional file 2: Table S4). These proteins were predominantly localized to the cytoplasm (42.1%), nucleus (19.7%) and extracellular region (10.1%) (Additional file 1: Fig. S2G). Pathway enrichment analysis revealed that spliceosome, DNA replication, cell cycle, base excision repair, NOD-like receptor signalling pathway, lysosome, and phagosome, among others, were overrepresented among upregulated proteins. Downregulated proteins were enriched mainly in metabolic pathways, including branched-chain amino acid (BCAA) degradation, carbon metabolism, oxidative phosphorylation, and tyrosine metabolism (Fig. 1D). The representative differentially expressed proteins of the above pathways and the corresponding hazard ratios are shown in Fig. 1E. Key prognostic proteins included spliceosome-associated SNRPB (FDR = 0.025) and SF3A1 (FDR = 0.027), endoplasmic reticulum stress regulators EIF2S1 (FDR = 0.022) and DNAJC10 (FDR = 0.034), which were elevated in ESCC and linked to poor patient outcomes. Conversely, oxidative phosphorylation components COX6B1 (FDR = 0.019) and UQCRC1 (FDR = 0.024), along with BCAA degradation enzymes DLD (FDR = 0.028) and BCKDHA (FDR = 0.043), were reduced in ESCC and associated with favourable prognosis (Fig. 1F). Combined with the published proteomic data of Liu’s [54] and Li’s [55] studies in ESCC, only 93 and 218 proteins were consistently upregulated and downregulated, respectively, in the three datasets, suggesting that different proteomic research strategies might have different specificities and sensitivities (Additional file 1: Fig. S2H, Additional file 2: Table S4).

Metabolomic landscape of ESCC highlights amino acid reprogramming

Metabolomic analysis on 136 ESCC patients demonstrated robust platform reproducibility, evidenced by overlapping base-peak chromatograms (BPC), tight quality control (QC) clustering in PCA, and CV ratios exceeding 90% (Additional file 1: Fig. S2I–M). An extensive BLAST search revealed a total of 2418 annotated metabolites in the positive and negative ionization modes (Additional file 3: Table S5). These metabolites mainly included lipids (20.2%), amino acids and derivatives (16.3%), benzene and derivatives (11.7%), and carboxylic acids and derivatives (8.7%) (Fig. 2A). OPLS-DA clearly distinguished between ESCC and NATs (Additional file 1: Fig. S2N, O). A total of 1081 metabolites displayed differential abundance between ESCC and normal tissue samples (629 with increased abundance and 452 with decreased abundance in ESCC) (VIP ≥ 1, fold change > 1.2, FDR < 0.05) (Fig. 2B). Moreover, we performed KEGG pathway-based analysis of differentially abundant metabolites in which the differential abundance (DA) score was calculated to represent the tendency of increased/decreased metabolites in a pathway relative to the normal tissue. The analysis revealed that among the 59 metabolic pathways in which at least three metabolites were detected, 17 had DA scores greater than 0.5, and 8 had DA scores less than − 0.5 (Fig. 2C). Interestingly, the majority of enriched metabolic pathways in ESCC patients primarily related to amino acid metabolism, while most downregulated pathways in these patients were involved in carbohydrate metabolism (Fig. 2C). Notably, amino acid metabolism, including phenylalanine metabolism, tryptophan metabolism, and tyrosine metabolism, was significantly increased in the ESCC samples. The biosynthesis of unsaturated fatty acids, the pentose phosphate pathway and glutathione metabolism were also enriched among the upregulated metabolites. In contrast, the metabolites of galactose metabolism and glycolytic gluconeogenesis were significantly decreased in ESCC patients (FDR < 0.05) (Fig. 2D). A metabolic network map further delineated perturbations in glycolysis and amino acid synthesis/degradation pathways in ESCC (Fig. 2E).

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

Integrated multi-omics analyses of ESCC

Integrated transcriptomic-proteomic analysis revealed moderate mRNA-protein correlation in tumours (Spearman’s ρ = 0.399) but weak correlation in NATs (ρ = 0.145, Additional file 1: Fig. S3A, B). There were 5938 common genes identified in the transcriptome and proteome, which can be divided into strongly correlated and weakly correlated genes. Genes with strong correlations were enriched mainly in glutathione metabolism, pyrimidine metabolism and arginine and proline metabolism. Conversely, genes displaying weak or negative correlations were mainly involved in ribosome, endocytosis, and spliceosome pathways, which were mostly consistent in ESCC and NATs (Additional file 1: Fig. S3A, B). We subsequently compared the prognostic power of individual genes (log-rank test, P < 0.01, FDR < 0.05) in the two datasets (transcriptomic and proteomic data). Prognostic analysis identified 411 transcript-specific and 441 protein-specific survival-associated genes. Of these, 95 genes were common to both datasets, with 51 genes exhibiting opposite associations. Only 44 genes and their corresponding proteins showed similar prognostic trends at the mRNA and protein levels, suggesting the role of these genes in the prognostic monitoring of ESCC patients (Additional file 1: Fig. S3C, Additional file 3: Table S6).

We subsequently associated somatic driver mutations with mRNA and protein expression, observing cis-acting (acting on the gene within 1 million base pairs of the mutation) and trans-acting (acting on other genes outside) events of significantly mutated genes (SMGs) in ESCC [10,11,12]. We detected obvious cis effects of TP53 and NOTCH1 mutations on increasing the expression of these genes at both the RNA and protein levels, which is consistent with the known conclusion that TP53 mutation can increase protein stability [56]. Cis-effects of PIK3CA and FAT1 mutations were also observed, with significant increases and decreases in the expression of these proteins, respectively. Trans-effects included NFE2L2 mutation-associated spliceosome activation and CDKN2A mutation-associated cell cycle dysregulation, suggesting the importance of these driver mutations for pathway activation. Moreover, NOTCH1-altered samples exhibited elevated levels of associated proteins, including PSMA2, PSMA3, PSMA7, PSMB1, PSMB2, PSMB8, and PSMD10 (Fig. 3A).

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

We also examined the regulatory effects of 17,850 somatic copy number alterations (CNAs) in ESCC on mRNA and protein expression. CNAs can affect the expression of these genes positively or negatively in cis or trans mode, respectively. As demonstrated in Fig. 3B, the cis-effects of CNAs on mRNAs and proteins (diagonal lines) were obvious, and the strong trans-regulatory effects of CNAs (vertical stripes) at the RNA and protein levels were highly consistent. A total of 9,702 CNA cis-affected mRNAs (53%) and 954 CNA cis-affected proteins (22%) were observed, with 725 overlapping genes showing significant cis-effects across both omics analyses (Fig. 3C, Additional file 3: Table S7). Among these 725 genes, 270 DEGs at the protein level were significantly enriched in the spliceosome, N-glycan biosynthesis, arginine and proline metabolism, and fatty acid degradation pathways (Fig. 3D). CNA analysis using Gistic2.0 revealed the most frequent gains in chromosomes 1q, 3q, 5p, 7p, 7q, 8q, 14q, and 20q, as well as losses on chromosomes 3p, 4p, 4q, 5q, 9p, 9q, 11p, 11q, 13q, 18q, and 21q (Additional file 1: Fig. S3D). These findings were consistent with previous studies conducted on ESCC and TCGA cohorts [9, 12]. We detected amplifications of several driver oncogenes, including NFE2L2, MYC, FAM84B, ERBB2, and CCND1, and deletions of key tumour suppressors, such as FAT1, CDKN2A/B, PTEN, RB1, and KMT2D (Additional file 1: Fig. S3E). Patients with 1p34 or 13q22 gain alone exhibited a poor prognosis, while those with a combined gain at both 1p34 and 13q22 had an even worse prognosis (P = 0.0035, Fig. 3E–G). On chromosome 1p34, we identified significant cis-acting CNA genes that had a significant impact on the overall survival time of the patients (Fig. 3H). PPCS, a cis-acting gene at the protein level, encodes phosphopantothenoylcysteine synthetase, responsible for the CoA biosynthesis pathway by converting phosphopantothenate to phosphopantothenoylcysteine. Through further analysis focused on proteins significantly correlated with PPCS, we observed that the proteins—especially those positively related to PPCS—were mainly enriched in the proteasome pathway. Moreover, the protein expression of PPCS was positively correlated with the proteasome pathway score (Fig. 3I; Spearman’s ρ = 0.360, P = 3.266*e − 05). These associated proteins were predominantly unfavourable prognostic factors for ESCC patients.

In regard to connecting metabolites with genomic features, we assessed relationships between 18 frequently mutated genes in ESCCs [11, 12] and metabolite levels via a linear regression model (accounting for potential confounders) (Additional file 1: Fig. S4A, Additional file 3: Table S8). The pairs of genes/metabolites were identified via the Recon3D database [40]. The cancer-related altered genes with the highest mutation frequency (such as TP53, NOTCH1, and FAT1) showed minimal metabolite associations. Notably, EP300 mutations were inversely correlated with the abundance of prostaglandin J2 (PGJ2) (FDR = 0.016), while PTEN mutations were inversely correlated with L-threonine abundance (FDR = 0.071) (Additional file 1: Fig. S4A, B), which was consistent with the findings in the database.

When copy number alterations were examined, the connections between CNAs and metabolites were generally weak, with only a small number of ESCC-specific CNA peaks showing a relationship with metabolite levels (Additional file 1: Fig. S4C, Additional file 3: Table S9). For example, 13q14.11 chromosomal region (RB1 located) copy number deletion correlated with elevated prednisone (FDR = 0.016, Additional file 1: Fig. S4D). In addition, the 19p13.3 region deletion, which includes the gene encoding GAMT/DOT1L, mediates methyl transport in organisms. The results revealed a positive correlation between copy number deletion in the 19p13.3 region and the abundance of S-adenosylhomocysteine (SAH) (FDR = 0.028, Additional file 1: Fig. S4D). In conclusion, investigating the relationships between genetic traits and metabolites may shed light on the mechanisms underlying metabolic reprogramming in ESCC.

Molecular subtyping based on integrated proteomics and metabolomics analysis

Based on integrated proteomics and metabolomics analyses, we classified the observed three subtypes observed in ESCC as S1 (n = 47), S2 (n = 32), and S3 (n = 47), each with distinct molecular and clinical features (Fig. 4A). Notably, patients identified with the S3 subtype exhibited the poorest overall survival (OS) (log-rank test, P = 1.575*e − 6) and progression-free survival (PFS) (log-rank test, P = 1.074*e − 6) outcomes compared to those with the S1 and S2 subtypes (Fig. 4B). Evaluation of the clinical characteristics of the three subtypes revealed a later clinical stage in S3 patients, and the proportion of patients with lymph node metastasis in S3 patients was significantly greater than that in S1 and S2 patients (Fig. 4C). With respect to subtype-specific proteins and metabolites, gene set enrichment analysis revealed that metabolism-related pathways, especially amino acid metabolism-related pathways (such as glutathione metabolism, arginine and proline metabolism, and tyrosine metabolism), were obviously upregulated in S3 ESCC patients. In contrast, immune-related pathways (such as the T-cell receptor signalling pathway, antigen processing and presentation, and complement and coagulation cascades) were obviously downregulated in S3 patients, suggesting the regulation of metabolic remodelling in the immune microenvironment in these tumours. The ESCCs of subtypes S1 and S2 were both characterized by elevated levels of immune-related pathways, and S1 was also characterized by an increased number of protein processing pathways (such as lysosome and proteasome processes) (Fig. 4A). Consistent with pathway findings, ssGSEA confirmed reduced immune infiltration and elevated tumour purity in the S3 subtype compared with the other two subtypes (Fig. 4D). Given the significant differences in overall survival and progression-free survival between the S1 and S3 subtypes, and considering the limited availability of clinical samples, we focused on these two subtypes to conduct multiplex immunofluorescence (MIF) experiments to further validate the differences in immune cell infiltration. MIF validated diminished infiltration of CD8⁺ T cells (CD8⁺), CD4⁺ T cells (CD4⁺), and NK cells (CD56⁺) in the S3 subtype compared to the S1 subtype (Fig. 4E).

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

Subsequently, an examination of the subtype diagnostic signatures that exhibited potential clinical application was conducted. This analysis identified four distinct proteins (SEPTIN5, HK3, SCFD1, SSR1) and two metabolites (creatine and 2’-DG) as unique signatures. Notably, these signatures demonstrated favourable predictive performance, with an area under the curve (AUC) of 0.909 in the training set and 0.915 in the test set (Fig. 5A, B). The expression of these signatures in each subtype is depicted in Fig. 5C. The prognostic value of these signatures in the subtype diagnostic model and their associations with clinical features are depicted in Fig. S5 (Additional file 1). Additionally, the subtype diagnostic model was validated by assessing the abundance of six signatures in an independent ESCC cohort (n = 52, cohort 2) via immunohistochemistry (IHC) and mass spectrometry imaging (MSI). The model subsequently predicted 40 patients as subtype S1/S2 and 12 patients as subtype S3 (Additional file 3: Table S10). IHC revealed significantly higher intensity of SEPTIN5 (P = 0.002) and SCFD1 (P = 2.19*e − 4) proteins, and significantly weaker intensity of HK3 protein (P = 0.010) in the S3 subtype (Fig. 5D). MSI analysis revealed that creatine was significantly more abundant in S3 than in the other two subtypes (P = 4.18*e − 7, Fig. 5E, F). Although not statistically significant, the expression trends of the additional SSR1 and 2’DG signatures in the S3 subtype were consistent with those in the training set. Moreover, Kaplan–Meier curves in cohort 2 confirmed significantly worse overall survival in S3 subtype versus S1/S2 subtype patients (log-rank test, P = 3.48*e − 7), demonstrating robust prognostic efficiency of our subtype diagnostic model in ESCC (Fig. 5G). Further correlation analysis between these signatures revealed a significant negative correlation between hexokinase 3 (HK3) expression and creatine abundance (P = 4.15*e − 4, Fig. 5H), which was also confirmed in validation cohort 2 (P = 3.87*e − 3, Fig. 5I). In fact, the combination of HK3 and creatine also had good predictive performance for S3 subtype diagnosis (AUC = 0.866 in the training set, AUC = 0.821 in the test set) (Fig. 5J). As expected, patients with high creatine and low HK3 expression had worse survival than those with low creatine and high HK3 expression (Fig. 5K, L).

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

The role of the subtype diagnostic markers creatine and HK3 in mediating the immunosuppressive microenvironment

Creatine is a nitrogenous organic acid that is abundant in the human body. Exogenous creatine is obtained from diet uptake, and endogenous creatine is synthesized from arginine (arginine also participates in the urea cycle) in two steps: arginine to guanidinoacetate catalysed by L-arginine: glycine amidino-transferase (AGAT/GATM) and creatine production by guanidinoacetate N-methyltransferase (GAMT). Upon being transported into cells, creatine is converted to phosphocreatine by creatine kinase (CK) and is involved in ATP metabolism as an energy buffer system. Interestingly, S3 subtype patients exhibited upregulated arginine to creatine enzymes and metabolites but suppressed the urea cycle components (Fig. 6A). It has been reported that creatine reprograms macrophage polarization by suppressing IFN-γ–JAK–STAT1-driven M1 activation while promoting IL-4–STAT6-mediated M2 polarization via SWI/SNF chromatin remodelling [57]. Consistently, our MIF results demonstrated that the infiltration of CD163⁺ M2-like tumour-associated macrophages (TAMs) was greater in the S3 subtype versus the S1 subtype (P = 0.038, Fig. 4E). Additionally, the expression of several molecules related to creatine-remodelling macrophages also differed among the three subtypes (Additional file 1: Fig. S6A). Further analysis revealed that creatine metabolism components, including transporters, metabolites, and enzymes, were negatively correlated with the infiltration of most immune cells. However, the involvement of JAK1/2 in M1 macrophage polarization, as shown in Fig. S6A (Additional file 1), showed positively correlated with the infiltration of most immune cells (Fig. 6A, B). Proteomics revealed creatine inversely correlated with complement activation and innate immune pathways (Spearman’s ρ > 0.3, FDR < 0.05), suggesting its role in suppressing antitumour immunity, which was consistent with the observation that a higher creatine content in the S3 subtype was associated with less immune cell infiltration (Additional file 1: Fig. S6B). In line with the poor prognosis of S3 subtype patients with increased creatine expression, many proteins and metabolites involved in creatine metabolism and macrophage polarization were also significantly and consistently related to the prognosis of ESCC patients (Additional file 1: Fig. S6C). The correlations between several stimulatory immune checkpoint molecules and creatine and the corresponding hazard ratios are shown in Fig. S6D, E (Additional file 1). More interestingly, the greater the number of lymph node metastases and the later the clinical stage, the greater the creatine expression (Additional file 1: Fig. S7A, B). Patients with higher creatine levels had lower immune scores and stromal scores (Additional file 1: Fig. S7C). Lymph node metastatic progression intensified the negative correlation between creatine and immune cell infiltration, and the correlation with immune checkpoints also increased (Additional file 1: Fig. S7D, E). MIF confirmed that creatine accumulation inversely correlated with cytotoxic immune cells (CD8⁺ T cells, CD4⁺ T cells and NK cells) but positively associated with M2-like TAMs (CD163⁺, Additional file 1: Fig. S7F), linking metabolic shifts to immune evasion. Further in vitro and in vivo studies demonstrated that creatine treatment at different concentrations had no effect on the proliferation, invasion or migration of mouse- or human-derived ESCC cells (Additional file 1: Fig. S8A–C). A high-creatine diet also failed to promote the growth of ESCC in subcutaneous xenograft models in severe immunodeficiency mice (Additional file 1: Fig. S8D–F).

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

Given that creatine and HK3 are two effective predictive markers for S3 subtypes and that they are significantly negatively correlated, we further explored the role of HK3 in tumour immunity and its relationship with creatine. Hexokinase 3 (HK3) is one of four hexokinase isoenzymes. HK3 is the first rate-limiting enzyme in the glycolysis process and is responsible for the phosphorylation of glucose to 6-phosphate glucose in the presence of ATP. HK3 demonstrated significant positive associations with immune score and immune cell infiltration across a wide range of proteomic (P = 3.37*e − 5) and transcriptomic data (P < 0.0001) for ESCC patients (Fig. 6C, D). Similarly, proteomics identified proteins, positively correlated with HK3 expression (Spearman’s ρ > 0.3, P < 0.05, FDR < 0.1), and were obviously enriched in immune response pathways—including T-cell receptor signalling pathway, Fc-gamma receptor signalling pathway, innate immune response, positive regulation of T-cell proliferation and IFN-γ-mediated signalling pathway—aligning with the immune-depleted phenotype in S3 subtype (Fig. 6E). Prognostically favourable HK3-related immune regulators were elevated in S1 but suppressed in S3 subtypes (Fig. 6F).

Next, we explored the effects of HK3 and creatine on the tumour microenvironment and investigated whether there was a synergistic effect. We fed wild-type and Hk3 knockout C57BL/6 mice a standard or high creatine diet. After two weeks of feeding, mEC25 cells were subcutaneously injected into the flanks of immunocompetent C57BL/6 mice. After continued feeding for 3–4 weeks, the spleen and tumour tissue of the mice were collected for flow cytometry analysis (Fig. 7A). We found that the tumour volume in the high-creatine diet and HK3 knockout alone groups was significantly greater than that in the control group. Moreover, Hk3 knockout in combination with a high-creatine diet synergistically promoted tumour growth (Fig. 7B, C). The ratio of pan-TAMs, which was defined as Cd45⁺F4/80⁺Cd11b⁺, improved upon high-creatine diet/Hk3 knockout treatment or in combination (Fig. 7D). We further observed that the effects of a high-creatine diet/HK3 knockout on increasing the percentage of M2-like TAMs (Cd45⁺F4/80⁺Cd11b⁺Cd206⁺ cells) and decreasing the percentage of M1-like TAMs (Cd45⁺F4/80⁺Cd11b⁺Cd86⁺Mhc-II⁺ cells) were strengthened by combined treatment (Fig. 7E, F). These results were consistent with the poor prognosis of S3-subtype ESCC patients and indicated that creatine and HK3 might affect the tumour microenvironment in a synergistic manner. Further study of 126 ESCC patients revealed that the creatine content in the HK3 high-expression group was significantly lower than that in the HK3 low-expression group (Additional file 1: Fig. S8G). In vitro experiments revealed that creatine treatment and Hk3 knockdown synergistically promoted the chemotaxis of RAW264.7 macrophages and IL- 4-induced M2-like RAW264.7 cells cocultured with mEC25 cancer cells (Fig. 7G, H). Additionally, creatine treatment and Hk3 knockdown synergistically promoted M2-like TAM polarization and the secretion of associated effector factors such as ARG1, IL-10 and TGF-β (Fig. 7I). These results suggested that creatine and HK3 might mediate the immunosuppressive microenvironment of ESCC by promoting M2 macrophage polarization in the S3 ESCC subtype.

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

Both creatine and HK3 are crucial regulators of energy metabolism. To investigate how creatine metabolism and HK3-mediated glycolysis reshape energy metabolism to adapt to the microenvironment, we further investigated the effects of Hk3 knockdown and creatine treatment on energy metabolism in RAW264.7 cells via targeted metabolomics analysis. Compared to the normal control, Hk3 knockdown reduced glycolysis derived from glucose-6-phosphate, pentose phosphate metabolism activity, and the tricarboxylic acid (TCA) cycle, reflecting Hk3's role as a hexokinase. Following Hk3 knockdown, supplementation with creatine significantly enhanced the TCA cycle and pentose phosphate metabolism, while glycolysis was suppressed. The increased synthesis of UDP-glucose and increased pentose phosphate pathway activity upon creatine supplementation led to greater production of NADPH (Additional file 1: Fig. S9A, B; Additional file 3: Table S11). It has been reported that an increase in NADPH inhibits intracellular ROS production. Previous study indicates that an increase in NADPH reduces intracellular ROS production [58], which is significantly associated with decreased M1-like polarization in macrophages. The increase in metabolites related to the TCA cycle could imply a key role for energy sources other than glycolysis in M2 macrophage polarization. Interestingly, supplementation with creatine also caused an increase in intracellular lactate levels, hinting at its potential effects on lactate metabolism. The metabolic reprogramming and energy transfer caused by the downregulation of Hk3 combined with creatine supplementation may contribute to driving macrophages toward M2-like polarization. M2 polarization of TAMs may contribute to the formation of an immunosuppressive microenvironment and a poor prognosis for patients through inhibiting the proliferation of CD8⁺ effector T lymphocytes, natural killer cells, etc. [59, 60]. In conclusion, our study highlights the pivotal roles of creatine and HK3 in mediating the immunosuppressive microenvironment, suggesting that targeting creatine metabolism might be an effective target for enhancing immunotherapy and targeted therapy for ESCC.

Discussion

Despite immunotherapy advances, ESCC remains greatly challenging due to its complex heterogeneity, which leads to varying therapeutic responses and drug resistance. Unlike other malignancies, ESCC lacks well-defined molecular subtypes to direct treatment strategies [61,62,63]. Molecular subtyping on the basis of microenvironment heterogeneity holds particular promise for predicting immunotherapy responders and identifying novel therapeutic targets. Recent studies have reported molecular subtyping of ESCC. For example, an analysis of whole-genome sequencing data from 508 ESCC patients revealed three prognostically significant molecular subtypes: NFE2L2 mutated, RTK-RAS-MYC amplified, and double negative [12]. Liu et al. defined ESCC into two molecular subtypes based on proteomic analysis and provided a potential therapeutic outlook for ESCC [55]. Multi-omics analyses, including genomic, epigenomic, transcriptomic, and proteomic analyses, have classified ESCCs into four subtypes: cell cycle pathway activation, NRF2 oncogenic activation, immune suppression (IS) and immune modulation (IM). A classifier with 28 features was developed to identify the IM subtype that could predict the response to anti-PD-1 therapy [13]. Nonetheless, current subtypes fail to encompass all patients and lack integration of metabolomics data, which is critical for understanding the impacts of metabolic remodelling on the immune microenvironment. Our study addresses these gaps through paired metabolomic and proteomic profiling of a large ESCC cohort, revealing three distinct subtypes (S1, S2, and S3) with prognostic significance. The S3 subtype demonstrated the worst clinical outcomes and an immunosuppressive microenvironment shaped by synergistic creatine accumulation and HK3 deficiency-driven macrophage polarization.

The immunosuppressive tumour microenvironment remains a major barrier to immunotherapy success [64]. Its development involves intricate interactions among diverse immunosuppressive cells, cytokines, and signalling molecules. Major immunosuppressive cells include myeloid-derived suppressor cells (MDSCs), regulatory T cells (Tregs), TAMs, tumour-associated neutrophils (TANs), and DCs [65]. Emerging strategies focus on microenvironment remodelling to enhance treatment efficacy. For instance, Gao et al. introduced an anti-PD-1 antibody-conjugated nanoplatforms to mitigate the immunosuppressive effects on Tregs and DCs [66]; Hu et al. recommended employing radiotherapy and PARP inhibition to reverse the immunosuppressive microenvironment caused by IDH1 mutations, thereby increasing the efficacy of immune checkpoint inhibitors [67]. The intersection of metabolic reprogramming and immunosuppression presents particularly promising therapeutic opportunities [68]. Additionally, some drugs target cancer metabolisms have been shown to assist in tumour immunotherapy. Clinical validation of metabolism-targeting agents (e.g., glutaminase-1 inhibitors with PD-1 blockade [69], IDO inhibitors in melanoma [70, 71]) underscores this potential. Given the scarcity of targeted metabolic drugs for ESCC, exploring the connection between metabolic changes and immune cell functionality, and employing metabolomics-based patient stratification, can potentially facilitate greater benefits from immunotherapy and other targeted treatments for ESCC patients. In this study, we unveiled the remodelling effect of creatine and HK3 on the immunosuppressive microenvironment and propose creatine metabolism targeting as a novel strategy to potentiate immunotherapy for ESCC.

Creatine in mammals is transported into the cell by SLC6A8 and is involved in ATP metabolism as an energy buffer system by creatine kinases. While the classic role of creatine metabolism in energy reserve and energy buffer is established, its two-sided role in cancer biology remains underexplored [72]. Emerging evidence suggests creatine enhances metastatic adaptability in distant organs [73], GATM-driven Smad2/3 activation promoting colorectal cancer liver metastasis [74], and stromal-induced creatine synthesis fuelling phosphocreatine/ATP production for pancreatic cancer progression [75]. Conversely, adipocyte creatine inhibition suppresses obesity-associated breast cancer growth [76]. Beyond tumour cells, creatine modulates immune landscapes to regulate tumour growth. The creatine-mediated ATP buffering system affects the antitumour activity of CD8⁺ T cells in melanoma and colon cancer [77]. CKMT1A-mediated creatine metabolism correlates inversely with CD8⁺ T cell infiltration and fibroblast abundance via glycolytic/gluconeogenesis pathway modulation [78]. In this study, we found that creatine was highly expressed in the S3 subtype, which has low infiltration of CD8⁺ T cells, CD4⁺ T cells, and NK cells and high infiltration of M2-type TAMs. Similarly, Ji et al. revealed that creatine reprogrammed macrophage polarization by suppressing M1 (IFN-γ) but promoting M2 (IL-4) effector functions, which was attributed to the role of suppressing IFN-γ-JAK-STAT1 transcription factor signalling while supporting IL-4-STAT6-activated arginase 1 expression [57]. M2-like macrophages participate in tumour immune escape by inhibiting the proliferation of cytotoxic T lymphocytes (CTLs), NK cells and other effector immune cells, resulting in poor patient prognosis [59, 60]. Therefore, creatine may be an important metabolic modulator of the tumour immune microenvironment, but the relevant mechanism is in its initial stage.

Notably, creatine and hexokinase 3 (HK3) emerged as co-predictors of the S3 subtype, exhibiting an inverse correlation in multiple cohorts. HK3 is the first rate-limiting enzyme in glycolysis and is responsible for the phosphorylation of glucose to 6-phospho-glucose in the presence of ATP. In liver cancer and cervical cancer, HK3 plays a protective role by reducing oxidative stress and increasing ATP levels [79]. In non-small cell lung cancer, HK3 is correlated with immune infiltrates and predicts the response to immunotherapy [80]. In clear cell renal cell carcinoma, HK3 is highly correlated with the abundance of immune cells and specifically stimulates the infiltration of monocytes/macrophages, promoting malignant biological processes and immune escape [81]. Mechanistically, macrophage polarization is regulated by metabolites involved in glycolysis, the TCA cycle, and fatty acid metabolism [82]. In this study, our experiments reveal that HK3 knockdown synergizes with creatine supplementation to promote tumour growth and macrophage M2-like polarization. Both HK3 and creatine are important molecules involved in ATP metabolism. In colorectal cancer, hexokinase and mitochondrial creatine kinase can bind to voltage-dependent anion channels (VDACs) in the outer membrane of mitochondria, thereby regulating the distribution of high-energy phosphate and reshaping energy metabolism to adapt to microenvironmental changes [83]. In glucose-deprived p53 wild-type mice, GAMT and creatine levels are significantly increased [84]. uMtCK enhances hexokinase-dependent glycolysis via the JNKMAPK/JUN signalling pathway and promotes cancer progression in gastric cancer [85]. In this study, metabolic analyses of RAW264.7 cells with Hk3 knockdown revealed decreases in glycolysis, the TCA cycle, and ATP concentrations. However, creatine supplementation significantly enhanced mitochondrial energy metabolism and the pentose phosphate pathway without affecting ATP generation. Combined Hk3 knockdown and creatine supplementation collectively induced alterations in the energy metabolism of macrophages toward M2 polarization. Therefore, creatine metabolism and glycolysis may jointly regulate the transfer network of energy metabolism that affects tumour-associated macrophage polarization, thereby affecting the immune microenvironment and antitumour response of other immune cells. While mechanistic details of creatine-HK3 crosstalk require further study, our data nominate creatine depletion or HK3 activation as actionable strategies to reprogram immunosuppressive microenvironment in S3 patients.

In addition, the current study broadened our understanding of the molecular characteristics of ESCC on the basis of metabolomics and proteomics combined with our previous genomics and transcriptomics data. Our proteomic landscape revealed the consistency between this study and other studies to some extent; for example, the spliceosome, DNA replication, cell cycle and lysosome pathways were upregulated, whereas metabolic pathways such as tyrosine metabolism pathways were downregulated in ESCC. However, some proteins were not identified in this study or had inconsistent expression patterns, which might be attributed to different proteomic research strategies with different specificities and sensitivities. Our metabolomics data revealed that tryptophan metabolism, tyrosine metabolism, and branched-chain amino acid metabolism were significantly associated with ESCC.

Conclusions

Our integrated analysis reveals molecular variation and defines ESCC subtypes with distinct metabolic-immune features, offering actionable targets to disrupt immunosuppression. Creatine accumulation and HK3 protein deficiency synergistically promoted M2 polarization of macrophages and participated in the formation of an immunosuppressive microenvironment. Targeted creatine metabolism might be an effective target for enhancing immunotherapy and targeted therapy for ESCC.