Identification and prognostic analysis of propionate metabolism-related genes in head and neck squamous cell carcinoma

Introduction: Head and neck squamous cell carcinoma (HNSCC) is a highly heterogeneous malignancy with poor overall prognosis. Recent studies have suggested that propionate metabolism-related genes (PMRGs) may play key roles in tumor progression and immune regulation, yet their functions in HNSCC remain unclear.

Methods: Transcriptomic data from 502 HNSCC tumor samples and 44 normal tissue samples were obtained from the UCSC Xena database as the training set. Two independent datasets (GSE41613 and GSE6631) from the GEO database were used for validation. Differentially expressed genes (DEGs), key module genes identified via weighted gene co-expression network analysis (WGCNA), and PMRGs were intersected to identify candidate genes. A prognostic model was constructed using Cox regression and LASSO analysis. Immune infiltration, somatic mutations, and drug sensitivity were compared between high- and low-risk groups. Gene expression was further validated by RT-qPCR using clinical samples.

Results: A total of 42 intersecting genes were identified, and four feature genes (PRKAA2, SLC7A5, GRIP2, CHGB) were selected to build the prognostic model. The model effectively stratified patients into high- and low-risk groups with significant survival differences in both the training and validation cohorts. The high-risk group exhibited marked differences in immune cell infiltration, immune checkpoint expression, and cancer immune cycle activity. Mutation burden and drug sensitivity also varied significantly between risk groups. A nomogram combining risk score and pathological N stage showed strong predictive performance.

Discussion: This study highlights the potential role of PMRGs in immune regulation and tumor progression in HNSCC. The proposed four-gene signature provides a novel tool for prognosis prediction and offers new insights for risk stratification and individualized therapy. Further multicenter validation and mechanistic studies are warranted.

1 Introduction

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer globally, with 5-year survival rates consistently ranging from 40% to 60% over the past few decades (1–3). Clinically, HNSCC is divided into HPV-positive and HPV-negative subtypes based on the presence of human papillomavirus (HPV), each with distinct etiologies, molecular profiles, therapeutic responses, and prognoses (4). HPV-positive HNSCC typically arises in the oropharynx, is more prevalent among nonsmokers, demonstrates relatively stable molecular features, and responds well to chemoradiotherapy, resulting in a favorable prognosis. In contrast, HPV-negative HNSCC is strongly associated with tobacco and alcohol use, displays considerable molecular heterogeneity, and is linked to poorer outcomes, including increased resistance to treatment and higher rates of local recurrence (4, 5). Current precision medicine strategies for HNSCC face two major challenges: the lack of reliable molecular biomarkers for prognostic prediction and significant individual variability in response to chemotherapeutic agents such as docetaxel and methotrexate (6). These issues highlight the urgent need for further exploration of molecular mechanisms to improve risk stratification and therapeutic approaches.

Recent research has emphasized the pivotal roles of tumor metabolic reprogramming and immune evasion. Metabolic reprogramming, for example, has been shown to influence the expression of immune checkpoint molecules such as PD-L1 (7). Tumor cells can increase PD-L1 expression through the activation of transcription factors like HIF-1α, thus suppressing T cell activity and enabling immune escape (8). Additionally, alterations in short-chain fatty acid (SCFA) metabolism, particularly propionate, have been implicated in tumorigenesis and progression (9, 10). Propionate, a key SCFA produced primarily through gut microbial fermentation of dietary fiber, not only contributes to energy metabolism but also plays pivotal roles in immunomodulation, epigenetic regulation, and cellular signaling (11). Growing evidence suggests that disturbances in propionate metabolism are closely associated with malignant progression and metastasis in various cancers (12). For instance, propionate promotes the differentiation of regulatory T cells (Tregs) and inhibits proinflammatory Th17 cells by activating G-protein-coupled receptors (GPR43/41) and suppressing HDAC activity, thereby fostering an immunosuppressive tumor microenvironment (TME) (13). Moreover, metabolites such as methylmalonic acid (MMA) can induce CD8⁺ T-cell exhaustion and enhance PD-L1 expression, further contributing to tumor immune evasion (14). In colorectal cancer and melanoma, disrupted propionate metabolism has been linked to the polarization of M2-type tumor-associated macrophages (TAMs) and the recruitment of myeloid-derived suppressor cells (MDSCs), suggesting a role in immune escape (15). Despite these findings, the biological functions and clinical significance of propionate metabolism-related genes (PMRGs) in HNSCC remain largely unexplored.

This study identified key genes associated with propionate metabolism in HNSCC and developed a prognostic model based on these genes. A comprehensive analysis of clinical features, immune cell infiltration, immune checkpoint expression, immune cycle dynamics, and drug sensitivity differences between high- and low-risk patient groups was performed. In summary, the findings of this study uncover potential therapeutic targets linked to propionate metabolism in HNSCC and offer novel insights that may aid in the development of precision treatment strategies for this challenging malignancy.

2 Materials and methods

2.1 Data source and tissues

Transcriptome sequencing data from 502 HNSCC tumor tissue samples and 44 normal tissue samples were retrieved from the UCSC Xena database (https://xenabrowser.net/datapages/) to serve as the training set. Two additional HNSCC datasets (GSE41613 and GSE6631) were sourced from the GEO database (https://www.ncbi.nlm.nih.gov/gds). The validation set included 97 oral tissue samples from patients with HPV-negative HNSCC from GSE41613 (platform GPL570). For expression verification, 22 tissue samples from patients with HNSCC and 22 normal tissue samples from GSE6631 (platform GPL8300) were utilized. A total of 603 PMRGs were obtained from the GeneCards database (https://www.genecards.org/). Real-time quantitative reverse transcription polymerase chain reaction (RT-qPCR) validation was conducted on tumors and adjacent normal tissues from 24 patients at the Department of Otolaryngology, Chongqing General Hospital. Histological evaluation was performed on each sample, and all participants provided written informed consent. The study was approved by the Ethics Committee of Chongqing General Hospital (Approval No. KY S2023-102-01).

2.2 Acquisition of intersecting genes

Gene expression data were standardized by converting probe IDs into gene identifiers and eliminating duplicate entries for the same gene in each sample to ensure a single representation per gene. Subsequently, differential expression analysis was performed using the “limma” package (v 3.58.1) (16) in the training set, identifying differentially expressed genes (DEGs) with a threshold of | Fold Change (FC)| ≥ 1 and adj. p< 0.05. Weighted gene coexpression network analysis (WGCNA) was conducted using the “WGCNA” package (v 1.70-3) (17) to identify the most relevant modules for HNSCC in the training set. Hierarchical clustering was initially performed to detect outliers, with any identified outlier samples excluded. The optimal soft threshold was determined based on the scale-free fit index (signed R²) and average connectivity (targeting a value close to 0). Genes were then grouped into modules using the hybrid dynamic tree-cutting algorithm. The correlation between these modules and the HNSCC phenotype was calculated, and the modules with the strongest correlations were defined as key modules. Genes within these key modules were identified as key module genes. Intersecting genes were derived by overlapping DEGs, key module genes, and PMRGs. To explore the biological functions and pathways involved in the intersecting genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler package (v 4.2.2) (18). Protein–protein interactions (PPI) among the intersecting genes were assessed using the STRING database (https://cn.string-db.org/).

2.3 Prognostic risk model

Univariate Cox regression analysis was performed on the intersecting genes in the training set to calculate the p-values, hazard ratios (HRs), and their 95% confidence intervals (CIs) for each gene (p< 0.05, HRs ≠ 1). Genes identified by univariate Cox regression were further analyzed using the Least Absolute Shrinkage and Selection Operator (LASSO) with the “glmnet” package (v 4.1-2) (19). Tenfold crossvalidation was conducted using the cv.glmnet function, and candidate genes were selected based on the lambda.min value that minimized the prediction error. These candidate genes were then subjected to multivariate Cox regression analysis (p< 0.05) and proportional hazards (PH) testing (p > 0.05) to identify feature genes. The risk score for each patient in the training set was calculated using the following formula: ${\displaystyle {\sum }_n^n=1(coefi*Xi)}$ . The median risk score was used to categorize the samples into high- and low-risk groups. Survival analysis was then conducted, and the Kaplan–Meier (K-M) curve was generated using the “survival” package (v 3.3-1) (20) (p< 0.05). Receiver operating characteristic (ROC) analysis was performed using the plotROC package (v 2.3.1) (21), and ROC curves for 1-, 3-, and 5-year survival were plotted, with the area under the curve (AUC) calculated (AUC > 0.7). Additionally, principal component analysis (PCA) was performed to evaluate the discriminative ability of the risk score in the training set. The same methodology was applied to validate the risk model in the validation set. The Wilcoxon test was used to assess differences in the expression of feature genes between HNSCC and control samples in both the training set and the validation set (GSE6631) (p< 0.05), with heatmaps generated to visualize the expression patterns.

2.4 Relationship between risk scores and clinical characteristics

Differential expression of feature genes across various clinical characteristics and risk groups was analyzed. The distribution of samples among each clinical characteristic group in the two risk groups was also examined. Additionally, differences in risk scores across clinical feature subgroups were evaluated, and survival differences between different risk subgroups within each clinical characteristic subgroup were computed.

2.5 Construction and evaluation of the nomogram model

Univariate and multivariate Cox regression analyses, based on risk scores, age, gender, stage, pathological T, pathological N, and grade, were performed using the “survival” package (v 3.3-1) to identify independent prognostic factors. The rms package (v 6.8-1) (21) was then employed to construct a nomogram based on the independent prognostic factors. The nomogram’s predictive performance was assessed using calibration and decision curves.

2.6 Differential expression analysis

To explore the differential gene expression between the high- and low-risk groups, differential expression analysis was performed using the DESeq2 package (v 1.34.0) (16) in the training set with the threshold set at |log₂FC| ≥ 1 and adj. p< 0.05. GO and KEGG enrichment analyses were conducted on the DEGs between the two risk groups using the clusterProfiler package (v 4.2.2) (18). Single-sample Gene Set Enrichment Analysis (ssGSEA) for KEGG pathways was performed across all samples in the training set, identifying pathways that differed between the high- and low-risk groups.

2.7 Somatic cell mutation, drug sensitivity, immune microenvironment, and immune cycle analyses

Somatic mutations in patients with HNSCC were analyzed and visualized using the maftool package (v 2.10.5). Mutation categories, types, and the frequency of the top 25 mutated genes were examined in both the high- and low-risk groups. Chemotherapeutic agents for HNSCC were obtained from the GDSC database (https://www.cancerrxgene.org). The IC₅₀ values for common chemotherapeutic and molecularly targeted drugs in each HNSCC sample were calculated using the R package pRRophetic (v 0.5) (22). Differences in IC₅₀ values between the high- and low-risk groups were assessed using the Wilcoxon rank-sum test. Subsequently, ssGSEA of 16 immune cell types, eight immune functions, 19 immune checkpoints, and seven immune cycles was performed for both groups in the training set using the GSVA package (v 1.42.0) (16). The estimate package (v 1.0.13) (18) was used to calculate stromal, immune, and ESTIMATE scores for each HNSCC sample in the training set, and differences in these scores were compared between the high- and low-risk groups.

2.8 RNA isolation, RT-PCR, semi-quantitative PCR, and qPCR

Total RNA was extracted from cell lines and tissues using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The RNA was quantified through spectrophotometry and stored at − 80°C. Primer sequences are listed in Table 1.

Table 1

Table 1. Primer sequences in the RT-qPCR experiment.

For qPCR, SYBR Green (Thermo Fisher Scientific, Hong Kong, China) was used according to the manufacturer’s instructions, with amplification performed on a 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). GAPDH served as the internal control. Gene expression levels were calculated using the 2^−ΔΔCt method, with all samples analyzed in triplicate.

2.9 Statistical analysis

Statistical analyses were performed using GraphPad Prism 9.0 (GraphPad Software Inc., CA, USA) and SPSS 23.0 (IL, USA). All experiments were conducted in triplicate, and data are presented as the mean ± standard deviation. Normality and equality of variance were assessed using the Shapiro–Wilk and Levene tests, respectively. For normally distributed data, comparisons between groups were made using Student’s t-test, with Welch’s correction for unequal variances. For non-normally distributed data, the Mann–Whitney U test was employed. The Wilcoxon rank-sum test was used to compare ssGSEA scores between groups. The Cox regression model was tested for PH assumptions, and survival analysis was conducted using the log-rank test.

3 Results

3.1 Intersecting genes were related to fatty acid metabolic processes

A total of 10,185 DEGs were identified in the training set, with 6,298 genes upregulated and 3,887 genes downregulated in HNSCC (Figures 1A, B). WGCNA revealed the green module, comprising 993 genes, as the most highly correlated with HNSCC (Cor = 0.43, adj. p = 2 × 10⁻²²) (Figure 1C). Forty-two intersecting genes were derived by overlapping the 10,185 DEGs, 993 key module genes, and 603 PMRGs (Figure 1D). GO analysis of these intersecting genes highlighted pathways such as fatty acid metabolic processes and protein-lipid complex binding (Figure 1E). KEGG pathway analysis further identified involvement in pathways such as alanine, leucine, and isoleucine degradation (Figure 1F), suggesting that these genes may influence HNSCC by modulating fatty acid metabolism. To investigate potential gene interactions, a PPI network was constructed. Genes such as ACADM and ACADS, ACHE and MAPT, as well as ACSS3 and AOX1, showed significant interactions (Figure 1G).

Figure 1

Figure 1. Acquisition of differentially expressed genes (DEGs) and key module genes. (A) Volcano plot of DEGs between HNSCC and normal samples (|log₂FC| > 0.5 and p-value< 0.05). (B) Heatmap of the top 20 DEGs between HNSCC and normal samples. (C) Correlation heatmap between gene modules and disease status. (D) Venn diagram showing the overlap of module genes, DEGs, and propionate metabolism-related genes (PMRGs) for screening PMRG-DEGs. (E) GO enrichment bubble plot. (F) KEGG enrichment bubble plot. (G) PPI network.

3.2 Prognostic risk models were constructed based on PRKAA2, SLC7A5, GRIP2, and CHGB

Univariate Cox regression analysis identified five genes (TAC1, PRKAA2, SLC7A5, GRIP2, CHGB) with p< 0.05 and HR ≠ 1 (Figure 2A). Four feature genes (PRKAA2, SLC7A5, GRIP2, CHGB) were selected through LASSO and multivariate Cox regression analysis (Figures 2B, C). Based on these feature genes, risk scores for patients with HNSCC in the training set were calculated. Patients were stratified into high- (n = 250) and low-risk (n = 251) groups based on the median risk score. As the risk score increased, mortality rates also increased (Figure 2D), with patients in the low-risk group exhibiting significantly longer survival (Log-rank test p< 0.0001) (Figure 2E). The AUC values for the 1-, 3-, and 5-year ROC curves of the risk model were all greater than 0.6, indicating strong model performance (Figure 2F). PCA demonstrated that the risk scores effectively distinguished between samples in the training set (Figure 2G). External validation in the GSE41613 dataset yielded consistent results with the training set (Figures 3A–D). The Wilcoxon test confirmed that the expression trends of feature genes in control and disease samples were consistent across both datasets, with SLC7A5 showing significant upregulation in HNSCC samples (p< 0.01) (Figures 3E–H).

Figure 2

Figure 2. Risk model construction and evaluation in the training set. (A) Forest plot of univariate Cox analysis. (B) Regression coefficient-lambda plot. (C) Forest plot of multivariate Cox analysis. (D) Risk score distribution in the training dataset. (E) Survival curves of high- and low-risk groups in the training set. (F) ROC curves for 1-, 3-, and 5-year survival based on the training set. (G) PCA dendrogram.

Figure 3

Figure 3. Validation of the risk model in the verification set. (A) Risk score distribution in the verification dataset. (B) Survival curves of high- and low-risk groups in the verification set. (C) ROC curves for 1-, 3-, and 5-year survival based on the verification set. (D) PCA dendrogram. (E, F) Heatmaps of model gene expression (training and verification datasets). (G, H) Box plots of feature gene expression in the training set and verification set (GSE6631). ns, p > 0.05; ^**p< 0.01; ^***p< 0.001.

3.3 Nomogram diagram could effectively predict the risk profile of patients with HNSCC

The expression of feature genes across different subgroups is shown in Figure 4A. The distribution of clinical characteristics in the high- and low-risk groups is presented in Figure 4B. Risk scores significantly differed between tumor grading and pathological stage T subgroups, but not between age, gender, tumor grading, tumor stage, and pathological stage N subgroups, indicating that risk scores are more closely associated with tumor grading and pathological stage T (Figure 4C). Significant survival differences between high- and low-risk groups were observed across 12 subgroups: age (≤ 60, > 60), gender (women, men), tumor grade (G2, G3), tumor stage (stage II, stage IV), pathological stage N (N1, N2), and pathological stage T (T2, T4) (Figure 4D).

Figure 4

Figure 4. Analysis of risk scores across clinical subgroups. (A) Heatmap of model gene expression across different clinical groups. (B) Distribution of clinical characteristics in high- and low-risk groups. (C) Boxplot of risk scores among different clinical characteristic subgroups. (D) Survival curves of high- and low-risk groups across different clinical characteristic subgroups.

Univariate and multivariate Cox regression analysis identified two independent prognostic factors: pathological stage N and risk score (Figures 5A, B). A nomogram was constructed based on these two factors (Figure 5C). The calibration curves for 1-, 3-, and 5-year survival showed slopes close to 1 (Figure 5D), indicating that the nomogram has high predictive accuracy. Furthermore, the 1-, 3-, and 5-year ROC curves for the nomogram demonstrated AUC values greater than 0.6 (Figure 5E), suggesting excellent prediction performance. In conclusion, the nomogram developed in this study exhibits favorable accuracy in predicting 1-, 3-, and 5-year overall survival (OS) in patients with HNSCC.

Figure 5

Figure 5. Nomogram construction and evaluation. (A) Forest plot of univariate Cox analysis. (B) Forest plot of multivariate Cox analysis. (C) Nomogram of independent prognostic factors. (D) Predicted probabilities of 1–5-year overall survival (OS) based on the nomogram. (E) ROC curves for 1-, 3-, and 5-year survival.

3.4 DEGs were related to immunity

A total of 1,336 DEGs were identified between the high- and low-risk groups, with 277 genes upregulated and 1,059 genes downregulated in the high-risk group (Figure 6A). GO analysis of these DEGs highlighted pathways such as adaptive immune response and immune system processes (Figure 6B). KEGG pathway analysis identified involvement in pathways such as primary immunodeficiency and the intestinal immune network for IgA production (Figure 6C), suggesting that these DEGs may influence risk scores through modulation of immune responses. The ssGSEA scores for seven of the 186 pathways showed significant differences between the two groups (Figure 6D), with the high-risk group exhibiting generally lower scores in these pathways.

Figure 6

Figure 6. Differential expression and pathway analysis. (A) Volcano plot of differentially expressed genes. (B) Circular plot of GO enrichment. (C) Circular plot of KEGG enrichment. (D) Heatmap of ssGSEA scores for the pathways.

3.5 High-risk high-mutation rate

In this study, 96.79% of samples in the high-risk group and 93.9% of samples in the low-risk group exhibited mutations in the top 25 most frequently mutated genes (Figure 7A). The most common mutation type was missense mutation (SNP), with TP53 showing the highest mutation frequency across samples (Figure 7B). Additionally, 12 drugs displayed significant differences in sensitivity between the high- and low-risk groups, such as dasatinib, lenalidomide, and lapatinib (Figure 7C).

Figure 7

Figure 7. Gene mutation analysis and drug sensitivity. (A) Waterfall plot of gene mutation analysis (high- and low-risk groups). (B) Gene mutation cartogram (high- and low-risk groups). (C) Boxplot of differential drug sensitivity analysis.

3.6 Significant differences in immune cells, immune checkpoints, and immune cycles

The heatmap of ssGSEA scores for 16 immune cell types is shown in Figure 8A. Except for macrophages, the remaining 15 immune cell types exhibited significant differences in scores between the high- and low-risk groups (Figure 8B). There was a negative correlation between risk scores and the scores of immune cells, with the strongest correlation observed between risk scores and CD8⁺ T cells (Figure 8C). Significant differences in immune scores and ESTIMATE scores were found between the two groups, while stromal scores showed no significant differences (Figure 8D). The ssGSEA scores for four immune functions were significantly different between the groups (Figure 8E), and 16 immune checkpoints also exhibited significant differences (Figure 8F). All seven cancer immune cycle scores differed significantly between the high- and low-risk groups (Figure 8G). Furthermore, these cancer immune cycle scores were negatively correlated with risk scores, with STEP 3 showing the strongest correlation (Figure 8H).

Figure 8

Figure 8. Immune cell and immune function analysis. (A) Heatmap of immune cell ssGSEA scores. (B) Boxplot of immune cell ssGSEA scores. (C) Lollipop diagram of correlation analysis between risk scores and immune cell scores. (D) Boxplot of stromal, immune, and ESTIMATE scores between high- and low-risk groups. (E) ssGSEA scores for immune function between high- and low-risk groups. (F) Boxplot of immune checkpoint inhibitor expression. (G) Boxplot of cancer immune cycle scores. (H) Correlation analysis between cancer immune cycle scores and risk scores.

3.7 Validation of the mRNA expression of four genes (PRKAA2, GRIP2, CHGB, SLC7A5) in HNSCC

To validate the expression changes of the feature genes in HNSCC, 24 pairs of tumor and adjacent noncancerous tissues were collected, and qPCR was performed for verification. The results showed no significant expression changes for PRKAA and GRIP2 in head and neck tumors (Figures 9A, B), whereas CHGB exhibited a noticeable upregulation, and SLC7A5 showed downregulation in head and neck tumors (Figures 9C, D).

Figure 9

Figure 9. RT-qPCR verification of model genes. (A) Expression of PRKAA. (B) Expression of GRIP2. (C) Expression of CHGB. (D) Expression of SLC7A5. Compared with Normal, **P<0.01, ***P<0.001.

4 Discussion

HNSCC is a highly heterogeneous malignancy whose development is strongly linked to HPV infection and immunometabolic reprogramming within the TME (23). Recent studies have highlighted the involvement of SCFAs, particularly propionate, as microbial metabolites that regulate energy metabolism and influence tumor progression through epigenetic modifications and immunomodulatory pathways (24). However, the exact mechanisms of PMRGs in HNSCC remain poorly understood. In this study, four characteristic genes associated with propionate metabolism in HNSCC—PRKAA2, SLC7A5, GRIP2, and CHGB—were identified through bioinformatics analysis, and their potential roles were explored, providing new theoretical insights for future research on HNSCC.

PRKAA2, also known as AMPKα2, encodes the catalytic α2 subunit of AMP-activated protein kinase (AMPK) (25). It regulates glucose metabolism, which affects tumor cell growth and energy supply (26). Notably, PRKAA2 expression is significantly elevated in hepatoblastoma (HB), where it acts as an oncogenic factor by promoting cell proliferation and inhibiting ferroptosis (27). In non-small cell lung cancer (NSCLC), PRKAA2 enhances tumor growth and suppresses ferroptosis via the SLC7A11/GSH/GPX4 pathway (28). These findings suggest that PRKAA2 may similarly influence tumor cell proliferation and survival in HNSCC.

SLC7A5 (LAT1) facilitates the cellular uptake of neutral amino acids, including leucine and glutamine (29). Its transport of leucine activates the mTORC1 signaling pathway, thereby promoting protein synthesis to support rapid tumor cell proliferation (30). Tumor cells can modulate SLC7A5 expression to alter immune cell function and evade immune surveillance (31, 32). Li et al. identified SLC7A5 as a potential prognostic biomarker in HNSCC associated with immune infiltration (33), suggesting that therapeutic targeting of SLC7A5 may offer a novel strategy for treatment.

GRIP2 encodes a PDZ domain-containing protein that binds GluR2 to anchor AMPA receptors within neuronal signaling complexes, playing pivotal roles in synaptic transmission and plasticity (34). Given the frequent dysregulation of signaling pathways in cancer cells (35), GRIP2 may influence HNSCC progression by modulating key tumorigenic pathways. Interestingly, GRIP2 has been linked to variations in innate CD8⁺ T cells (36), suggesting its potential immunomodulatory effects in HNSCC progression. Thus, GRIP2 may regulate both tumor signaling pathways and immune cell function, making it a promising therapeutic target.

CHGB is a highly conserved eukaryotic protein involved in secretory regulation (37). While CHGB genetic variants have been associated with cardiovascular disease risk (38) and the protein regulates ion channels to maintain secretory granule homeostasis (37), its role in cancer remains poorly understood and warrants further investigation.

Drug sensitivity analysis identified 12 compounds, including dasatinib, lenalidomide, and lapatinib, with significantly different IC₅₀ values between high- and low-risk groups, suggesting their potential clinical applications. Dasatinib, a multi-target tyrosine kinase inhibitor, may enhance treatment response in high-risk patients by inhibiting SRC family kinases and exerting immunomodulatory effects (39–41). The immunomodulator lenalidomide could improve the TME and increase sensitivity to chemoradiotherapy (42, 43). Lapatinib, an oral tyrosine kinase inhibitor targeting the EGFR/HER2 pathways, may provide precision therapy for specific molecular subtypes (44). These observed differences in drug sensitivities support the rationale for molecular classification and personalized treatment strategies in HNSCC. Validation through in vitro experiments and clinical cohorts is essential, alongside exploration of combination therapies with existing treatments, such as immune checkpoint inhibitors, to refine and optimize precision treatment regimens.

Significant differences in the expression of 16 immune checkpoint genes, including PD-L1, were identified between risk groups. Previous studies have shown that HNSCC cells often overexpress PD-L1, which binds to PD-1 on T cells, thereby suppressing their activation and function, enabling immune evasion (45–47). This immunosuppression is a key mechanism driving HNSCC progression (48). Furthermore, PD-L1 overexpression is associated with poorer prognosis in patients with HNSCC (49), likely due to reduced survival rates from PD-L1-mediated immune suppression. The elevated expression of PD-L1 in high-risk patients observed in this study supports these immune escape mechanisms and offers valuable insights for understanding prognostic differences and developing novel immunotherapies.

In summary, this study identified four characteristic genes associated with propionate metabolism through bioinformatics analysis and established a risk model based on these genes. These findings provide new insights for prognostic assessment and the development of innovative therapeutic strategies for HNSCC. However, several limitations must be acknowledged. The current sample size necessitates further validation through multicenter studies with larger cohorts to confirm the clinical applicability of the model. Additionally, while these metabolic genes have been identified as potential therapeutic targets, their precise mechanisms in modulating the immune microenvironment require further functional studies and clinical trials. Future research should refine this risk stratification system and investigate metabolism-targeted combination therapies to develop more precise treatment strategies for patients with HNSCC.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving humans were approved by the Ethics Committees of Chongqing General Hospital. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

SZ: Conceptualization, Funding acquisition, Writing – original draft, Writing – review & editing. YJ: Conceptualization, Investigation, Supervision, Writing – review & editing. PX: Data curation, Formal Analysis, Methodology, Resources, Software, Validation, Visualization, Writing – review & editing. ZL: Investigation, Writing – review & editing. LJ: Investigation, Writing – review & editing. WY: Investigation, Supervision, Writing – review & editing. XL: Investigation, Writing – review & editing. XA: Investigation, Writing – review & editing. JH: Investigation, Methodology, Writing – review & editing. RL: Investigation, Writing – review & editing. HM: Investigation, Writing – review & editing. HF: Investigation, Supervision, Writing – review & editing. YY: Supervision, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study was supported by a joint project between the Chongqing Science and Technology Bureau and Municipal Health Commission (2023MSXM120).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Mody MD, Rocco JW, Yom SS, Haddad RI, and Saba NF. Head and neck cancer. Lancet. (2021) 398:2289–99. doi: 10.1016/S0140-6736(21)01550-6

PubMed Abstract | Crossref Full Text | Google Scholar

2. Ruffin AT, Li H, Vujanovic L, Zandberg DP, Ferris RL, and Bruno TC. Improving head and neck cancer therapies by immunomodulation of the tumour microenvironment. Nat Rev Cancer. (2023) 23:173–88. doi: 10.1038/s41568-022-00531-9

PubMed Abstract | Crossref Full Text | Google Scholar

3. Leemans CR, Snijders PJF, and Brakenhoff RH. The molecular landscape of head and neck cancer. Nat Rev Cancer. (2018) 18:269–82. doi: 10.1038/nrc.2018.11

PubMed Abstract | Crossref Full Text | Google Scholar

4. Liu L, Xie Y, Yang H, Lin A, Dong M, Wang H, et al. HPVTIMER: A shiny web application for tumor immune estimation in human papillomavirus-associated cancers. Imeta. (2023) 2:e130. doi: 10.1002/imt2.v2.3

PubMed Abstract | Crossref Full Text | Google Scholar

5. Krsek A, Baticic L, Sotosek V, and Braut T. The role of biomarkers in HPV-positive head and neck squamous cell carcinoma: towards precision medicine. Diagnostics (Basel) 14. (2024) 14:1448. doi: 10.3390/diagnostics14131448

PubMed Abstract | Crossref Full Text | Google Scholar

6. Johnson DE, Burtness B, Leemans CR, Lui VWY, Bauman JE, and Grandis JR. Head and neck squamous cell carcinoma. Nat Rev Dis Primers. (2020) 6:92. doi: 10.1038/s41572-020-00224-3

PubMed Abstract | Crossref Full Text | Google Scholar

7. Muralidharan S, Sehgal M, Soundharya R, Mandal S, Majumdar SS, Yeshwanth M, et al. PD-L1 activity is associated with partial EMT and metabolic reprogramming in carcinomas. Curr Oncol. (2022) 29:8285–301. doi: 10.3390/curroncol29110654

PubMed Abstract | Crossref Full Text | Google Scholar

8. Shi Z, Hu C, Zheng X, Sun C, and Li Q. Feedback loop between hypoxia and energy metabolic reprogramming aggravates the radioresistance of cancer cells. Exp Hematol Oncol. (2024) 13:55. doi: 10.1186/s40164-024-00519-1

PubMed Abstract | Crossref Full Text | Google Scholar

9. Mann ER, Lam YK, and Uhlig HH. Short-chain fatty acids: linking diet, the microbiome and immunity. Nat Rev Immunol. (2024) 24:577–95. doi: 10.1038/s41577-024-01014-8

PubMed Abstract | Crossref Full Text | Google Scholar

10. Tejero J, Lazure F, and Gomes AP. Methylmalonic acid in aging and disease. Trends Endocrinol Metab. (2024) 35:188–200. doi: 10.1016/j.tem.2023.11.001

PubMed Abstract | Crossref Full Text | Google Scholar

11. Gomes AP, Ilter D, Low V, Drapela S, Schild T, Mullarky E, et al. Altered propionate metabolism contributes to tumour progression and aggressiveness. Nat Metab. (2022) 4:435–43. doi: 10.1038/s42255-022-00553-5

PubMed Abstract | Crossref Full Text | Google Scholar

12. Ramesh V, Gollavilli PN, Pinna L, Siddiqui MA, Turtos AM, Napoli F, et al. Propionate reinforces epithelial identity and reduces aggressiveness of lung carcinoma. EMBO Mol Med. (2023) 15:e17836. doi: 10.15252/emmm.202317836

PubMed Abstract | Crossref Full Text | Google Scholar

13. Du Y, He C, An Y, Huang Y, Zhang H, Fu W, et al. The role of short chain fatty acids in inflammation and body health. Int J Mol Sci 25. (2024) 25:7379. doi: 10.3390/ijms25137379

PubMed Abstract | Crossref Full Text | Google Scholar

14. Li S, Duan Y, Luo S, Zhou F, Wu Q, and Lu Z. Short-chain fatty acids and cancer. Trends Cancer. (2024) 11:154–68. doi: 10.1016/j.trecan.2024.11.003

PubMed Abstract | Crossref Full Text | Google Scholar

15. Chen L, Zhou X, Wang Y, Wang D, Ke Y, and Zeng X. Propionate and butyrate produced by gut microbiota after probiotic supplementation attenuate lung metastasis of melanoma cells in mice. Mol Nutr Food Res. (2021) 65:e2100096. doi: 10.1002/mnfr.202100096

PubMed Abstract | Crossref Full Text | Google Scholar

16. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

17. Langfelder P and Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | Crossref Full Text | Google Scholar

18. Yu G, Wang LG, Han Y, and He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

19. Ramsay IS, Ma S, Fisher M, Loewy RL, Ragland JD, Niendam T, et al. Model selection and prediction of outcomes in recent onset schizophrenia patients who undergo cognitive training. Schizophr Res Cognit. (2018) 11:1–5. doi: 10.1016/j.scog.2017.10.001

PubMed Abstract | Crossref Full Text | Google Scholar

20. Lei J, Qu T, Cha L, Tian L, Qiu F, Guo W, et al. Clinicopathological characteristics of pheochromocytoma/paraganglioma and screening of prognostic markers. J Surg Oncol. (2023) 128:510–8. doi: 10.1002/jso.v128.4

PubMed Abstract | Crossref Full Text | Google Scholar

21. Sachs MC. plotROC: A tool for plotting ROC curves. J Stat Softw. (2017) 79:2. doi: 10.18637/jss.v079.c02

PubMed Abstract | Crossref Full Text | Google Scholar

22. Gibson JM and Doniach I. Correlation of dose of x-radiation to the rat thyroid gland with degree of subsequent impairment of response to goitrogenic stimulus. Br J Cancer. (1967) 21:524–30. doi: 10.1038/bjc.1967.62

PubMed Abstract | Crossref Full Text | Google Scholar

23. Xie C, Ji N, Tang Z, Li J, and Chen Q. The role of extracellular vesicles from different origin in the microenvironment of head and neck cancers. Mol Cancer. (2019) 18:83. doi: 10.1186/s12943-019-0985-3

PubMed Abstract | Crossref Full Text | Google Scholar

24. Martin-Gallausiaux C, Marinelli L, Blottière HM, Larraufie P, and Lapaque N. SCFA: mechanisms and functional importance in the gut. Proc Nutr Soc. (2021) 80:37–49. doi: 10.1017/S0029665120006916

PubMed Abstract | Crossref Full Text | Google Scholar

25. Shariq M, Khan MF, Raj R, Ahsan N, and Kumar P. PRKAA2, MTOR, and TFEB in the regulation of lysosomal damage response and autophagy. J Mol Med (Berl). (2024) 102:287–311. doi: 10.1007/s00109-023-02411-7

PubMed Abstract | Crossref Full Text | Google Scholar

26. Zhang Q, Hong Z, Zhu J, Zeng C, Tang Z, Wang W, et al. miR-4999-5p predicts colorectal cancer survival outcome and reprograms glucose metabolism by targeting PRKAA2. Onco Targets Ther. (2020) 13:1199–210. doi: 10.2147/OTT.S234666

PubMed Abstract | Crossref Full Text | Google Scholar

27. Xie Y, Cui Z, Fang S, Zhu G, Zhen N, Zhu J, et al. Anti-ferroptotic PRKAA2 serves as a potential diagnostic and prognostic marker for hepatoblastoma. J Gastrointest Oncol. (2023) 14:1788–805. doi: 10.21037/jgo-23-110

PubMed Abstract | Crossref Full Text | Google Scholar

28. Wei Z, Zhou Z, Zhang Y, Wang J, Huang K, Ding Y, et al. PRKAA2 promotes tumor growth and inhibits ferroptosis through SLC7A11/GSH/GPX4 pathway in non-small cell lung cancer. Biotechnol Appl Biochem. (2024). doi: 10.1002/bab.2710

PubMed Abstract | Crossref Full Text | Google Scholar

29. Lu X. The role of large neutral amino acid transporter (LAT1) in cancer. Curr Cancer Drug Targets. (2019) 19:863–76. doi: 10.2174/1568009619666190802135714

PubMed Abstract | Crossref Full Text | Google Scholar

30. Bharadwaj R, Jaiswal S, Velarde de la Cruz EE, and Thakare RP. Targeting solute carrier transporters (SLCs) as a therapeutic target in different cancers. Dis 12. (2024) 12:63. doi: 10.3390/diseases12030063

PubMed Abstract | Crossref Full Text | Google Scholar

31. Branchi V, Hosni R, Kiwitz L, Ng S, van der Voort G, Bambi N, et al. Expression of the large amino acid transporter SLC7A5/LAT1 on immune cells is enhanced in primary sclerosing cholangitis-associated cholangiocarcinoma and correlates with poor prognosis in cholangiocarcinoma. Hum Pathol. (2024) 153:105670. doi: 10.1016/j.humpath.2024.105670

PubMed Abstract | Crossref Full Text | Google Scholar

32. Zhang C, Wang Y, Guo X, Wang Z, Xiao J, and Liu Z. SLC7A5 correlated with Malignancies and immunotherapy response in bladder cancer. Cancer Cell Int. (2024) 24:182. doi: 10.1186/s12935-024-03365-7

PubMed Abstract | Crossref Full Text | Google Scholar

33. Li C, Chen S, Jia W, Li W, Wei D, Cao S, et al. Identify metabolism-related genes IDO1, ALDH2, NCOA2, SLC7A5, SLC3A2, LDHB, and HPRT1 as potential prognostic markers and correlate with immune infiltrates in head and neck squamous cell carcinoma. Front Immunol. (2022) 13:955614. doi: 10.3389/fimmu.2022.955614

PubMed Abstract | Crossref Full Text | Google Scholar

34. Mao L, Takamiya K, Thomas G, Lin DT, and Huganir RL. GRIP1 and 2 regulate activity-dependent AMPA receptor recycling via exocyst complex interactions. Proc Natl Acad Sci U.S.A. (2010) 107:19038–43. doi: 10.1073/pnas.1013494107

PubMed Abstract | Crossref Full Text | Google Scholar

35. Nisar S, Hashem S, Macha MA, Yadav SK, Muralitharan S, Therachiyil L, et al. Exploring dysregulated signaling pathways in cancer. Curr Pharm Des. (2020) 26:429–45. doi: 10.2174/1381612826666200115095937

PubMed Abstract | Crossref Full Text | Google Scholar

36. Jo Y, Balmer L, Lee B, Shim JA, Ali LA, Morahan G, et al. Variants of innate CD8(+) T cells are associated with Grip2 and Klf15 genes. Cell Mol Immunol. (2020) 17:1007–9. doi: 10.1038/s41423-019-0357-3

PubMed Abstract | Crossref Full Text | Google Scholar

37. Yadav GP, Wang H, Ouwendijk J, Cross S, Wang Q, Qin F, et al. (CHGB) is dimorphic and responsible for dominant anion channels delivered to cell surface via regulated secretion. Front Mol Neurosci. (2023) 16:1205516. doi: 10.3389/fnmol.2023.1205516

PubMed Abstract | Crossref Full Text | Google Scholar

38. Zhang K, Rao F, Wang L, Rana BK, Ghosh S, Mahata M, et al. Common functional genetic variants in catecholamine storage vesicle protein promoter motifs interact to trigger systemic hypertension. J Am Coll Cardiol. (2010) 55:1463–75. doi: 10.1016/j.jacc.2009.11.064

PubMed Abstract | Crossref Full Text | Google Scholar

39. Hasinoff BB and Patel D. Mechanisms of the cardiac myocyte-damaging effects of dasatinib. Cardiovasc Toxicol. (2020) 20:380–9. doi: 10.1007/s12012-020-09565-7

PubMed Abstract | Crossref Full Text | Google Scholar

40. Lindauer M and Hochhaus A. Dasatinib. Recent Results Cancer Res. (2018) 212:29–68. doi: 10.1007/978-3-319-91439-8_2

PubMed Abstract | Crossref Full Text | Google Scholar

41. Yu GT, Mao L, Wu L, Deng WW, Bu LL, Liu JF, et al. Inhibition of SRC family kinases facilitates anti-CTLA4 immunotherapy in head and neck squamous cell carcinoma. Cell Mol Life Sci. (2018) 75:4223–34. doi: 10.1007/s00018-018-2863-3

PubMed Abstract | Crossref Full Text | Google Scholar

42. McDaniel JM, Pinilla-Ibarz J, and Epling-Burnette PK. Molecular action of lenalidomide in lymphocytes and hematologic Malignancies. Adv Hematol. (2012) 2012:513702. doi: 10.1155/2012/513702

PubMed Abstract | Crossref Full Text | Google Scholar

43. Gribben JG, Fowler N, and Morschhauser F. Mechanisms of action of lenalidomide in B-cell non-hodgkin lymphoma. J Clin Oncol. (2015) 33:2803–11. doi: 10.1200/JCO.2014.59.5363

PubMed Abstract | Crossref Full Text | Google Scholar

44. Nolting M, Schneider-Merck T, and Trepel M. Lapatinib. Recent Results Cancer Res. (2014) 201:125–43. doi: 10.1007/978-3-642-54490-3_7

PubMed Abstract | Crossref Full Text | Google Scholar

45. Tahara M, Lim DW, Keam B, Ma B, Zhang L, Wang C, et al. Management approaches for recurrent or metastatic head and neck squamous cell carcinoma after anti-PD-1/PD-L1 immunotherapy. Cancer Treat Rev. (2025) 136:102938. doi: 10.1016/j.ctrv.2025.102938

PubMed Abstract | Crossref Full Text | Google Scholar

46. Mishra PS, Sidhu A, Dwivedi G, Mulajker DS, and Awasthi S. Determining PD-L1 expression in head and neck squamous cell carcinoma using immunohistochemistry. Indian J Cancer. (2022) 59:474–9. doi: 10.4103/ijc.IJC_920_19

PubMed Abstract | Crossref Full Text | Google Scholar

47. Zheng A, Li F, Chen F, Zuo J, Wang L, Wang Y, et al. PD−L1 promotes head and neck squamous cell carcinoma cell growth through mTOR signaling. Oncol Rep. (2019) 41:2833–43. doi: 10.3892/or.2019.7053

PubMed Abstract | Crossref Full Text | Google Scholar

48. Chen Y, Ding X, Bai X, Zhou Z, Liu Y, Zhang X, et al. The current advances and future directions of PD-1/PD-L1 blockade in head and neck squamous cell carcinoma (HNSCC) in the era of immunotherapy. Int Immunopharmacol. (2023) 120:110329. doi: 10.1016/j.intimp.2023.110329

PubMed Abstract | Crossref Full Text | Google Scholar

49. Wang Q, Zhao Y, Chen Y, Chen Y, Song X, Zhang L, et al. High PD-L1 expression associates with low T-cadherin expression and poor prognosis in human papillomavirus-negative head and neck squamous cell carcinoma. Head Neck. (2023) 45:1162–71. doi: 10.1002/hed.27329

PubMed Abstract | Crossref Full Text | Google Scholar