Integrated multi-omics reveal lactate metabolism-related gene signatures and PYGL in predicting HNSCC prognosis and immunotherapy efficacy

Abstract

Background

Head and neck squamous cell carcinoma (HNSCC) treatment faces significant clinical challenges. Lactate metabolism plays a crucial role in the initiation of many cancers and the tumor microenvironment (TME). However, the prognostic significance of lactate metabolism-related genes (LMRGs) and the role of TME in HNSCC require further elucidation.

Methods

We built a prognostic multigene signature with LMRGs and systematically correlated the risk signature with immunological characteristics and immunotherapy efficacy. Next, a series of single-cell sequencing analyses were used to characterize lactate metabolism in TME. Finally, single-cell sequencing analysis, immunofluorescence analyses, and a series of in vitro experiments were used to explore the role of PYGL in HNSCC. Potential drugs targeting PYGL were screened using AutoDock 4.2.

Results

A prognostic multigene signature based on LMRGs was developed, which effectively stratified patients into high- and low-risk groups, with significant differences in overall survival (OS) and progression-free survival (PFS). Patients in the low-risk group exhibited reduced lactate metabolism, higher CD8 + T cell infiltration, and improved response to immunotherapy. Single-cell sequencing revealed that tumor cells had the most active lactate metabolism compared to other cells in the TME. PYGL, identified as the most critical prognostic gene, was highly expressed in tumor-associated macrophages and played a role in inhibiting M1 macrophage polarization. Knockdown of PYGL led to reduced lactate levels, and its expression was inversely correlated with CD8 + T cell infiltration. Furthermore, PYGL was involved in copper-dependent cell death, highlighting its potential as a therapeutic target. Drug screening identified elesclomol, which showed promising results in PYGL-knockdown cells.

Conclusions

The study established a robust LMRGs-based prognostic model that not only predicts patient survival but also correlates with the immune microenvironment in HNSCC. PYGL emerged as a key biomarker with significant implications for both prognosis and therapeutic intervention. Its role in regulating lactate metabolism and immune suppression suggests that targeting PYGL could enhance the efficacy of immunotherapies. This research provides a foundation for future clinical strategies aimed at improving outcomes in HNSCC by modulating the tumor’s metabolic and immune landscapes.

Background

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common malignancy worldwide [1]. There are about 600,000 new HNSCC patients and 350,000 cancer-related deaths each year [2]. Despite tremendous advancements in medical technologies for early detection of HNSCC, over half of all patients are detected after the disease has progressed to an advanced phase [3]. Although the efficacy of HNSCC treatments (surgery, radiation, and chemotherapy) has improved in recent years, the 5-year survival rate for HNSCC patients remains poor, especially for those with late-stage disease [4, 5]. The discovery of disease-related prognostic biological markers in HNSCC is therefore essential in enhancing the likelihood of early detection and the development of innovative treatment approaches to improve the survival probability of patients.

Lactate, a byproduct of anaerobic glycolysis, has traditionally been considered merely a waste metabolite or endpoint in cancer metabolism [6]. Nonetheless, this assumption has been challenged since new research suggests that lactate is a crucial player in a variety of carcinogenic mechanisms, including immunosuppression, metabolism, angiogenesis, and metastasis [7]. Research indicates that the tumor microenvironment (TME) significantly influences immune responsiveness to malignancies [8]. As cancer advances, nutrients and oxygen required for anti-tumor immunological cell activity are exhausted, and by-products of tumor metabolic activity including adenosine and lactate concentrations in the TME, repress the anti-tumor immune reaction [9]. A 2019 meta-analysis found that serum LDH levels can predict clinical response to immune checkpoint inhibitors (ICIs) in non-small cell lung cancer patients, with elevated pre-treatment LDH levels associated with significantly shorter progression-free survival (PFS) and overall survival (OS) durations [10]. Lactate production, on the other hand, is a complex process involving numerous processes. Hundreds of cellular biological changes are implicated in this process, which is remarkably controlled by a group of genes known as lactate metabolism-related genes (LMRGs) [11, 12]. Consequently, as compared to a single gene, a model that incorporates numerous LMRGs that are suggested to perform essential functions in HNSCC could improve the prediction performance of the prognosis for patients with HNSCC. These critical LMRGs might possibly be utilized as therapeutic targets for treating HNSCC.

This study systematically evaluated the correlations between LMRGs, prognosis, and TME in HNSCC. A scoring system was developed using these genes to predict OS and immunotherapy efficacy. The predictive performance of the risk model was validated using an integrated Gene-Expression Omnibus (GEO) cohort of 97 patients and an independent FuJianZhongLiu (FJZL) cohort of 192 patients. To understand the long-term impact of LMRGs on TME, the associations between risk model, immune cell infiltration, and immunotherapy were meticulously investigated. Single-cell sequencing data characterized lactate metabolism across various cell types in the TME. Through in vitro experiments and sequencing analysis, we examined the regulation of lactate metabolism and the immune microenvironment of PYGL, which is crucially enriched in macrophages and CD8 + T cells.

Methods

Data collection

We obtained mRNA expression profiles and clinical data for 501 HNSCC patients from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). An independent validation dataset was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/gds) under the accession number GSE41613. An additional transcriptomic validation dataset was obtained from 192 patients diagnosed with nasopharyngeal carcinoma, a subtype of head and neck squamous cell carcinoma, treated at Fujian Cancer Hospital between January 2015 and January 2018. Eligible participants were those aged 18 or older, with normal blood, kidney, and liver function, no history of other malignancies, and who underwent standard radiotherapy or chemoradiotherapy. The study was approved by the Ethics Committee of Fujian Cancer Hospital (Fuzhou, China; Approval No. K2022-084-01). Each patient provided written informed consent for participation in this study. Nasopharyngeal carcinoma tissue samples were collected via biopsy and preserved in liquid nitrogen for future RNA extraction. Extracted RNA samples were subsequently sent to Shenzhen Huada Gene Technology Co. for transcriptome sequencing (RNA-seq). To eliminate batch effects, we utilized the R function “sva” [13]. The clinical proteomic tumor analysis consortium (CPTAC) databases provided information on how PYGL was expressed translationally in tumor and normal samples. Additionally, 244 LMRGs were acquired from the Molecular Signatures Database (MSigDB) by obtaining “Lactate“ (Table S1) [14].

Differentially expressed LMRGs and establishment of prognostic signature predicated

We employed the “limma” package in R [15] to evaluate the expression data and identify genes that were differentially expressed between tumor and normal samples. The analysis criteria were set with a log2 fold change (FC)| >1 and an adjusted p-value < 0.05. Subsequently, we constructed a prognostic prediction model in the training cohort using the least absolute shrinkage and selection operator (LASSO), based on univariate Cox analysis [16]. Each HNSCC patient’s risk score is predicted using the equation given below: $risk\,score\, = \,\sum\nolimits_{i = 1}^n {\left( {exp\, \times coef} \right)}$, where exp signified the gene’s expression level, and coef signified the gene’s coefficient (Data S2). For the training and validation, each patient’s risk score was computed. In this case, the threshold level was determined to be the median risk score.

Gene set enrichment analysis

Gene Set Enrichment Analysis (GSEA) was performed on the MSigDB collections (c2.cp.kegg, c5.go.bp.v7.2.symbols.gmt, and h.all.v6.2.symbols.gmt) to evaluate whether the specified gene sets showed statistically significant differential expression between high-risk and low-risk groups. We used the normalized enrichment score (NES) and false discovery rate (FDR) to categorize the gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and hallmark pathways enriched in different phenotypes. Additionally, a validated set of 31 genes related to cell cycle progression (CCP) was used to estimate cell proliferation rates [17]. Subsequently, we selected six energy metabolism-related pathways, including hypoxia, lactate, glycolysis, lipid metabolism, pentose phosphate pathway, and oxidative phosphorylation, to explore their spatial activities in our HNSCC samples. The corresponding references for these gene signatures were obtained from previous studies [18].

Analyses of the TME’s immunological properties

By evaluating tumor samples utilizing the R program “ESTIMATE,” we successfully determined their tumor purity, immune, and estimate scores. Charoentong et al.‘s research provided initial data on 80 immunomodulators, including chemokines, receptors, and MHC (Table S2) [19]. Seven steps form the tumor immunity cycle, which represents the anticancer immune response (Table S3) [20]. Based on the individual specimens’ gene expressions, Xu et al. utilized single sample gene set enrichment analysis (ssGSEA) [21] to examine these steps. Subsequently, various algorithms were developed to calculate TIIC infiltration levels in the TME using TCGA RNA-seq data (Table S4) [17, 21,22,23,24,25]. In addition, we utilized prior investigations to identify the TIICs’ effector genes (Table S5) [26].

Next, using data from Auslander’s work, we obtained 10 inhibitor immune checkpoints molecules with therapeutic promise [27]. Gene sets indicating T cell-inflamed gene expression profile (GEP) or immunological cytolytic activity (CYT) were sourced from earlier studies [28, 29]. Finally, we retrieved T-cell receptor (TCR) and B-cell receptor (BCR) Shannon entropy data from the work of Vésteinn Thorsson et al. [30].

To elucidate the relationship between immunological signatures and the effectiveness of immunotherapy, we obtained the Tumor Immune Dysfunction and Exclusion (TIDE) score, T cell exclusion score, and T cell dysfunction score from the online website (http://tide.dfci.harvard.edu) [31]. We obtained patient immunophenoscores (IPS) from The Cancer Immunome Atlas. Patients’ responses to immunotherapy have been predicted using IPS [19]. Additionally, we gathered multiple immunotherapy cohorts (melanoma-GSE100797, melanoma-PRJEB23709, and melanoma-GSE91061) from the tumor immunotherapy gene expression resource (TIGER) website (http://tiger.canceromics.org).

Machine learning analysis

We employed the R package “ranger,” a weighted variant of the random forest method, to analyze the impact of selected gene expression on OS in patients and to calculate each gene’s variable importance score (VIS). Sliding window sequential forward feature selection (SWSFS) method identified the most significant genes in the risk signature [32].

Single-cell level analysis and spatial transcriptome analysis

The HNSCC single-cell transcriptome data (GSE103322, GSE139324) were obtained from the GEO database [33, 34]. The expression matrix was normalized by “NormalizeData” of the “Seurat” package. Using “FindVariableFeatures”, we identified the top 2,000 highly variable genes. Cell annotation followed a previous study [33]. Using the same formula, each cell from GSE103322 was assigned a risk score. The software CellChat v1.1.3 analyzed ligand-receptor interactions to determine whether cell-cell communication took place. Cell groups with fewer than 20 cells were excluded from the analysis. Pairwise tests were used to assessed the statistical significance of communication probability values. Transcription factors (TFs) were predicted from the 400 cells using the SCENIC package [35]. A volcano plot was used to visualize the TF activity scores between the high and low risk score strata. Tumor immune single-cell hub (TISCH, http://tisch.comp-genomics.org/) was utilized to assess PYGL expression [39]. Additionally, the SpatialDB online tool (https://www.spatialomics.org/SpatialDB/) [36], a database for spatially resolved transcriptomes, was employed to analyze the spatial expression and co-localization of PYGL and macrophage marker CD68 in breast cancer, prostate cancer, and melanoma.

Cell culture and cell transfection

The HNSCC cell line murine squamous cell carcinoma cell (SCC7) was kindly provided by Prof.Zhu from the Shanghai Ninth Peoples Hospital (Shanghai, China). THP1 monocytes were purchased from Wuhan Saiweier Biotechnology Co., Ltd. PYGL shRNA knockdown cells were obtained by lentivirus infection, which was synthesised by Shanghai Jikai Company. The Pygl-targeting siRNAs (siPygl) and a negative control siRNA (siNC) synthesized by OBIO Technology (Shanghai, China). Lactate was obtained from Sigma-Aldrich (St. Louis, MO, USA).

Macrophage polarization

THP-1 cells were differentiated into M0 macrophages using 50 ng/ml PMA for 48 h, followed by stimulation with 100 ng/ml LPS and 20 ng/ml IFN-γ for an additional 48 h to induce M1 polarization. M0 macrophages were stimulated with 20 ng/ml IL-4 for 48 h to induce M2 polarization. RAW267.4 cells were stimulated with 100 ng/ml LPS for 48 h to induce M1 polarization and stimulated with 20 ng/ml IL-4 for 48 h to induce M2 polarization.

Clinical samples and multiplex fluorescent immunohistochemical (mfIHC)

To examine the link between PYGL expression and CD8 + T lymphocytes in the HNSCC milieu, we utilized a tissue microarray (HoraC060PG01) of 48 oral cancer (OC) patients provided by Shanghai Outdo Biotechnology. OC accounted for about 40% of HNSCC, which is a common cancer in HNSCC. The Ethics Committee of Fujian Provincial Tumor Hospital approved the use of OC tissues for this study. The mfIHC was prepared as described in a previous study [37]. Primary antibodies were used as follows: anti-PYGL rabbit antibody (15851-1-AP, Proteintech, Wuhan, China), anti-PD-L1 rabbit antibody (ab205921, Abcam, US), and anti-CD8 rabbit antibody (66868-1-IG, Proteintech, Wuhan, China). Each antibody was treated using the Opal TSA fluorophores. Lactate levels were measured using pan anti-Kla (PTM-1401) by immunohistochemical (IHC).

Chemotherapy sensitivity analysis and molecular Docking study

Using Pearson’s correlation analysis, we examined the relationship between PYGL expression and drug sensitivity, which was obtained through the CellMiner interface (https://discover.nci.nih.gov/cellminer) [38]. The 3D structure of PYGL was sourced from RCSB PDB (https://www.rcsb.org/). The binding site of PYGL was then identified using Binding Site Detect in AutoDock 4.2 [39] and virtual screening and docking were performed to identify potential drugs targeting PYGL.

Statistical analysis

All statistical data analyses shown in the figures were performed using R (version 3.6.0). The Wilcoxon rank-sum test assessed differences in continuous variables between the two groups, while the chi-square test compared categorical data. The log-rank test was conducted to assess the prognostic significance of categorical factors. Statistical significance was determined for all analyses if the two-paired p-value was less than 0.05. Statistical significance was defined as *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.

Fig. 1

Flowchat

Results

Identification of differentially expressed LMRGs and construction of the prognosis prediction model

Figure 1 presents a schematic workflow of our study. The TCGA identified 5968 differentially expressed genes (DEGs) in 501 tumor tissues and 44 adjacent non-tumor tissues. By overlapping the DEGs and LMRGs in HNSCC, we identified 23 LMRGs that were differentially expressed (Fig. 2A-B). We utilized univariate Cox regression and Lasso regression analysis on 16 LMRGs to construct an optimal risk model (Supplementary Figs. 1 A-B). The HNSCC specimens were divided into high-risk and low-risk groups according to their median risk score. Kaplan-Meier curve analysis revealed significantly poorer OS and PFS in the high-risk group (Fig. 2C, Supplementary Figs. 1 C). The GEO (GSE032062) database and the FJZL cohort were used to validate these findings (Fig. 2D-E). ROC analysis of TCGA data demonstrated that the risk signature reliably predicted survival, with AUC values of 0.65, 0.68, and 0.70 for one, three, and five years, respectively (Fig. 2F). Similar AUC values were observed in the GEO and FJZL cohorts. Univariable and multivariable Cox analyses demonstrated an independent association between the high-risk signature (HR = 3.929; 95% CI = 2.290–6.741; P < 0.001) and poor OS (Supplementary Figs. 1D-E). A nomogram was created to estimate 1-, 3-, and 5-year OS rates by incorporating the risk score, age, and clinicopathological parameters (Fig. 2G). The constructed nomogram’s performance in the TCGA-HNSCC cohort closely matched that of an ideal model (Fig. 2H). The ROC curves of the nomogram for predicting 3-year and 5-year survival rates were 0.726 and 0.718, respectively (Supplementary Figs. 1 F).

Fig. 2

Construction of LMRGs prognostic model and survival analysis. (A) The Venn diagram shows the intersection between LMRGs and DEGs. (B) The heatmap depicts the relative expression of 23 genes involved in lactate metabolic activities that were shown to vary between tumor and normal specimens. (C-E) In the TCGA, GEO, and FJZL cohorts, low-risk group patients had a favorable OS rate as opposed to those in the high-risk group. (F) The ROC curve for HNSCC patients’ OS over 1, 3, and 5 years in the TCGA, GEO, and FJZL cohorts. (G) Construction of a nomogram based on riskscore and clinical characteristics. (H) Calibration curves of the nomogram for predicting OS in TCGA-HNSCC patients. (I) Correlation analysis of clinical features and risk score

Evaluation of the prognostic signature’s relationship to clinical parameters

Additionally, we analyzed the risk signature’s predictive value in TCGA HSNSCC patients with a variety of clinical and pathological features (Table S6). As illustrated in Fig. 2I, risk scores increase with the T3-4 stage, HPV- status, perineural invasion (PNI) positivity, and SD/PD (Fig. 2I). Furthermore, the risk score demonstrated significant prognostic differentiation performance among patients with TNM stage, HPV status, PNI status, lymphovascular invasion (LVI) status, and clinical outcomes (Supplementary Figs. 1G-L).

Gene set enrichment analysis

Several immune-associated pathways were enriched in low-risk phenotypes, including T cell activation modulation, antigen processing, and presentation, positive regulation of immune effector mechanisms, and the JAK-STAT signaling pathway (Supplementary Figs. 2 A-B). Meanwhile, we found that high-risk phenotypes were enriched in glycolysis, glycogen synthesis process, pyrimidine metabolism, and MYC pathway (Supplementary Figs. 2 C), suggesting that high-risk score was involved in cancer metabolism and oncogenic biological processes. The CCP scores supported this finding, as patients in the high-risk group exhibited elevated CCP scores (Fig. 3A). The high-risk group exhibited significantly elevated activities in hypoxia, glycolysis, and lactic acid pathways, suggesting intense oxidative and aerobic glycolysis in proliferating cancer cells (Fig. 3B).

Fig. 3

The correlation of risk signature and immune microenvironment of HNSCC. (A) Patients in the high-risk group have high CCP scores. (B) Difference between low- and high-risk groups at six energy metabolism pathways. (C) In HNSCC, 53 immunoregulators are differentially expressed (MHC, receptors, and chemokines). (D) Difference between low- and high-risk groups at distinct stages of the cancer-immune cycle. (E) Relationship of the risk score with infiltration levels of severe TIICs, as determined by seven separate algorithms

The relationship between the risk model and inflamed TME

In our results, low-risk HNSCC patients exhibited better immune and estimate scores compared to high-risk HNSCC patients but had lower tumor purity (Supplementary Figs. 3 A). The risk signature showed a significant negative correlation with numerous immunomodulators (Fig. 3C). The activities associated with the cancer-immune cycle involve both direct and systemic actions by chemokines and other immunoregulators [40]. Most steps of the cycle, such as priming and activation (Step 3) and immune cell transportation to tumors (Step 4), were enhanced in the low-risk group (Fig. 3D).

The risk signature was inversely correlated with DC cells, CD8 + T cells, and B cells across seven distinct algorithms (Fig. 3E). Similar findings were observed in 20 tumor-infiltrating immune cells (TIICs) analyzed by ssGSEA algorithms (Supplementary Figs. 3B). We identified that gene biomarkers of these immune cells were upregulated in the low-risk group (Fig. 4A, Supplementary Figs. 3 C-E). Studies indicate that tumor-infiltrating B cells and tumor-associated tertiary lymphoid structures (TLS) are pivotal in enhancing immunotherapy responsiveness [41]. TLS primarily develop in inflamed tissues, and our findings indicate an inverse correlation between risk score and TLS score (Fig. 4B). Inhibitory immune checkpoints, such as PD-1 and PD-L1, showed elevated expression levels in the inflamed TME [42]. The risk signature was inversely correlated with the majority of ICIs (Fig. 4C).

As previously stated, solid tumors can be categorized into two types: hot and cold tumors, with hot tumors being more sensitive to immunotherapeutic regimens [43]. We employed an unsupervised clustering method to categorize 501 HNSCC specimens from the TCGA cohort into hot and cold tumor samples, based on hot tumor signature genes (Fig. 4D, Supplementary Figs. 3 F). In hot tumors, the low-risk group showed significantly higher expression compared to cold tumors (Fig. 4E), indicating its involvement in the clinical responsiveness to immunotherapy. The risk signature showed a strong correlation with CYT and GEP in HNSCC, with both being elevated in the low-risk group (Fig. 4F).

Furthermore, we analyzed TCR and BCR repertoires from the TCGA HNSCC cohort. Shannon entropy measurements indicated that TCR and BCR diversity values were highest in the low-risk groups (Fig. 4G). According to our findings, CYT, the percentage of tumor-infiltrating lymphocytes (TILs), T-cell co-stimulation levels, inflammation-promoting factors, HLA, and checkpoint levels decreased in the TCGA cohort as the risk increased (Fig. 4H). Another potential predictor of TME is carcinogenic dedifferentiation, characterized by cancer stem cell (CSC) features. The mRNAsi and EREG-mRNAsi are quantitative measure that assesses the “stemness” characteristics of tumor cells by analyzing their gene expression profiles. EREG-mRNAsi not only reflects the stemness of tumor cells but also highlights the role of growth factor signaling in maintaining these characteristics. High mRNAsi and EREG-mRNAsi values are often associated with increased tumor aggressiveness, metastasis, and reduced immune infiltration. Both mRNAsi and EREG-mRNAsi indicators were significantly elevated in the high-risk group (Fig. 4I). We found significant enrichment in inflammation and modulation of inflammation pathways in the low-risk score group (Supplementary Figs. 3G). In summary, the low-risk score group may have formed an inflammatory TME.

Fig. 4

Association of the risk score with immunological subtypes and indicators associated with immunotherapy. (A) Association of risk score with CD8 + cell effector genes. (B) Association of risk score with TLS. (C) Relationship of risk score with 10 inhibitory immune checkpoints, including TIGIT, IDO1, PD-1, LAG-3, CTLA-4, PD-L1, and TIM-3. (D) Heatmap plot demonstrated enrichment of hot tumor signature genes (CXCR4, CXCL10, CXCL9, CD4, CXCR3, CD3E, CD8B, CXCL11, PDCD1, CD274, CD8A, and CCL5) in hot tumor specimens. (E) Risk score was considerably reduced in hot tumors, demonstrating its involvement in clinical responsiveness to immunotherapeutic regimens. (F) Differences in GEP and CYT between high- and high-risk groups. (G) Differences in TCR and BCR diversity values between high- and low-risk groups. (H) Differences in the enrichment scores of most immunotherapeutic-positive gene signatures between high- and low-risk groups. (I) Differences in the tumor stemness index between high- and low-risk groups

The investigation revealed that the low-risk group showed significantly higher responsiveness to immunotherapy compared to the high-risk group (Fig. 5A-C). Patients with low-risk scores demonstrated higher responsiveness to PD-1 and CTLA4 antagonists based on IPS and three melanoma cohorts, suggesting they would benefit from immunotherapy (Fig. 5D-G).

Fig. 5

The response of immunotherapy of low- and high-risk groups. (A) Differences in TIDE score between high- and elevated-risk groups in the TCGA cohort. (B) Differences in TIDE score between high- and low-risk groups in the GEO cohort. (C) The anticipated immunotherapy (TRUE/FALSE) response rate to anti-PD-L1 in high- and low-risk groups in the TCGA cohort. (D) Differences in IPS between high- and low-risk groups. (E-G) Differences riskscore in melanoma- GSE91061, melanoma-PRJEB23709, and melanoma- GSE100797

Single-cell immune landscapes and cellular communication

19 distinct cell types were identified in the GSE103322 cohort using two-dimensional spatial visualization through sTNE analysis (Fig. 6A). Major cell subsets were identified as malignant cells, endothelial cells, immune cells, and fibroblasts (Fig. 6B). The lactate metabolism risk score was highest in malignant cells, followed by fibroblasts (Fig. 6C). Malignant cells exhibited significantly elevated activities in oxidative phosphorylation, glycolysis, and lactic acid pathways, suggesting a hypermetabolic state (Fig. 6D). To investigate lactate metabolism differences in infiltrating immune cells within the TME, 5905 HNSCC cells from 21 specimens were clustered into 10 cell types: Myofibroblast, Malignant, Treg, CD8 T cell, Fibroblast, Endothelial, Macrophage, B cell, Mast, Dendritic, and Myocyte (Fig. 6E). The detailed annotation results are shown in Supplementary Figs. 4 A. Violin plots depict the riskscore of each sample (Supplementary Figs. 4B). Samples were categorized into high and low groups based on median values. The low-risk group contained higher levels of infiltrating Tregs, higher CD8 + T cell ratios, and lower levels of malignant cells (Fig. 6F, Supplementary Figs. 4 C). Subsequently, we conducted functional exploration. The primary pathways enriched for differential genes between high- and low-risk groups involved intercellular adhesion and immune cell activation (Fig. 6G). Active pathways differed between high- and low-risk groups, with LIFR, GRN, CCL, GAS, LIGHT, and SPP1 pathways active in the high-risk group, and IL-1, TGFb, and IL10 pathways active in the low-risk group (Supplementary Figs. 4D, Fig. 6H). Figure 5F illustrates the interaction between malignant cells and other cells via ligand-receptor binding. Intercellular communication involving WNT and SPP1 was upregulated in the high-risk group compared to the low-risk group. The SPP1 signaling pathway was enriched in malignant cells in the high-risk group, whereas SPP1 was enriched in macrophages in the low-risk group (Fig. 6I). TF dysregulation significantly influences tumor progression. We subsequently identified upstream transcription factors involved in lactate metabolism. High risk group exhibited activated TFs (including FLI1, IRF4, RUNX3, FOXP3 extended, ZNF467 extended, IKZF1) and inhibited TFs (including MYC extended, HTATIP2, SREBF2, NFE2L2; Fig. 6J-K, Supplementary Figs. 4E-F).

Fig. 6

Immune landscapes and cellular communication at the single-cell level. (A) t-distributed stochastic neighbor embedding (t-SNE) of 19 clusters in GSE103322 cohort. (B) t-SNE of major cell subsets in GSE103322 cohort. (C) The lactate metabolism riskscore in major cell subsets by violin plot. (D) Differences in six energy metabolism pathways among major cell subsets. (E) t-SNE of 10 cell types in GSE103322 cohort. (F) The proportion of cell composition in high- and low-risk groups. (G) The major pathways enriched for differential genes between the high- and low-risk groups. (H) The up-regulated and down-regulated signaling between malignant cells and other cells. (I) SPP1 signaling pathways in high- and low-risk groups. (J) Differential TFs between high- and high-risk groups. (K) Top 10 TFs target genes between high- and high-risk groups

PYGL expression, survival, and metabolism analyses

For better clinical application, PYGL, which was the most important prognostic gene in 16 LMRG risk signatures, was screened by the SWSFS algorithm and Ranger algorithm (Fig. 7A). Correlation analysis indicates that PYGL is significantly positively correlated with the risk score, showing stronger association compared to other lactate-related genes (Supplementary Figs. 5 A-B). Furthermore, PYGL has a close relationship with LDHA, a key molecule in lactate metabolism, highlighting the importance of PYGL in this metabolic pathway (Supplementary Figs. 5 C). Therefore, PYGL can, to some extent, serve as a representative marker for the risk score. In the TCGA cohort, PYGL expression was considerably reduced in paired paracancerous tissues as opposed to tumor samples (Fig. 7B). Additionally, through the TIMER website, we discovered that PYGL is expressed at high levels in most cancers (Supplementary Figs. 6 A). Moreover, CPTAC dataset and our IHC cohort both revealed that tumor tissue PYGL expression level was greater than that of normal tissue (Fig. 7C-E). In the TCGA and GEO cohort, individuals with HNSCC who had lower PYGL expression levels exhibited substantially lower OS times as opposed to those who had elevated PYGL expression levels (Fig. 7F). Subsequently, the function of PYGL in HNSCC was explored. Figure 7G demonstrates a significant association between PYGL and the clinical T stage (T1-2 vs. T3-4) in the TCGA sample.

Next, we study the relationship between PYGL and lactate metabolism. The expression level of PYGL is positively associated with global lactylation in our IHC cohort (Fig. 7H). We discovered that hypoxia, glycolysis, and lactic acid pathways positively correlated with PYGL, indicating high PYGL was significantly involved in energy metabolism in tumor progression (Fig. 7I). To further verify this result, we performed PYGL knockdown in SCC7 cells (Fig. 7J). Lactate concentration in the medium supernatant decreased following PYGL knockdown (Fig. 7K). Moreover, PYGL and LDHA expression was remarkably increased at mRNA (Fig. 7L-M) under hypoxic conditions. PYGL and HIF1A expression was significantly increased at protein levels under hypoxic conditions (Fig. 7N). Collectively, the data suggest that hypoxia induces PYGL expression and lactate accumulation.

Fig. 7

PYGL expression and survival Analyses. (A) By Ranger and SWSFS algorithm, PYGL was identified as the most important hub gene in the risk signature. (B) PYGL expression in paired 44 paracancerous and tumor tissues from the TCGA cohort. (C) PYGL protein expression in tumor and normal samples in the CPTAC cohort. (D-E) IHC staining images and PYGL protein levels in adjoining non-tumor tissue and OC tissue. (F) Kaplan-Meier curves of OS for PYGL in the TCGA and GEO cohorts. (G) The association of the PYGL expression with the clinical T stage. (H) Analyses of correlation between PYGL and lactylation level. (I) Analyses of correlation between PYGL and six energy metabolism score. (J) Wb was performed to detect the efficiency of PYGL-shRNA transfection. (K) Lactate concentration was detected in medium supernatant after transfection with PYGL shRNA for 48 h. (L-M) PYGL and LDHA mRNA expression in SCC7 cells exposed to hypoxia for 0–24 h examined by qRT-PCR. (N) HIF1A and PYGL protein expression in SCC7 cells 0–24 h examined by western blotting

Effect of PYGL on infiltrating immune cells of TME

The TIICs of HNSCC specimens were utilized to additionally study the interaction between PYGL expression and the TME. We obtained PYGL expression data for single-cell subtypes from TISCH. PYGL was expressed by macrophages or malignant cells in BRCA, Glioma, and HNSC (Fig. 8A, Supplementary Figs. 6B-C). In detail, the scatter plot displays the distribution of PYGL within each celltype in HNSC single-cell sequencing database (GSE103322, GSE139324; Fig. 8B-C). We identified PYGL expression in M2 or undefined macrophages in GSE135324. Pseudotime trajectory analysis revealed the temporal sequence of cell differentiation, positioning M1 and M2 macrophages predominantly at the terminal stage of the differentiation pathway (Fig. 8D). PYGL was predominantly located in the early macrophages and impeded M1 polarization (Fig. 8E). Spatial transcriptional data from SpatialDB were utilized to illustrate the spatial overlap of PYGL and macrophage biomarker CD68 in BRCA, PRCA cancer, and melanoma tissues (Supplementary Figs. 7 A). PYGL and CD68 markers exhibited similar spatial distributions, suggesting potential co-expression. We quantified the gene expression of various macrophage markers in M0, M1, and M2 macrophages using THP-1 and RAW267.4 cells to further characterize PYGL and polarized macrophages (Fig. 8F-G, Supplementary Figs. 7B). RT-PCR analysis revealed a significant decrease in PYGL expression in M1 macrophages compared to M0 and M2 macrophages. Additionally, we observed that CD86 protein expression was significantly increased in PYGL-knockdown THP-1 M0 and RAW267.4 M0 cells (Fig. 8H). Furthermore, RNA-seq analysis of PYGL-knockdown THP-1 M0 macrophages revealed an upregulation of HLA-DPA1, HLA-DQB1, HLA-DRA, and HLA-DRB5 in the knockdown group. These findings suggest that PYGL may influence M1 macrophage polarization through MHC class II-dependent mechanisms. (Supplementary Figs. 7 C).

Fig. 8

Associations between PYGL and macrophages. (A) PYGL expression in HNSCC single-cell clusters obtained from TISCH online tool. (B) The expression of PYGL in tSNE plot of GSE103322. (C) The distribution of immune cell clusters and PYGL in UMAP plot of GSE139324. (D) Pseudotime trajectory analysis of macrophages types based on celltype (left), pseudotime (right). (E) Pseudotime trajectory analysis of PYGL in different macrophages types. (F-G) CD68, TGF-β, CCL22, PYGL, CXCL10, and CXCL9 mRNA expression in THP-1 cells by qRT-PCR. (H) CD206, PYGL, and CD86 protein expression in THP-1 M0 and RAW267.4 cells by western blotting

In addition to macrophages, we investigated other contributors to immunosuppression. The findings confirmed a strong association between CD8 + T cells, B cells, and PYGL (Fig. 9A-B). Moreover, immunohistochemistry (IHC) testing indicated that an elevation in the PYGL protein expression levels contributed to a decrease in the CD8 protein level (indicating tumor-infiltrating CD8 + T cells), confirming the results at the transcriptional level (Fig. 9C). Furthermore, multicolor immunofluorescence analysis revealed that CD8 and PD-L1 were co-localized in HNSCC tissues. Following that, we discovered that PYGL was expressed at the lowest levels within the inflamed phenotype, which was negatively linked to CD8 + T-cell infiltrates (Fig. 9D-E). Given the role of chemokines and their receptors in immune cell mobility within the TME, particularly the recruitment of CD8 + T cells into tumors, we found a strong inverse correlation between PYGL expression and the expression of CXCR3, CXCR9, CCL4, CCR5, and CCR2 in HNSCC (Fig. 9F-G).

The clinical responsiveness of ICBs was subsequently investigated. IHC staining demonstrated a reduction in the protein level of PD-L1 with the increase of PYGL protein expression (Fig. 9H). An inverse correlation between PYGL and the expression of LAG3, PDCD1 (PD-1), CD274 (PD-L1), and CTLA-4 was validated in the TCGA cohort (Fig. 9I, Supplementary Figs. 8B). Furthermore, in the TCGA cohort, PYGL was shown to be positively linked to the enrichment scores of the majority of immunotherapeutic-positive gene profiles (Supplementary Figs. 8 C). Elevated TIDE scores were shown to be remarkably linked to higher PYGL expression levels (Fig. 9J). The final mulberry plot suggested that PYGL could represent the lactate metabolism risk model (Fig. 9K), monitoring clinical efficacy and Screening patients who benefit from immunotherapy.

Fig. 9

Associations between PYGL, the immune infiltration, and the responsiveness to immunotherapy in HNSCC. (A) Relationship between PYGL and infiltration levels of CD8 + T cell in the TCGA cohort. (B) Association of PYGL with infiltration levels of B cell in the TCGA cohort. (C) Inverse link between PYGL and CD8 in OC tissue microarray (n = 48) by IHC. (D) The HNSCC samples were classified into two groups depending on the existence or absence of CD8 + T cells into inflamed and non-inflamed phenotypes. Immunofluorescence showed lower PYGL positivity and higher CD8 and PD-L1 positivity in inflammatory TME. (E) Mean densitometric analysis revealed that PYGL was inversely linked to CD8 + T cell infiltrates. (F) Relationship between PYGL and chemokines. (G) Relationship between PYGL and chemokine receptors. (H) Inverse correlation between PYGL and CD8 in OC tissue microarray (n = 48) by IHC. (I) Correlation analysis of PYGL with CD274 (PD-L1) in the TCGA cohort. (J) Analyses of correlation between PYGL and TIDE score. (K) Mulberry plot showing the relationship between risk signature, PYGL, immune subtype, and predicted immune efficacy

Identifying potential drugs targeting PYGL

We investigated the expression of PYGL in NCI-60 cell lines to identify genes predictive of drug sensitivity and analyzed its correlation with drug response. The findings indicated a significant correlation between PYGL and the IC50 of eight drugs (Fig. 10A), guiding the treatment of patients with high PYGL expression. AutoDock 4.2 was utilized for virtual screening to identify potential PYGL-targeting drugs. Strong interactions were observed between PYGL and raloxifene (docking score = -5.83), elesclomol (docking score = -8.223), and irofulven (docking score = -5.497) (Fig. 10B-D). Elesclomol, an anticancer drug, targets mitochondrial metabolism and induces copper-dependent cell death [44]. We found that elesclomol combined with copper ions caused increased cell death in PYGL-knockdown cells. Furthermore, PYGL exhibited a negative correlation with LIPT1, a cuproptosis regulatory protein, as validated at the protein level. Hence, elesclomol shows potential as a treatment for HNSCC.

Fig. 10

Correlation of PYGL and drug response. (A) Scatter plots of the top six kinds of associations between PYGL and chemotherapy drug sensitivity. (B-D) AutoDock-derived structure of raloxifene (B), elesclomol (C), and irofulven (D) bound to PYGL. (E) Cell viability is measured in SCC7 PYGL- knockdown cells treated with elesclomol combined with copper ions (40 nm, 60 nm, 80 nm, 100 nm, and 120 nm) or vehicle. (F) Correlation analysis of PYGL with LIPT1 in the TCGA cohort. (G) LIPT1 protein expression in SCC7 PYGL- knockdown cells examined by western blotting

Discussion

This study used RNA-seq data from TCGA and 16 differentially expressed LMRGs to develop a prognostic signature for risk classification and clinical outcome prediction. A highly significant link between the risk model and clinicopathological characteristics is discovered. Immune-associated pathways were shown to undergo substantial enrichment in the low-risk group. The high-risk group was shown to be inversely linked to many immunoregulators and to have a TME that was not inflamed. Patients with HNSCC who have low-risk scores might gain from ICB treatment. Machine learning analysis demonstrated that PYGL was a hub prognostic gene related to HNSCC patients, and played a significant function in TME.

Lactate gathers in tumors as a result of metabolic abnormalities in cancerous cells [45]. Ping et al. discovered that intracellular lactate levels in TIL from gastric carcinomas were dramatically elevated [46]. Elevated lactate levels were shown to be inversely linked to the count of TH1 cells and CTLs in the tumor, indicating that the immunological potential inside the TME has been changed and impaired [6]. Lactate is a dominant metabolite in the TME that suppresses the activation and proliferation of dendritic cells, NK cells, and CD8 + T cells, thereby affecting the metabolic activity of adaptive and innate immune responses [47]. It is commonly accepted that an inflamed TME is linked to a better patient prognosis and a higher sensitivity to ICB [48,49,50,51]. This investigation found significantly higher CD8 + T cell infiltration levels in the low-risk group compared to the high-risk group. Chemokines regulate the function of CD8 + T lymphocytes in tumors. Chemokines like CCR5, CXCL9, and CXCL10 are associated with the extent of CD8 + T cell infiltration in melanoma [52]. STAT proteins are important for immune cell development and function since they are the primary signaling proteins for inflammatory markers [53]. The low-risk group showed significant enrichment in genes regulating immune pathways, such as the JAK-STAT pathway and the inflammatory response, alongside higher CYT and GEP scores—both strong indicators of antitumor immune activity. Brand et al. proposed a mechanism for lactate-induced immune inhibition, demonstrating that tumor acidosis and lactic acid suppress the nuclear factor of activated T cells (NFAT) [54]. This suppression, affecting tumor-infiltrating NK cells and CD8 + T lymphocytes, leads to a decrease in IFNγ production [11]. The primary immune evasion mechanisms in the T-cell-inflamed tumor microenvironment include: (i) elevated PD-L1 expression leading to T cell suppression through PD-L1/PD-1 interaction; (ii) upregulation of IDO; (iii) CCL22-driven recruitment of regulatory T cells (Tregs); and (iv) selection of tumor cells with reduced immunogenicity [49]. Our findings suggest that the low-risk group exhibited higher levels of ICIs and CCL22 compared to the high-risk group. ICB immunotherapy seems to be most effective in those with low-risk scores.

Our study underscores the pivotal role of lactate metabolism in shaping the tumor microenvironment (TME). Single-cell analysis revealed distinct metabolically active cell populations that modulate immune responses and drive tumor progression. Notably, the correlation between reduced lactate metabolism, enhanced CD8⁺ T cell infiltration, and improved prognosis highlights the potential of targeting metabolic pathways to optimize immunotherapy efficacy.

PYGL is a gene that encodes the glycogen phosphorylase distributed in cells and is associated with cell metabolic processes. PYGL has been suggested as a hypoxia signal in breast cancer that has a prognostic value [55]. A meta-analysis of over 2,000 cases revealed elevated PYGL levels in multiple cancer types, including clear cell renal cell carcinoma, seminoma, and brain cancer, particularly within the hypoxic TME [56]. In the TME, hypoxia is a prominent feature, and its impacts on immune cells have been postulated to be essential contributors to the process of tumor immune evasion [57]. Our findings suggest that PYGL knockdown promotes M1 macrophage polarization and is strongly associated with PD-L1 expression in tumor cells and infiltrating CD8⁺ T cells, indicating its potential role in modulating the immune microenvironment and influencing immunotherapy response in HNSCC. Additionally, we noticed that PYGL was involved in the copper-dependent cell death process. PYGL competitively binds to elesclomol, preventing the chelation of elesclomol with copper ions to form a complex. Knockdown of PYGL enhances copper ion accumulation, leading to the upregulation of LIPT1, a key enzyme involved in protein lipoylation. Ultimately, this process triggers cuproptosis.

Although we developed a predictive signature and provided further insights for improved HNSCC treatment, our work has several shortcomings. First, the size of the sample may be inadequate. The results must be confirmed by other independent cohorts to demonstrate the therapeutic significance of the risk model. Second, given the significant differences in the origin and molecular properties of the different HNSCC subtypes, a more nuanced approach based on the specific origin and potential molecular characteristics of each subtype is essential. In the future, we will also further analyze each subtype to find better biomarkers to provide precise personalized treatment for patients. In addition, more functional studies are required to investigate the fundamental processes of PYGL in TME. The effects of PYGL on the proliferation and metastasis potential of HNSCC and the related pathways involved need further exploration.

Conclusion

In this study, we finally constructed a lactate metabolism-related risk model to accurately predict the prognosis of patients with HNSCC as well as the efficacy of multiple immunotherapies. In addition, we described the lactate metabolism of cells in TME at the single-cell level. Finally, PYGL, a lactate metabolism-related molecule in the model was shown to promote the malignant progression of HNSCC cells and also play an essential role in TME.