The workflow and standardization of the dataset in our study are depicted in Supplementary Figure S1. We compared the expression levels of 289 LMRGs across the TCGA-HNSC cohort, with detailed results summarized in Supplementary Table S1. Employing a correlation algorithm, we analyzed the LMs for each HNSC patient. Based on the median LMs value, patients were stratified into two distinct groups: LMs-high and LMs-low. We examined the correlation between the LMs status and clinical variables (such as gender, age, M stage, N stage, and T stage). Our analysis revealed no significant associations between LMs grouping and gender, age, or T stage. However, intriguing correlations were observed with the M stage and N stage, as depicted in Figure 1A. Notably, LMs-high patients were predominantly observed in the M1 stage and exhibited a higher proportion in the N2 stage compared to LMs-low patients. Furthermore, we generated Kaplan-Meier curves for the LMs-high and LMs-low groups (Figure 1B). These curves revealed a significant difference in overall survival, with patients in the LMs-low group exhibiting a more favorable prognosis compared to those in the LMs-high group (p = 0.015).

Moreover, we conducted an investigation to identify genes that were differentially expressed (DEGs) between the LMs-low and LMs-high groups. The results of this analysis are vividly depicted in the volcano plot (Figure 1C), which highlights 2334 genes that were significantly upregulated and 1810 genes that were downregulated in the LMs-high group compared to the LMs-low group. To perform Weighted Gene Co-expression Network Analysis (WGCNA), we identified 20 gene modules, among which the MEblue module exhibited the highest correlation with LMs (Figure 1D). We then carried out an intersection analysis between DEGs from the TCGA-HNSC dataset and LMRGs, resulting in 96 LMRDEGs (refer to Figure 1E). Subsequently, we performed another intersection analysis between LMRDEGs and MEblue genes, leading to the identification of 45 key genes (Figure 1F; Supplementary Table S2). To visually represent the differential expression of these key genes between the LMs-low and LMs-high groups, we generated a heatmap (Figure 1G).

Furthermore, we conducted a comparative analysis of the expression levels of critical genes across three distinct datasets: TCGA-HNSC, GSE6631, and GSE107591. In the TCGA-HNSC dataset, we observed a substantial increase in the expression of 45 key genes in the LMs-high group compared to the LMs-low group. Similarly, within the GSE6631 dataset, the LMs-high group exhibited significantly higher expression of 11 genes (C1QBP, COX4I1, COX5A, COX6B1, LDHB, NDUFA2, NDUFB3, NDUFS3, NDUFS4, SOD1, UQCRB) when compared to the LMs-low group. Likewise, in the GSE107591 dataset, 12 genes (C1QBP, COX4AI1, CYC1, MRPS34, NDUFA12, NDUFA6, NDUFA8, NDUF4F2, NDUFB9, NDUFV1, SOD1, TXN2) displayed significantly increased expression in the LMs-high group relative to the LMs-low group (Supplementary Figure S2).

The initial investigation focused on evaluating the prognostic significance of 45 crucial genes in HNSC patients. This assessment was achieved through the utilization of univariate Cox regression (Figure 2A) and Kaplan-Meier analysis to establish the association between gene expression and molecular subtypes. Employing a screening threshold of P < 0.05, we identified COX8A (Figure 2B), LDHB (Figure 2C), NDUFA1 (Figure 2D), NDUFA4 (Figure 2E), and NDUFB3 (Figure 2F) as key genes associated with prognosis. Subsequently, we employed unsupervised consensus clustering to explore the underlying structure of the HNSC patient cohort. This analysis revealed strong intragroup correlations and low intergroup correlations when k = 2. This finding suggests that HNSC patients can be categorized into two distinct subtypes: cluster 1 (n = 342) and cluster 2 (n = 160) (Figures 2G–I). Notably, principal component analysis (PCA) exhibited significant disparities in lactic acid metabolism-related transcriptome profiles between cluster 1 and cluster 2 (Figure 2J). Furthermore, we examined the distribution of the five key genes (COX8A, LDHB, NDUFA1, NDUFA4, NDUFB3) across the different immune feature subtypes identified by consensus clustering. Significant variations in the expression patterns of these genes were observed between cluster 1 and cluster 2 (Figure 2K).

We analyzed gene mutations in five key genes (COX8A, LDHB, NDUFA1, NDUFA4, NDUFB3) using samples from the TCGA-HNSC dataset and the cBioPortal database. The findings (Figure 3A) showed that NDUFB3 had the highest total mutation count, accounting for 2.4% of the total samples in the TCGA-HNSC dataset, mainly presenting as significant amplification. We also observed a significantly high correlation (p < 0.01) for all pairwise combinations of the 5 key genes (Figure 3B).

To further investigate their prognostic significance in HNSC, these five key genes were subjected to multivariate Cox regression analysis (Figure 3C; Supplementary Table S3). The results demonstrated that patients with elevated expression levels of COX8A, LDHB, NDUFA4, and NDUFB3 had a poorer prognosis. Furthermore, we conducted prognostic calibration analysis at discrete time points: one-year (Figure 3D), three-year (Figure 3E), and five-year (Figure 3F). These analyses utilized the nomogram derived from the multivariate Cox regression model to generate corresponding calibration curves (Figures 3D–F). Notably, our DCA results indicated that the clinical predictive effect of the constructed Cox regression prognostic model was most pronounced at the five-year mark, followed by the three-year and one-year intervals (Figures 3G–I).

In order to assess the correlation between the risk score of the Cox model and the LMs, scatter plots were generated, revealing a moderate correlation (r = 0.684) (Figure 3J). Moreover, the KM curve (Figure 3K) demonstrated a worse prognosis for HNSC patients with high-risk scores compared to those with low-risk scores. Additionally, the ROC curve (Figure 3L) showed that the Cox risk score had a certain predictive ability for the overall survival of HNSC patients, with improved prediction as time increased. Finally, we visualized the distribution of lactate metabolism-related risk scores among HNSC patients by generating a heatmap based on the Cox risk scores (Figure 3M).

The TCGA-HNSC dataset samples were categorized into high-risk and low-risk groups based on the median Cox risk score. To identify DEGs associated with the Cox risk score, we conducted a differential expression analysis on the expression profile data of the TCGA-HNSC dataset. The DEGs for our study were selected based on the criteria of p.adjust <0.05 and the top 20 genes with positive and negative log fold change (logFC). These genes included GHRH, GSTA1, GAGE12J, UCN3, WIF1, CPLX2, NTS, DDC, CES1, SHCBP1L, MT3, FDCSP, HMX1, USH1C, RPRML, GAL, CLPSL1, OR11H4, DEFB126, GSTA3, NEU2, FLG2, TGIF2LX, MAGEB1, KPRP, AL354761.1, LCE2B, PCDHGC5, KRT76, KRTAP9-8, RXFP3, GCOM1, DSG1, LCE6A, LCE1A, SMR3B, PLA2G4D, LCE2C, INSYN2B, and FRG2B. In order to visualize the results of the differential analysis, we generated a volcano plot (Figure 4A) and a heat map (Figure 4B).

To gain insights into the biological processes, molecular functions, cellular components, and biological pathways associated with these 40 DEGs, we performed GO and KEGG enrichment analyses. The GO functional enrichment analysis (Figure 4C) and KEGG enrichment analysis revealed that the 40 DEGs were significantly enriched in various biological processes, including epidermis development, skin development, epidermal cell differentiation, keratinization, keratinocyte differentiation, lipid catabolic process. In terms of cellular components, they were associated with presynapse, axon terminus, neuron projection terminus, distal axon, cornified envelope, and terminal bouton. The molecular functions identified neuropeptide receptor binding, hormone activity, G protein-coupled receptor binding, and neuropeptide hormone activity. Additionally, the pathways identified were Neuroactive ligand-receptor interaction and Drug metabolism - other enzymes. (Supplementary Table S4).

To evaluate the impact of gene expression levels on the categorization of high and low risk score groups as determined by the Cox regression model, we conducted Gene Set Enrichment Analysis (GSEA). The GSEA results revealed significant enrichment of DEGs in various pathways (Figure 4D). Notably, these pathways included Gluconeogenesis (Figure 4E), Pyruvate metabolism (Figure 4F), Cellular Response to Hypoxia (Figure 4G), and PI3K/AKT Signaling in Cancer (Figure 4H), among others, within the TCGA-HNSC dataset (Supplementary Table S5).

Utilizing the STRING database, we elucidated the interactions among differentially expressed genes (Supplementary Figure S3A). Subsequently, we identified central genes by intersecting the highest-ranked genes obtained from five disparate algorithms (Degree, MNC, MCC, EPC, and Closeness), as shown in Supplementary Figure S3B. Moreover, the GeneMANIA algorithm was employed to visualize the interaction network of genes functionally associated (Supplementary Figure S3D). In addition, we predicted the miRNA interactions for the eight central genes, as depicted in Supplementary Figure S3C, which led to the construction of an mRNA-miRNA interaction network. This network consists of six central genes (NTS, GHRH, LCE2B, DDC, PCDHGC5, and MT3) and 53 miRNA molecules, encompassing a total of 54 mRNA-miRNA interaction pairs. (Supplementary Figure S3C; Supplementary Table S6).

In order to evaluate immune infiltration and scoring in the TCGA-HNSC dataset, a range of immune and stromal scores were examined based on the Cox risk score. The findings indicated significant variations (P < 0.05) among different groups within this dataset for the ESTIMATE score (Figure 5A), Immune score (Figure 5B), and Stromal score (Figure 5C).

Subsequently, an analysis was conducted to compare immune cell infiltration between the low risk and high risk score groups as identified by Cox regression analysis (Figure 5D). Additionally, a correlation heat map was generated to assess the relationship between immune cell populations (Figure 5E). Furthermore, an evaluation of the correlation between immune cells and hub genes (GAL, NTS, GHRH, UCN3, LCE2B, DDC, PCDHGC5, MT3) was performed (Figure 5F). The results demonstrated significant differences (p < 0.05) in infiltration for 21 immune cell types (Central memory CD4 T cell, Central memory CD8 T cell, Effector memory CD4 T cell, Effector memory CD8 T cell, Gamma delta T cell, Immature B cell, Memory B cell, Regulatory T cell, T follicular helper cell, Type 1 T helper cell, Type 17 T helper cell, Type 2 T helper cell, Activated dendritic cell, Eosinophil, Mast cell, MDSC, Monocyte, Natural killer cell, Natural killer T cell, Neutrophil, Plasmacytoid dendritic cell) between the low and high risk score groups.

The somatic mutation profiles of patients in the Cox regression analysis-based low-risk and high-risk score groups were examined in the TCGA-HNSC dataset. The findings revealed that the primary types of somatic mutations observed in both groups were Missense Mutation and Nonsense Mutation, with a higher prevalence of missense mutations. Additionally, single-nucleotide polymorphisms (SNPs) were identified as the most common type of mutation observed in HNSC patients within the Cox low-risk and high-risk score groups. In the low-risk score group, the most frequent single nucleotide variant (SNV) among HNSC patients was C > T, followed by C > A (Figure 6A). Conversely, in the high-risk score group, the most prevalent SNV among HNSC patients was C > T, followed by C > G (Figure 6B).

The Copy Number Variation (CNV) analysis demonstrated an increased number of CNVs at 11q13.3 (Figure 6C) and a decreased copy number at 9p21.3 in the mutant group of the TCGA-HNSC dataset (Figure 6D). Among the 8 hub genes analyzed for CNV, GHRH exhibited the highest amplification frequency, while PCDHGC5 had the highest deletion frequency (Figure 6E).

Furthermore, an evaluation of microsatellite instability (MSI) and tumor mutation burden (TMB) was conducted among the low-risk and high-risk score groups of HNSC patients based on Cox regression analysis in the TCGA-HNSC dataset. The sensitivity to immune therapy across different groups of HNSC patients was assessed using the TIDE algorithm. The analysis indicated no statistically significant differences in TMB between the various groups (Figure 6G). However, a statistically significant difference (P < 0.05) in MSI (Figure 6F) and TIDE scores (Figure 6H) was observed among the different groups of HNSC patients in the dataset.

To delve deeper into the clinical relevance of the risk score linked to lactic acid metabolism, a comparative analysis was conducted to examine the distribution of gender and disease stage between the high-risk and low-risk groups (Figures 7A, B). Subsequently, a prognostic model was constructed for patients with HNSC based on the risk score obtained from Cox regression analysis, which considered mitochondrial energy metabolism as well as clinical and pathological features including age and disease staging. The resulting model was visualized using a column line chart (nomogram) (Figure 7C). The accuracy of the model was evaluated using calibration curves, which demonstrated a strong agreement between the estimated one-year (Figure 7D), three-year (Figure 7E), and five-year (Figure 7F) overall survival (OS) values and the actual observations of patients. Additionally, the ROC curve (Figure 7G) showed that the mode had a certain predictive ability for the overall survival of HNSC patients.

The expression levels of eight hub genes were examined in different types of head and neck cancers, including tongue cancer, hypopharyngeal cancer, laryngeal cancer, and oropharyngeal cancer. The findings revealed that, in comparison to their levels in normal tissues, the genes GAL, GHRH, LCE2B, PCDHGC5, and MT3 were notably upregulated in tongue cancer, while NTS was downregulated. For hypopharyngeal cancer, GHRH, UCN3, LCE2B, DDC, and PCDHGC5 were upregulated. All eight key genes displayed upregulation in laryngeal cancer. Regarding oropharyngeal cancer, NTS, LCE2B, DDC, PCDHGC5, and MT3 demonstrated increased expression levels (Figure 8).

We sourced RNA-Seq transcriptomic data from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov), comprising 535 tumor samples and 59 normal samples, for HNSC patients (Colaprico et al., 2016). Concurrently, the corresponding clinical data were retrieved from the UCSC Xena database (http://genome.ucsc.edu). A total of 289 genes associated with lactic acid metabolism (LMRGs) were collected from the GeneCards database (Stelzer et al., 2016) and Molecular Signatures Database (MSigDB) (Liberzon et al., 2015). To compare the differences in tumor immune therapy between high-risk score and low-risk score groups, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu) based on the Cox regression risk score. The copy number variation (CNV) segments underwent GISTIC 2.0 analysis using the Hiplot website (https://hiplot-academic.com/advance/gistic2). For the visualization of the protein-protein interaction (PPI) network and selecting hub genes, we utilized the Cytoscape software and the String website (Gao et al., 2013).

To assess the prognostic relevance of the Lactic Acid Metabolism Score (LMs), we utilized the R package GSVA and ssGSEA algorithm to calculate LMs based on the expression matrix of LMRGs (Hänzelmann et al., 2013). In order to categorize patients into high and low risk groups based on the expression of the five key genes, we developed a prognostic risk model using Cox regression analysis. This model incorporates the expression levels of all five genes to generate a risk score for each patient. Patients were then divided into two groups based on the median risk score: those with a risk score above the median were classified as the high-risk group, while those below the median were classified as the low-risk group. Next, we identified Differentially Expressed Genes (DEGs) associated with the high-LMs and low-LMs groups in the TCGA-HNSC dataset by employing the DESeq2 R package, adhering to stringent criteria of |logFC|>0.5 and p.adjust<0.0001 (Love et al., 2014). Furthermore, we obtained co-expression modules linked to the high-LMs and low-LMs groups using the WGCNA package (Zhang and Horvath, 2005). The pivotal genes were identified by intersecting the DEGs with LMRGs and those genes demonstrated the highest correlation with LMs within the co-expression modules. Finally, univariate Cox regression analysis was conducted to pinpoint key genes that were significantly associated with patient prognosis.

Key genes predictive of prognosis were identified through univariate Cox regression analysis. The validation of key gene expression on prognosis was conducted by scrutinizing Kaplan-Meier (KM) survival curves. Subsequently, patients were categorized into distinct molecular subtypes based on the prognostic key genes, employing the Consensus Cluster Plus R package. To analyze the data, Principal Component Analysis (PCA) was performed using the “ggplot2” R package (Ringnér, 2008).

The expression profiles of key genes were interrogated through the application of the Cox regression model. From the Cox multifactorial model, the top 20 differentially expressed genes (DEGs) associated with high and low-risk score groups were selected for further analysis. The GOSemSim package computed the top 20 DEGs linked to risk stratification in the Cox multifactorial model. Subsequently, the protein-protein interaction (PPI) network for these top 20 DEGs was constructed utilizing the String website, while the mRNA-miRNA interaction network was established using the MiRDB database (Gao et al., 2013).

The clusterProfiler package was employed to perform enrichment analysis on the DEGs. The gene set “c2.cp.all.v2022.1.Hs.symbols.gmt [All Canonical Pathways] (3050)” was acquired from the MSigDB database. Significantly enriched pathways were chosen based on stringent criteria, including a p.adjust <0.05 and a false discovery rate (FDR) value (q.value) < 0.25 (Liberzon et al., 2015).

The estimate package was utilized to assess the differences in the immune score, stromal score, and ESTIMATE score between the high-risk score and low-risk score groups, as determined by the Cox multivariate model employing the ESTIMATE algorithm. The infiltration levels of immune cells in each sample were analyzed using Single-sample Gene Set Enrichment Analysis (ssGSEA) (Barbie et al., 2009). For the analysis of gene mutations in HNSC samples derived from the TCGA database, the “MAFTOOLS” R package was employed (Mayakonda et al., 2018). Additionally, copy number variation analysis of HNSC samples from the TCGA database was carried out using the “TCGAbiolinks” R package (Mermel et al., 2011).

To evaluate the predictive efficacy of overall survival (OS), the Cox multivariate risk score was combined with patients’ clinical and pathological characteristics. Subsequently, a clinical prediction nomogram was constructed by incorporating the risk score model and clinical and pathological features. The performance of the nomogram was assessed by comparing the predicted values from the nomogram with the observed survival rates, which generated calibration curves. The calculation formula for the risk score is as follows: Risk S c o r e = ∑ i C o e f f i c i e n t   g e n e i ∗ m R N A   E x p r e s s i o n   g e n e i

TRK Lysis Buffer (20–30 mg tissue/700 µL TRK Lysis Buffer) was used to extract total RNA from 12 pairs of HNSC tumors and adjacent tissues. All of the specimens were approved by the Research Ethics Committee of the Third Xiangya Hospital of Central South University (XMXH-2022-3282). Following the instructions of the ReverTra Ace qPCR RT Kit, cDNA synthesis was carried out using the Enzyme Mix containing reverse transcriptase and the Primer Mix. KOD SYBR^® qPCR Mix was utilized for RT-qPCR. Data analysis was performed using the 2^−ΔΔCT values. The primer sequences for the hub genes are detailed in Supplementary Table 1.

The previous section has presented an extensive explanation of the statistical methods employed in processing transcriptome data. R software (Version 4.1.2) was utilized to perform all statistical analyses, with a significance level set at p < 0.05.