The genomics data and clinical information of 528 HNSCC samples and 43 adjacent normal tissues were procured from the public TCGA (https://cancergenome.nih.gov/). The selection criteria were used as follows: (1) histologically confirmed HNSCC and (2) complete clinical and OS data. Lastly, 479 patients with the corresponding clinical information, including age, gender, stage, HPV subtype, and radiation therapy were collected for further analysis. The mutation data (e.g., somatic mutation and copy number variation data) was downloaded from the UCSC Xena (https://gdc.xenahubs.net/). Twenty-one m⁶A regulators were collected based on published literature. Next, the differential expression of the 21 m⁶A regulators was presented in a heatmap. Nonsynonymous mutation and synonymous mutation counts were defined as tumor mutation burden. The GSE65858 (N = 267) from GEO was used as the validation cohort. The detailed information of clinical data and 21 m⁶A regulators are shown in Supplementary Tables S1 – S3 .

To elucidate the biological function of the m⁶A regulators in HNSCC, ConsensusClusterPlus package based on Euclidean distance and Wards linkage was employed to classify the patients into different distinct m⁶A subtypes (31). The “PCA” package was used to investigate gene-expression arrays among distinct m⁶A subgroups.

We utilized the gene set variation analysis (“GSVA”) package to investigate the biological processes among different m⁶A subgroups (32). The well-defined biological pathways and functions were derived from the Hallmarker gene set “c2.cp.kegg.v7.4.symbols.gmt” and “c5.go.v7.4.symbols.gmt” (download from MSigDB database v7.4) and IMvigor210CoreBiologies package (33, 34). The “ClusterProfiler” package was used to determine the Gene Ontology (GO) annotation of m⁶A-related genes (the cutoff value were q-value <0.05 and p-value <0.05) (35).

The relative abundance and activity levels of 23 immune cell types, obtained from published signature gene lists, were quantified using the single sample gene set enrichment analysis (ssGSEA) in R package GSVA (36, 37). In this study, the innate immune cells (including natural killer (NK) cells, CD56dim NK cells, CD56bright NK cells, dendritic cells (DCs), plasmacytoid dendritic cells (pDC), immature DCs (iDC), neutrophils, mast cells, and macrophages) and the adaptive immune cells (including B cells, T cells, CD8 T cells, T follicular helper (TFH), Th1, Th2, Th17, and Treg cells) comprised these signatures. In addition, we also used ssGSEA to explore the relationship between different m⁶A subtypes and immune-related pathways (such as cytolytic activity, T-cell co-stimulation, inflammation promoting, and parainflammation) in HNSCC expression profile of TCGA. The biosimilarity of the infiltrating immune cells and immune-related functions were estimated by the Gaussian fitting model.

Immunophenoscore (IPS) is an effective predictor of response to immune therapy via characterizing the determinant factors of cancer immunogenicity and antigenomes (37). The major histocompatibility complex (MHC)-related molecules, checkpoints or immunomodulators (CP), effector cells (EC), and suppressor cells (SC) developed the IPS scoring scheme. The sum of the four classes, calculated by averaging the Z-scores, was defined as the IPS. To predict immune checkpoint blockade response (ICB), we utilized the tumor immune dysfunction and exclusion (TIDE) method to model tumor immune evasion mechanisms, including the dysfunction of T-cell dysfunction in tumors with high infiltration of cytotoxic T lymphocytes (CTLs) and the prevention of T cell in tumors with exclusion of CTLs (38). For patients with higher TIDE score, cancers more likely to occur immune escape in these patients' body, thus ICB treatment might bring these patients less and short-lasting clinical benefits. The ESTIMATE algorithm was used to evaluate the tumor cellularity and tumor purity, which were composed of the TIME, based on expression matrixes. The analysis method is integrated in the “ESTIMATE “ package (39). We extracted these gene expression data from RNASeqV2 data to predict different infiltration levels of immune cells and the proportion of stromal cells. Tumor purity is the summation of stromal score and immune score from individual cases. The tumor sample with higher immune scores and lower tumor purity indicated that it had an abundance of immune cell infiltration.

The mutation annotation format (maf) file was analyzed using MutSigCV algorithm to identify significant SMGs based on the significance threshold, and the maf data were processed using the “maftools” package (40). MutSigCV measures the significance of nonsilent somatic mutations in a gene based on the background mutation rates by silent mutations (41). The false discovery rates (q-values) were then calculated, and genes with statistical significance (q-values ≤0.1) were set as SMGs ( Supplementary Table S4 ). We then utilized the waterfall plot to visualize the mutation information of these significant SMGs in the TCGA cohort. Furthermore, we applied Fisher’s test to detect mutually exclusive or co-occurring ratio of m⁶A regulators. Mutational signatures were determined using the genomic data by adopting ExtractSignatures function that applies the Bayesian nonnegative matrix factorization-based framework (42). The optimal number of mutational signatures for the TCGA cohort could be detected by the SignatureEnrichment function and then it automatically assigned a given signature to each sample.

Patients were grouped into the three m⁶A clusters based on consensus clustering algorithm to identify differentially expressed genes (DEGs) associated with the m⁶A modification. The “limma” package was implemented to determine DEGs between three m⁶A clusters (43). The significance filtering cutoff of DEGs were set as the significance-adjusted p-value <0.001.

The overlapped DEGs identified from DEGs were used to perform the univariate Cox regression. The consensus clustering algorithm was utilized to define the number of gene clusters. The prognosis-related genes were extracted for further analysis. We then curated the final genes determined to conduct principal component analysis (PCA), and principal component 1 and 2 were extracted to construct the m⁶A score (44, 45). This method has an advantage of mainly focusing on positively correlated (or negatively correlated) genes. We then define the m⁶A score of each patient by adopting a similar formula based on the previous studies:

To determine the TMB of each patient, we also counted the nonsynonymous and synonymous mutation counts in the TCGA cohort (46). The association with TMB and m⁶A score was evaluated by Spearman’s method based on survival curve.

Mariathasan et al. constructed a panel of signatures that stored genes associated with various biological pathways, including (1) immune-checkpoint; (2) CD8 T-effector signature; (3) epithelial-mesenchymal transition (EMT), including EMT1, EMT2, and EMT3; (4) pan-fibroblast TGFb response signature (Pan-F-TBRS); (5) Fanconi anemia pathway; (6) homologous recombination; (7) base excision repair; (8) WNT target; (9) DNA damage repair; (10) mismatch repair; (11) nucleotide excision repair; (12) DNA replication; (13) antigen processing; (14) cell cycle regulation; (15) FGFR3-related genes; and (16) cell cycle (34, 47). We performed the Spearman’s method to explore the correlation between m⁶A score and these biological pathways.

We systematically performed a search for the ICB cohorts in the public databases, which could be available for detailed genomic and clinical information. Three independent anti-PD-L1 cohorts, IMvigor210 cohort (patients with metastatic urothelial cancer treated with atezolizumab) (34), Riaz et al. cohort (patients with metastatic melanoma treated with nivolumab) (48), and GSE78220 cohort (patients with metastatic melanoma treated with pembrolizumab) (49), were finally downloaded to analyze the predictive value of the m⁶A score for immunotherapy. The raw gene expression data of all cohorts were normalized.

We used the largest public pharmacogenomics database, Genomics of Drug Sensitivity in Cancer (GDSC), to predict the sensitivity of different drugs between high– and low–m⁶A score subgroups (50). The prediction process used was the “pRRophetic” package where the half-maximal inhibitory concentration (IC₅₀) was estimated by ridge regression model based on gene expression profiles (51).

The statistical analyses were generated by using R version 4.1.0. To compare more than two groups, statistical significance was estimated by the Kruskal-Wallis test. Student’s t-test was used to compare the difference between two subgroups (52). Kaplan-Meier analysis generated the differences between m⁶A subgroups and prognosis via the “survminer” package. To determine the optimal cutoff values of each cohort, we used the “surv-cutpoint” function from the “survival” package. We adopted Cox regression to calculate the hazard ratios (HR) of m⁶A regulators and m⁶A-related genes. The multivariate Cox regression was used to evaluate the independent prognostic factors. The “forestplot” package was employed to show the results of Cox regression analysis for m⁶A score in the GEO cohort and TCGA cohort. We assessed the specificity and sensitivity of m⁶A score through drawing receiver operating characteristic (ROC) curve by using “pROC” and “‘timeROC” package. Also, the Spearman’s method was used to compute the correlation coefficient. All comparisons were presented by the p-values (two-tailed), whereby <0.05 indicated statistical significance.

We firstly identified 21 m⁶A regulators (including eight “writers,” 11 “readers,” and two “erasers”) in the TCGA cohort. Figure 1A and Supplementary Figure S1A summarize the significant biological processes and functions of 21 m⁶A regulators conducted by Metascape database. Then, the waterfall plot presented the incidence of copy number variations and the ratio of somatic mutations of 21 m⁶A regulators. A total of 72 of the 479 (15.03%) patients experienced mutations, mainly including missense mutation, splice site, and nonsense mutation. In Figure 1B , we found that KIAA1429 exhibited the highest mutation frequency, followed by LRPPRC and YTHDC2, while YTHDC1, YTHDF2, IGF2BP2, HNRNPC, METTL14, and RBM15B did not show any mutations. The results of mutation co-occurrence examined the significant relationship between IGF2BP3 and FTO, RBM15 and YTHDF1, LRPPRC and YTHDF2 ( Supplementary Figure S1B ). Further investigation revealed the CNV frequency of 21 m⁶A regulators. Most m⁶A regulators showed the prevalent deletions in copy number, while IGF2BP2, YTHDC1, and CBLL1 had a widespread frequency of CNV amplification ( Figure 1C ). Figure 1D shows the location of CNV of all m⁶A regulators on chromosomes. We further demonstrated that the expressions of ALKBH5, METTL3, YTHDF2, and YTHDC2 were significantly downregulated in tumor samples, and in contrast the expression of CBLL1, METTL14, IGF2BP2, IGF2BP3, KIAA1429, YTHDF1, and YTHDC1 were significantly upregulated in tumor samples ( Figure 1E ). Compared with normal tissues, m⁶A regulators (such as CBLL1 and YTHDF1) with amplificated CNV demonstrated markedly higher expression, and YTHDF2 and YTHDC2 with prevalent CNV deletions were markedly decreased in the tumor ( Figures 1C, E ). Spearman’s method presented the correlation among these m⁶A regulators ( Supplementary Figure S1C ). We found that IGF2BP2 showed no significant correlation with some m⁶A regulators (RBM15B, YTHDC2, RBM15, YTHDF2, and METTL14). We then ascertain the prognostic value of 21 m⁶A regulators using the Cox regression. The Cox regression revealed that YTHDC2 was a protective factor, significantly associated with prolonged overall survival rate, while HNRNPA2B1 was a risk factor ( Supplementary Figures S1D, E ). Based on these results, we demonstrated that m⁶A regulators had significant heterogeneity of genomic and transcriptomic alteration landscape between normal and HNSCC samples.

The TCGA dataset with available survival and clinical information were enrolled into the training cohort. The regulator network comprehensively depicted the whole interactions of 21 m⁶A regulators and their prognostic significance ( Figure 2A ). We found that not only eraser genes were all risk factors, while some of the writer and reader genes were favorable factors. Moreover, we demonstrated that the connection among 21 m⁶A regulators were positively correlated. These results indicated that cross-talk among the 21 regulators probably play critical roles in the formation of different m⁶A modifications and pathogenesis and progression in individual tumors. Based on the hypotheses, we utilized unsupervised clustering to classify samples into different m⁶A clusters. Moreover, we could completely distinguish one m⁶A cluster from other clusters based on PCA ( Figure 2B ). Accordingly, three distinct m⁶A clusters were eventually identified, including 128 cases in m⁶A cluster A, 247 cases in m⁶A cluster B, and 121 cases in m⁶A cluster C ( Figure 2C ; Supplementary Figures S2A, B ).

Among these clusters, m⁶A cluster A, m⁶A cluster B, and m⁶A cluster C, patients in m⁶A cluster A had an advantage in overall survival rate, whereas m⁶A cluster C revealed the poorer prognosis in the TCGA cohort (p = 0.022). In the validation cohort (GEO cohort), the identical analyses obtained similar results (p = 0.049, Figure 2D ; Supplementary Figure S2C ).

In the TCGA cohort, multivariate Cox regression further demonstrated that patients in m⁶A cluster C had worst overall survival rate after adjusting clinical parameters [m⁶A cluster C vs. m⁶A cluster A, HR, 1.68 (95% CI, 1 to 2.8), p = 0.049, Supplementary Figure S4A ]. However, there was no statistical significance between m⁶A cluster C and prognostic outcome in the GEO cohort [m⁶A cluster C vs. m⁶A cluster A, HR, 1.47 (95% CI, 0.88 to 2.47), p = 0.143, Supplementary Figure S4B ]. We also noticed that the 21 m⁶A regulators showed different significances between the three m⁶A clusters. In detail, KIAA1429 and FTO were significantly elevated in m⁶A cluster A; CBLL1, IGF2BP2, and IGF2BP3 were significantly elevated in m⁶A cluster B; and WTAP, ALKBH5, and RBM15 were significantly elevated in m⁶A cluster C ( Supplementary Figures S2B, C ).

To explore the biological functions and pathways underlying these m⁶A clusters, we performed GSVA enrichment analysis against the GO and KEGG gene sets ( Supplementary Figures S3A, B ). As shown in the GSVA analysis, m⁶A cluster A was markedly enriched in immune activation-related pathways. Intriguingly, m⁶A cluster C was markedly associated with carcinogenic pathways, such as DNA replication, nucleotide excision repair, and mismatch repair pathways. Whereas, m⁶A cluster B was highly enriched in both carcinogenic and stromal-related signaling pathways.

The heatmap visualized the infiltration levels of 23 immune cells among three m⁶A clusters ( Figure 3A ). Antitumor lymphocyte cells, such as activated CD8+ T cells, and NK cells, were mainly enriched in the m⁶A cluster A. However, regulatory T cells and type 1/2/17 T helper cells were mainly enriched in the m⁶A cluster B. To our surprise, innate immune cells including natural killer cell, macrophage, eosinophil, mast cell, and MDSC were increased in the m⁶A cluster C. To explore the subsets of immune cell in TIME, CIBERSORT package was further used to characterize the immune cell infiltration based on the expression file. We observed the consistent result in the Figure 3B . Previous studies revealed a novel immune phenotype, immune-excluded phenotype, with an abundance of immune cells, retained in the tumor stroma rather than in the parenchyma. Therefore, we speculated that the m⁶A cluster B with higher stromal score exhibited an ineffective antitumor immune response ( Figure 3E ). Cancer-related pathway analyses demonstrated that the m⁶A cluster B was related to TGF-β and WNT-target pathways, which further corroborated with our hypothesis ( Supplementary Figures S4C, D ). In Figure 3C , we found that m⁶A cluster A exhibited the highest immune scores, followed by m⁶A cluster B and m⁶A cluster C. Conversely, m⁶A cluster C had a higher tumor purity than m⁶A cluster B and m⁶A cluster A, suggesting that tumors in m⁶A cluster B and m⁶A cluster A are surrounded by more immune cells and stromal cells ( Figure 3D ).

Then, we examined the association between 21 m⁶A regulators and immune cells via Spearman’s method. We focused on the regulator HNRNPA2B1, an independent prognostic risk factor based on the above results ( Supplementary Figures S1D, E ), which was negatively correlated with numerous immune cells ( Supplementary Figure S5A ). The ESTIMATE showed that low-expression subgroup of HNRNPA2B1 exhibited higher immune score, which confirmed the above findings ( Supplementary Figure S5B ).

We also found that low-expression subgroup of HNRNPA2B1 exhibited a significant increased among 23 immune cells ( Supplementary Figure S5C ). The low-expression subgroup of HNRNPA2B1 also exhibited elevated expression of HLA molecules ( Supplementary Figure S5D ). Subsequent function enrichment analyses found that low-expression subgroup of HNRNPA2B1 exhibited an obvious enhancement in immune activation including T-cell costimulation and type I/II IFN responses, which hinted that the expression of HNRNPA2B1 might affect the efficacy of immunotherapy ( Supplementary Figure S5D ). Thus, we investigated two anti-PD-L1 immunotherapy cohorts (IMvigor210 cohort and GSE78220 cohort). In the IMvigor210 cohort, patients with low expression of HNRNPA2B1 had prolonged overall survival rate ( Supplementary Figure S5E ). In the GSE78220 cohort, there was no significant survival trend ( Supplementary Figure S5F ). Therefore, we speculated that HNRNPA2B1-mediated m⁶A methylation modification may enhance the antitumor response via promoting the activation of immune cells.

To identify the biological behaviors (e.g., genetic alterations and expression perturbations) of these m⁶A clusters, we fixed attention on the m⁶A-related transcriptional expression alterations across three m⁶A clusters in HNSCC. The Venn diagram determined 4,269 overlapping differentially expressed genes (DEGs) ( Figure 4A ). A total of 311 DEGs related to prognosis were considered the representative m⁶A-related genes ( Supplementary Table S5 ). GO enrichment analysis revealed that the biological processes related to RNA transcription and modification were significant functions ( Figure 4B ). Similar to the above analysis, unsupervised clustering method based on the expression of these 311 DEGs separated patients into three stable gene clusters (gene clusters A–C) in the TCGA cohort ( Supplementary Figure S6A ). Figure 4C demonstrates that three m⁶A gene cluster had different clinicopathological features. We found that patients in m⁶A gene cluster C exhibited advanced clinical stage. In addition, patients receiving radiotherapy were mainly concentrated in the m⁶A gene cluster A, while patients with negative HPV subtype were represented by the m⁶A gene clusters B and C.

The survival analysis further indicated that the three m⁶A gene clusters had significant prognostic differences in HNSCC samples. m⁶A gene cluster A was proven to be related to better prognostic outcome, while patients in m⁶A gene cluster C was associated with poorer outcome ( Figure 4D ). The Cox regression determined m⁶A gene cluster C (vs. m⁶A gene cluster A) as an independent risk factor after considering age, gender, stage, HPV subtype, and radiotherapy [HR, 1.52 (95% CI, 1.01 to 2.28), p = 0.045; Figure 4E ]. Supplementary Figure S6B observes the different expression levels of the 21 m⁶A regulators, which were consistent with our expected results.

Accordingly, the above results showed that the m⁶A regulators played a nonnegligible role in regulating prognosis and TIME. However, these analyses were only based on the overall population and could not interpret the heterogeneity and complexity of m⁶A regulators individually. Based on these identified m⁶A-related genes, we developed a scoring scheme, considered m⁶Ascore, to quantify individual patients.

The alluvial diagram visualized the quantification changes of patients ( Figure 5A ). These results illustrated that m⁶A gene clusters B and C were linked to higher m⁶A score, whereas m⁶A gene cluster A exhibited lower m⁶A score. Notably, m⁶A cluster C showed the highest m⁶A score, followed by m⁶A cluster B, while m⁶A cluster A revealed the lowest m⁶A score ( Supplementary Figure S7A ). Furthermore, we conducted the analysis of Spearman’s correlation to illustrate the patterns of m⁶A regulators. The heatmap indicated that m⁶A score was positively correlated with WNT target signatures, cell cycle signatures, and EMT pathways ( Figure 5B ).

There was an inverse trend between the m⁶A score and the immune score (R = −0.35, p = 2.3e−16) and stromal score (R = −0.08, p = 0.076), which demonstrated the crosstalk between m⁶A score and TIME ( Supplementary Figure S7B, C ). Compared with the high m⁶A score, patients in low–m⁶A score subgroup had higher relative level of immune checkpoint molecules and CD8 effector cells. However, high m⁶A scores were significantly associated with stromal pathways ( Supplementary Figure S7D ).

Furthermore, we determined the prognostic value of m⁶A score in predicting patients’ survival outcome. Based on the cutoff value of 3.3615, we divided patients into low– or high–m⁶A score subgroups. As expected, patients with low–m⁶A score were associated with a prominent prognosis (p < 0.001, Figure 5C ), and the ROC validated the predictive accuracy of the m⁶A score (AUC = 0.634, Supplementary Figure S7E ). Integrating the clinical information (e.g., age, gender, stage, HPV subtype, and radiotherapy), multivariate Cox regression confirmed that the high m⁶A score was an independent prognostic factor for evaluating survival outcome (high m⁶A score vs. low m⁶A score; HR, 0.61 [95% CI, 0.44 to 0.86], p < 0.01, Supplementary Figure S7G ).

We also investigated the relationship between the m⁶A score and level and found that the expression level of PD-L1 was elevated in the low–m⁶A score subgroup than in the high m⁶A score subgroup ( Figure 5E ). The constructed m⁶A score was validated in the GEO cohort by integrating clinical genomic information. The m⁶A score displayed the potential predictive value in GEO cohort (AUC = 0.672, Supplementary Figure S7F ), and patients in low–m⁶A score subgroup had a better survival outcome (p = 0.044; Figure 5D ). Multivariate Cox regression also confirmed that the m⁶A score was an independent prognostic biomarker in GEO cohort [HR, 0.52 (95% CI, 0.3 to 0.92), p = 0.024; Supplementary Figure S7H ]. We then further analyzed the distribution of somatic mutated gene between low– and high–m⁶A score subgroups. As shown in Figure 5G , high m⁶A score subgroup presented more tumor somatic mutations than the low–m⁶A score group. Increasing studies have demonstrated and there was a link between the TMB and immunotherapy responses. Consequently, we further explored the distribution of TMB in low– and high–m⁶A score subgroups. We confirmed that the low–m⁶A score group had lower TMB frequencies ( Figure 5F ). The m⁶A score was markedly positively correlated with TMB (R = 0.16, p = 0.00041; Figure 5H ). In addition, we found that patients with low TMB frequencies demonstrated a survival benefit (p = 0.012; Figure 5I ), while patients with low m⁶A score showed significant survival advantages among patients with low TMB frequencies ( Figure 5J ).

The prognostic value of m⁶A score subjected to various clinical parameters was also estimated. We found that patients in low m⁶A score had a better survival outcome than those in m⁶A score among different subgroups ( Supplementary Figure S8 ). Furthermore, the OS of patients with radiotherapy in the high– and low–m⁶A score groups was superior, but patients with low m⁶A score benefited significantly more than those with high m⁶A score from radiotherapy. Accordingly, patients with low m⁶A score were more likely to benefit for survival from radiotherapy than those with high m⁶A score.

TIDE and IPS served as novel imunotherapeutic predictors and are strongly suggested to evaluate the response of immunotherapy to patients. We revealed that TIDE was significantly decreased in the low–m⁶A score subgroup, and IPS was significantly elevated in the low–m⁶A score subgroup (IPS: p = 0.0014, Supplementary Figure S9A ; TIDE: p = 0.0035, Supplementary Figure S9B ). In detail, the levels of the four groups were significantly increased in the low m⁶A score group ( Supplementary Figures S9C–F ).

We investigated whether the m⁶A score could predict immunotherapy response to ICB treatment based on three cohorts. Among IMvigor210 cohort and Riaz et al. cohort, patients with low m⁶A score exhibited clinical benefits markedly (IMvigor210 cohort, p < 0.001, Figure 6A ; Riaz et al. cohort, p = 0.048, Figure 6J ). In GSE78220 cohort, the survival curve presented an opposite result due to the small number of samples ( Figure 6H ). The immunotherapeutic advantages and anti-PD-1/L1 response to patients were confirmed in the low– and high–m⁶A score subgroups ( Figures 6B, C, I, K ).

We investigated the difference between m⁶A score and three immune phenotypes in the IMvigor210 cohort and found that lower m⁶A score was remarkably associated with excluded immune phenotype, indicating that immune checkpoint inhibitors are less effective for these patients ( Figure 6D ). Furthermore, we revealed that EMT were significantly activated in tumors with low m⁶A score ( Figure 6E ). Figure 6F indicates that individuals with a combination of high m⁶A score and low neoantigen burden showed a poorer prognosis. The ROC curves implied that m⁶A score was a robust biomarker to assess clinical prognosis of patients under immunotherapy ( Figure 6G ). In summary, our work strongly demonstrated that m⁶A regulators was significantly correlated with TIME and mediated prognostic response to immunotherapy.

Considering the frequent use of chemotherapy in the treatment of HNSCC, we further explored the response of patients with 138 different types of drugs. In detail, the GDSC dataset was used to predict the IC₅₀ of the selected drugs based on the “pRRophetic” package. A total of 54 drugs demonstrated obviously lower IC₅₀ in the low–m⁶A score group ( Supplementary Figure S10 ). Based on the guidelines of the National Comprehensive Cancer Network (53) and Chinese Society of Clinical Oncology (54), we summarized all the drugs used for the treatment of head and neck tumor (including cisplatin, methotrexate, cetuximab, afatinib, capecitabine, oxaliplatin, carboplatin, docetaxel, nivolumab, camrelizumab, gemcitabine, nimotuzumab, 5-FU, paclitaxel, pembrolizumab, toripalimab, and nedaplatin). However, only paclitaxel presented obvious lower IC₅₀ in the low–m⁶A score group. The finding suggested that patients with low m⁶A score were more sensitive to the treatment of paclitaxel than those with high m⁶A score in HNSCC.