N1-methyladenosine (m 1 A) regulator-mediated methylation modification modes and tumor microenvironment infiltration characteristics in head and neck squamous cell carcinomas

: Background: Recent researches have investigated the biological importance of RNA N1-methyladenosine (m 1 A) modifications in oncogenesis and progression of head and neck squamous cell carcinoma (HNSCC). However, whether m1A modifications also have latent effect in tumor microenvironment (TME) generation and immune regulation in HNSCC is unknown. Methods: We evaluated the m 1 A modification patterns and related to these modification patterns with TME cell infiltration features in 1041 HNSCC samples by bioinformatics approach. Results: The m 1 A score is an independent prognostic indicator. HNSCC patients with low m 1 A score group with poor overall survival (OS) was mainly characterized by stroma activation, lack of sufficient immune infiltration, and exhibited an immune- desert TME phenotype. Low m 1 A scores were also correlated with increased tumor mutation burden (TMB), and HNSCC patients with high TMB and low m 1 A scores had the worst OS. In addition, anti-CTLA-4 combined with anti-PD1 treatment was more effective in the high m 1 A score subgroup than in the low m 1 A score subgroup. Conclusions: This study revealed that m 1 A modifications play a non-negligible role in developing the TME versatility and complexity of HNSCC. Assessing m 1 A modification patterns in HNSCC helps improve our comprehension of its TME infiltration profile and guides more effective immunotherapeutic approaches. a highly synergistic manner. Consequently, an in-depth understanding of the regulation of TME cell infiltration properties by multiple m 1 A modulators can contribute to a better comprehension of TME immune regulation. This present research comprehensively evaluates the relationship between m 1 A modification modality and TME cell infiltration traits via combining transcriptional and genomic data of 1041 HNSCC cases from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. We constructed two diverse m 1 A modification patterns. Surprisingly, TME characteristics were very similar to the formerly reported


Introduction
Head and neck squamous cell carcinomas (HNSCC), develop from the mucosal epithelium of larynx, pharynx, and oral cavity, are the most ordinary malignancies of the head and neck [1]. HNSCC is the sixth most ordinary cancer worldwide, with a high mortality and morbidity rate, nearly 890,000 new cases, and 450,000 deaths in 2018[2,3]. According to the Global Cancer Observatory estimate (IARC https://gco.iarc.fr/today/home): The incidence of HNSCC was sustainable rising and is expected to enhance by 30% (or 1.08 million new cases per year) in 2030 [2,4]. These risk factors include exposure to environmental pollutants, smoking, alcohol consumption, and infection with viruses, namely EBV and HPV. HNSCC is featured through genetic instability and frequent gain or loss of chromosomal regions [5]. HNSCC is mainly treated with surgical resection, followed by adjuvant radiation or chemotherapy combinations of radiation (known as chemoradiotherapy (CRT) or chemoradiation) according to cancer stage. Survival rate of HNSCC has enhanced over the past thirty years. The enhanced survival rate is partly ascribed to the appearance of HPV-related HNSCC, a population with an improved prognosis rather than an improvement in multimodal therapy. However, some types of cancer have stagnant survival rates, such as laryngeal cancer [1]. Thus, there is an urgent need to investigate the molecular mechanism of HNSCC and determine new approaches to treat and intervene in HNSCC patients.
Epigenetics plays an essential role in the progression of cancer. It is mainly involved in DNA methylation, non-coding RNA, nucleosome remodeling and histone modifications, etc [6]. Many RNA modifications have been recognized in diverse RNAs. N1methyladenosine (m 1 A) was identified in RNA isolated from rat liver in the early 1960s [7] and later it was discovered to be present in tRNA [8], rRNA [9], and mRNA [10]. Until a few years ago, the function of m 1 A in eukaryotic rRNAs was virtually untouched. Nevertheless, its diverse biological functions have only recently begun to be developed. In human, the modification of m 1 A at position 58 of tRNA molecule is a dynamic and varies process that is written through the TRMT6/TRMT61 heterodimer ("writers"), and demethylated through ALKBH1 ("erasers") [11]. Several previous studies have indicated that aberrant expression and genetic modifications of the m 1 A regulator are related to disruptions in a variety of biological processes, including malignant progression of cancer and abnormalities in immune regulation [12][13][14]. A more profound comprehension of the genetic variants and expression disorders behind tumor heterogeneity will further recognize therapeutic targets based on m 1 A RNA methylation.
In recent years, with the enhancing comprehension of the distinct and complexity of the tumor microenvironment (TME), the vital role of immune cell subsets in cancer genesis and metastasis has been gradually recognized [15,16]. The tumor component consists of TME, including stromal cells and cancer cells, like resident fibroblasts and macrophages, infiltrating immune cells, bone marrow-derived cells, secreted factors, etc. [17]. Immunotherapies represented as immune checkpoint inhibitors (ICIs), for example, programmed cell death-1/ ligand-1(PD-1/L1) and cytotoxic T-lymphocyte-related antigen-4 (CTLA-4), have durable and practical clinical efficacy in some cancer patients. Unfortunately, most cancer patients have little or no clinical benefit and fall far short of clinical need [18][19][20]. However, ICIs, such as anti-PD-1/L1, have transformed the treatment of HNSCC. FDA granted the first immunotherapy approvals -for the anti-PD-1 ICIs pembrolizumab and nivolumab -for the treatment of recurrent HNSCC in 2016, based on results of KEYNOTE-012 [21] and for first-line treatment of metastatic or inoperable HNSCC in 2019, according to the results of KEYNOTE-048 [22,23], which changed the face of HNSCC treatment and clinical trial opportunities. Estimating immune infiltration based on TME signatures is critical to evaluate the response to current ICIs and explore new immunotherapeutic approaches [24]. Therefore, identifying different tumor immune phenotypes through a comprehensive assessment of the heterogeneity and complexity of the TME environment could significantly improve the ability to guide and predict immunotherapeutic responses. In addition, several biomarkers may be identified that will prove to be highly effective in identifying patient responses to immunotherapy and will contribute to identify new therapeutic targets [25].
Recent researches have elucidated the interaction between m 1 A modification and TME immune cell infiltration, but RNA degradation mechanisms cannot fully interpret it. The deletion of YTHDF1 in typical dendritic cells increases the cross-presentation of tumor antigens and cross-stimulation of CD8+ T cells in vivo. Inhibition of YTHDF1 also improves the therapeutic effect of anti-PD-L1 receptor blockers [26]. However, because of technical constraints, the above researches are indeed limited to one or two m 1 A modulators and cell types, whereas the antitumor action is characterized by some anti-cancer inhibitors interacting in a highly synergistic manner. Consequently, an in-depth understanding of the regulation of TME cell infiltration properties by multiple m 1 A modulators can contribute to a better comprehension of TME immune regulation.
This present research comprehensively evaluates the relationship between m 1 A modification modality and TME cell infiltration traits via combining transcriptional and genomic data of 1041 HNSCC cases from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. We constructed two diverse m 1 A modification patterns. Surprisingly, TME characteristics were very similar to the formerly reported immune phenotypes: immune excluded and immune inflammation, respectively. [27]. Furthermore, we developed a scoring system to quantify the pattern of m 1 A modifications in individual patients and forecast patients' clinical response to ICI therapy. This finding suggests that m 1 A modifications play a significant role in forming a personalized TME profile for HNSCC and guiding therapeutic intervention strategies. The workflow of our study was manifested in Figure 1.

Data sources and screening process
The mRNA expression levels, somatic mutation data, and clinical information (including survival status, overall survival time, age, gender, grade, stage, TNM stage, smoking status, alcohol consumption, and Human papillomavirus (HPV) infection status) of HNSCCC were retrospectively obtained from the TCGA (https://www.cancer.gov/aboutnci/organization/ccg/research/structuralgenomics/tcga) database and NCBI-GEO database (https://www.ncbi.nlm.nih.gov/geo/ . We download the series matrix file(s) and platforms files (GPL10150 for the GSE47443 cohort and GPL10558 for the GSE65858 cohort), then obtain the gene expression and clinical information using a Perl script. We used the R-package "limma" to merge the gene expression data of TCGA-HNSCC, GSE47443, and GSE65858 cohorts, then rectify the batch impact from non-biological technical biases by using the "ComBat" algorithm of the R-package "sva" [32]. The baseline information of TCGA-HNSCC, GSE47443, and GSE65858 cohorts was compiled in Table S1. The copy number variation (CNV) of TCGA-HNSCC was downloaded from the UCSC Xena database (http://xena.ucsc.edu/) [33].

Comprehensive analysis of nine m 1 A regulators
The R package "maftools" is applied to analyze and visualize Mutation Annotation Format (MAF) files of somatic mutation data. We obtained the CNV frequencies of the nine m 1 A regulators in each HNSC patient from the CNV file of HNSC by a perl script. The R package "Rcircos" was used to visualize the location of nine m 1 A regulators in human chromosomes. Besides, we used the GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) a web-based analysis platform for gene set cancer analysis [34], to explore the m 1 A regulators involved in the regulation cancer-associated pathways (activation and inhibition) in the HNSCC. Screening of m 1 A regulators, significantly associated with prognosis of HNSCC patients, using univariate Cox regression analysis (P < 0.05). We used Kaplan-Meier (KM) analysis and log-ranking test to investigate the association of m 1 A regulator expression with overall survival (OS) in HNSCC patients (P<0.05). OS is the time between the original diagnosis and the date of death (for any reason) [35]. Survival time and survival status information was obtained from the clinical data of the above samples. The median expression level of each m 1 A regulator gene was selected as the cut-off value, and HNSCC patients were divided into high-risk subgroups and low-risk subgroups. The R-packages "limma" were used to identify prognosis-related m 1 A regulators and drew the survival curve according to differential expression of m 1 A regulators by R package "survminer" and "survival" (P<0.05). Correlation analysis between expression of prognosis-related m 1 A regulators and infiltration levels of 24 immune cells [36] was performed using single-sample gene-set enrichment analysis (ssGSEA), which is gene set variation analysis (GSVA) package built-in algorithm [37]. In addition, we applied the R package "RColorBrewer", "psych", "igraph", and "reshape2" to build a prognostic network plot of prognosis-associated m 1 A regulators.

Consistent clustering of nine m 1 A RNA methylation regulators
Based on previously published literature associated with m 1 A methylation modification, we gathered and analyzed nine known m 1 A-regulated genes to identify different m 1 A methylation modification patterns [14]. Eight regulators were derived from merged GSE47443, GSE65858, and TCGA-HNSCC datasets for confirming distinct m 1 A modification modes mediated through m 1 A regulators. The nine m 1 A regulators, showed in Fig-ure2A, included three writers (TRMT6, TRMT61A, TRMT10C), four readers (YTHDC1, YTHDF1, YTHDF2, YTHDF3), and two erasers (ALKBH1, ALKBH3). We performed consensus clustering according to the expression of eight m 1 A RNA methylation regulators by the R-package "ConsensusClusterPlus" (50 repetitions, Pearson correlation resample rate is 0.8, cluster algorithm= PAM (Partitioning Around Medoids) to investigate diverse m 1 A modification patterns. The optimal number of clusters is selected on the basis of the area under the cumulative distribution function (CDF) curve, the correlation of the clusters, and the correlated change in the number of samples in the clusters [38].

Functional enrichment analysis by GSVA
We performed GSVA to explore the different biological processes in Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between different m 1 Aclusters by the R package "GSVA" [ (http://www.gseamsigdb.org/gsea/msigdb/index.jsp). The R packages "limma", "pheatmap", and "GSEA-Base" were utilized to merge data, perform differential analysis, clustering, and visual analysis results.

Estimation of immune cell infiltration by ssGSEA
The ssGSEA algorithm was utilized to measure the relevant abundance of each cell infiltration in the TME of HNSCC. Genomic markers for each TME-infiltrating immune cell types were obtained from previous researches [42][43][44][45]. The corresponding number of each immune cell type is represented by the enrichment fraction in the ssGSEA analysis, which is normalized to the distribution of 0-1 units. Table S2 shows the detailed genes for each TME infiltrating cell type.

Identification of differentially expressed genes (DEGs) between Function enrichment analysis and different m 1 A modification phenotypes
The previous consistent clustering algorithm has divided patients into two different m 1 A modification phenotypes, and then we identified the DEGs associated with m 1 A modifications under distinct m 1 A patterns. We used the empirical Bayesian method based on the R package "limma" to evaluate DEGs between distinct modification clusters in HNSCC cases (adjust. P-value＜0.05). We used the R package "VennDiagram" drawing the Venn diagram to indicate the number of DEGs among distinct m 1 A subtypes and determined the intersection DEGs. We performed Go functional annotation and KEGG pathway analysis based on the R packages "DOSE", "ggplot2", "clusterProfiler", "stringi", "enrichplot", and "org.Hs.eg.db" to explore the biological processes involved in the regulation of intersection DEGs. GO analysis included Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF), adjust. P-value＜0.05 was used as the cut-off.

Construction of m 1 A gene signature scoring system
The prognosis was analyzed using an unsupervised clustering method, and HNSCC patients were classified into two groups for further analysis. A consistent clustering algorithm is used to identify the number and stability of gene clusters. Subsequently, principal component analysis (PCA) was undertaken to recognize m1A-related gene features, and principal components 1 and 2 were selected as feature scores. We then used a formula similar to the one used in the previous study to ascertain the m1A score. [14,17,46,47]: m 1 A score=Σ(PC1a+PC2a) a stand for the expression of the final identified m1A phenotype-related genes.

The association between m 1 A score and tumor mutational burden (TMB)
The correlation analysis between m 1 A scores and TMB was mainly performed by the R packages "ggpubr" and "reshape2". Mutation characteristics of m 1 A-modified genes in high and low m 1 A scoring subgroups described by the waterfall function of the R package "maftools" [48]. Survival curves were plotted by applying the R package "survival" and " survminer " according to the high and low TMB status and m 1 A score.

Survival analysis and correlation analysis of m 1 A score and other elements
The median m 1 A score was selected as the cut-off value and divided the HNSCC cohort into high and low-m 1 A score subgroups. R-packages "survminer" and "survival" were applied to plot the survival curve based on the high and low-m 1 A scores. We used R package "ggalluvial", "dplyr", and "ggplot2" to analyze and visualize the association among m 1 Aclusters, gene clusters, m 1 A score, and clinical pathological characteristics (including age, survival status, gender, stage, and TNM-stage). Then, the R package "corrplot" is utilized to perform the correlation analysis between the immune cell infiltration and the m 1 A score. To further observe whether there were differences in m 1 A scores between m 1 Acluster subtypes and gene clusters, we performed a differential analysis of m1A scores. The R package "ggpubr" and "limma" were applied to analyze and visualize the results. In addition, the R package "ggplot2", "ggalluvial", and"dplyr" were applied for analyzing and visualizing the relationship between m 1 A score and clinicopathological features (including gender, stage, HPV infection status, survival status, alcohol status, radiation therapy status, and smoking status).

Correlation analysis of m 1 A score and immune checkpoint gene expression
We select immune checkpoints, including PDCD1(PD-1), CTLA-4, CD274 (PD-L1), LAG3, TIGIT, and TIM-3, for further analysis [49-51]. The R package "limma" and "ggpubr" were applied to identify the relationship between m 1 A score and immune checkpoint gene expression. We performed an immunotherapy correlation analysis further to explore the association between m 1 A score and immunotherapy. The immunotherapy score for HNSCC was obtained from The Cancer Immunome Database (TCIA, https://tcia.at/home) [52,53]. Then, we utilized the R package "ggpubr" to analyze and visualize the results through the violin map. To further explore the relationship between m 1 A score in immunotherapy, we further performed the survival analysis between different m 1 A score groups in GSE78220 cohort, a metastatic melanoma cohort treated with anti-PD-1 checkpoint inhibition therapy (pembrolizumab).

Statistical analysis
The statistical analyses in this research were developed via Rx64-4.0.3. For quantitative data, the statistical significance of normally distributed variables is accessed through Student's t-test. They were estimated through Student's t-test, while non-normally distributed variables are estimated via Wilcoxon rank-sum test. For comparison of more than two groups, the Kruskal-Wallis test and one-way analysis of Variance (ANOVA) are used as nonparametric and parametric methods. The R package "maftools" was applied to analyze and visualize Mutation Annotation Format (MAF) files of somatic mutation data [48]. All statistical P values are two-sided, with p < 0.05 as statistically significant.

The landscape of genetic variation of m 1 A regulators in HNSCC
In the present research, we explored the roles of nine m 1 A RNA methylation regulators in HNSCC. These m 1 A regulators act as dynamic reversible processes that identify, remove and add m1A modification sites and change critical biological processes (Fig-ure2A). We identified the prevailing rate of somatic mutations in nine m 1 A regulators in HNSCC. Among the 506 samples, 22 (4.35%) samples undergo genetic mutations of m 1 A modulators, including missense mutation, splice site, nonsense mutation, frame shift insert, and multi-hit. YTHDF3 indicated the highest mutation frequency, followed closely by YTHDF1, ALKBH1, ALKBH3, and YTHDF3, while other m 1 A RNA methylation regulatory genes do not demonstrate any mutations in HNSCC cases ( Figure 2B). Further analysis of the nine m 1 A regulators revealed the general existence of CNV mutations. YTHDC1, ALKBH3, TRMT61A, ALKBH1, TRMT6, and YTHDF1were mainly focused on the amplification in copy number. On the contrary, TRMT10C, YTHDF3, and YTHDF2 had an overall frequency of CNV loss ( Figure 2C). CNV alteration site of m 1 A regulators on 23 chromosomes was showed in Figure2D. The correlation analysis between nine m 1 A regulators and cancer-related pathway indicated that YTHDF3 and YTHDF2 were mainly participated in the active regulation of RTK and TSC/mTOR pathway, respectively. YTHDF1, YTHDC1, TRMT61A, TRMT6, TRMT10C, and ALKBH1 were mainly participated in the regulation of cell cycle pathway activity (Figure2E). The interaction map of m 1 A regulators and pathway showed in Figure2F. Further analysis demonstrated that m 1 A regulator expression with CNV amplification has a higher proportion of active cell cycle pathway, such as YTHDF1, YTHDC1, TRMT61A, TRMT6, and ALKBH1. While CNV-loss m 1 A regulators were significantly associated with the active RTK and TSC/mTOR pathway, such as YTHDF3 and YTHDF2 (Figure2E, 2F). These results demonstrated that the alterations of CNV might be a prominent factor contributing to the disruption of m 1 A regulators' expression. Interestingly, the expression of most of the m 1 A regulators in HNSCC cancer tissues was significantly higher than that in the corresponding normal tissues in both unpaired and paired analyses ( Figure 3A, B). The above results indicate that there are significant differences and associations between normal and HNSCC tissues at the genomic and m 1 A regulatory gene transcription levels. Therefore, it can be concluded that alternation in the expression and genetic variation of m 1 A regulators might play an essential role in regulating the progression of HNSCC.

Survival analysis of m 1 A regulators in HNSCC and immune cell correlation analysis
Univariate cox analysis indicated that only ALKBH1is the significantly high-risk prognosis-associated m 1 A regulator (Hazard ratio (HR) >1, P＜0.05). Four m 1 A regulators are significantly associated with the OS of HNSCC in the KM survival analysis. Among these m 1 A regulators, only ALKBH1 are significant high-risk m 1 A regulators for both univariate cox and KM survival analysis (all P＜0.05, HR >1). The detailed information is shown in Table 1. The high-risk m 1 A regulators, including ALKBH1, ALKBH3, TRMT6, and YTHDF2 had a detrimental effect on OS for HNSCC patients, where patients with higher expression indicated a worse survival rate than those with low expression level (all P＜0.05, Figure 3C-F). Besides, we used the Spearman method based on the ssGSEA algorithm to identify the correlation between the immune cell infiltration level and the above four m 1 A regulators with prognostic significance. We found that ALKBH1 expression is significant positive correlation with the infiltration level of Th2 cells and T helper cells; negative correlation with DC, iDC, Mast cells, Neutrophils, pDC, Treg, NK CD56dim cells, Cytotoxic cells, Th17 cells, TFH, NK cells, T cells, Eosinophils, and Macrophages (P ＜0.05, Figure 3G). ALKBH3 expression is significant negative correlation with DC, Eosinophils, iDC, Macrophages, Mast cells, Neutrophils, Tcm, TFH, Treg, Th1 cells, Tem, NK cells, T cells, T helper cells, and B cells (P＜0.05, Figure 3H). TRMT6 expression is significant positive correlation with the infiltration level of Th2 cells; negative correlation with Cytotoxic cells, DC, iDC, Neutrophils, NK CD56dim cells, pDC, T cells, TFH, Th17 cells, Treg, Th1 cells, Macrophages, NK cells, aDC, Mast cells, and Tem (P＜0.05, Figure 3I). YTHDF2 expression is significant positive correlation with the infiltration level of T helper cells, Th2 cells, Tcm, aDC, and Th17 cells; negative correlation with iDC, Macrophages, pDC, DC, Neutrophils, and Mast cells (P＜0.05, Figure 3J). The detailed immune cell correlation analysis results listed in Table S3-S6.

Consistent clustering of eight m 1 A RNA methylation regulators
The m 1 A regulatory network demonstrated the interaction of eight m 1 A regulators, the connections between the regulators, and their prognostic significance in HNSCC patients. Except for YTHDC1 has a negative correlation with ALKBH3 and TRMT61A. The other m 1 A regulators were positively correlated with each other ( Figure 4A). The above results suggest that the linkage of these interactions among the regulators of erasers, readers, and writers might play essential roles in the composition of multiple m 1 A modification modes and are associate with cancer occurrence and progression. Based on this assumption, we divided the samples into different m 1 A modification patterns according to the expression of the eight m 1 A regulators using consensus clustering analysis. Therefore, we identified two different clusters of modification patterns based on merge the GSE47443 and GSE65858 cohorts, including 405 patients in m 1 A cluster A and 364 patients in cluster B. m 1 Acluster-A demonstrated remarkable survival advantage. In contrast, m 1 Acluster-B had the worst OS (P＜0.05, Figure 4B, FigureS1). We also observed an extraordinary distinction in the m 1 A regulator expression between diverse m 1 A modification patterns and clinicopathological features of HNSCC. The heatmap indicated that the high-risk m 1 A regulators that we previously screened for survival significance were mainly highly expressed in m 1 Acluster-B subtypes, such as ALKBH1, ALKBH3, TRMT6, and YTHDF2, as the significant high-risk m 1 A regulators with a detrimental effect on OS for HNSCC patients, significantly elevated their expression in the m 1 Acluster-B subtype (Figure4C).

Functional enrichment analysis by Gene set variation analysis (GSVA)
We performed a GSVA analysis to determine the biological processes in these distinct m 1 A modification patterns (Figure5). The m 1 Acluster-B was mainly enriched in mRNA methylation, modification, and RNA polymerase-associated biological processes (Fig-ure5A). We also performed GSVA enrichment analysis to identify KEGG pathways analysis among these different m 1 A modification patterns (Figure5B). The m 1 Acluster A was mainly involved in the regulation pathway of metabolism-related, intestinal immune network for IGA production, neuroactive ligand receptor interaction. etc.m 1 Acluster B prominently participated in the regulation pathway of ubiquitin mediated proteolysis, oocyte meiosis, basal transcription factors, nucleotide excision repair, RNA degradation, etc. (Fig-ure5B).

Characterization of TME cell infiltration under different m 1 A modification patterns
The ssGSEA-based analysis of TME cell infiltration by different m 1 A modification patterns indicated that the two m 1 A modification patterns had significantly diverse TME cell infiltration characteristics. m 1 Acluster A was categorized as immune-inflamed phenotype, typified through adaptive immune activation and immune cell infiltration. m 1 A cluster A had the highest number of immune cell species with high immune cell infiltration levels, including Activated B-cell, Activated CD4 T-cell, Activated CD8 T-cell, Activated dendritic cell, Eosinophil, Immature B-cell, Immature dendritic cell, MDSC, Macrophage, Mast cell, Monocyte, Natural killer T-cell, Natural killer cell, Plasmacytoid dendritic cell, Regulatory T-cell, T follicular helper cell, Type.1 T helper cell, and Type.17 T helper cell ( Figure 5C). m 1 Acluster B was classified as an immuno-desert phenotype, characterized through immunosuppression. The cell content of neutrophil and Type.2 T helper cell is the highest in the m 1 Acluster B ( Figure 5C).

The DEGs between distinct m 1 A modification phenotypes and function enrichment analysis
Although HNSCC patients have been classified into two m 1 A clusters based on a consistent clustering algorithm for m 1 A-regulated gene expression, the underlying genetic changes and expression perturbations in these clusters are still poorly understood. Based on the above questions, we further explored the changes in potential m 1 A-related transcript expression of the two m 1 A clusters in HNSCC. 131 genes were selected as the m 1 Aassociated gene signature of the two m 1 Aclusters. PCA analysis indicated significant diversities in the m 1 A transcriptional profile among the two distinct m 1 Aclusters ( Figure  5D). GO enrichment analysis of these genes signatures was significantly over-expressed in the BP associated with regulating neutrophil activation involved in immune response, neutrophil degranulation, T cell activation, etc. CC is related to cell−substrate junction, focal adhesion, vacuolar membrane, et al., and MF is associated with ATPase activity, cadherin binding, and DNA−binding transcription factor binding, et al. The top-10 Go terms were visualized by Figure5E. The detailed GO analysis results showed in Table S7.  Table S8.
Based on the 131 most typical m 1 A phenotype-related genes, we performed an unmonitored consensus clustering analysis and obtained two stable gene clusters (FigureS2).
These categorizations divided HNSCC patients into two distinct m 1 A gene clusters with diverse clinicopathologic characteristics and were named geneCluster A and geneCluster B ( Figure 6A). This suggests that two different m 1 A methylation modification modes indeed exist in HNSCC. We observed that these m 1 A phenotype-associated signature genes are primarily lower expressed in geneCluster A, and higher expressed in gene cluster B ( Figure 6A). In addition, we also compared the expression levels of eight m 1 A regulators in the two gene clusters. We observed a significant distinction in the expression of eight m 1 A regulators between the two gene clusters, which is consistent with the anticipated results of the m 1 A methylation modification modes. The expression of TRMT6, TRMT61A, YTHDC1, YTHDF1, YTHDF2, YTHDF3, ALKBH1, and ALKBH3 in geneCluster A is higher than that in geneCluster B ( Figure 6B). Further survival analysis demonstrated a remarkable prognostic distinct among the two m 1 A gene clusters in HNSCC cases. 497 HNSCC patients have clustered in geneCluster A signature, which is proven to be associated with poor OS, whereas HNSCC patients in geneCluster B (272 patients) have the better OS (Figure6C).

Exploration of the generation of m 1 A gene signature score and its clinical significance
Although our findings confirm the impact of m 1 A methylation modifications on prognosis and immune cell infiltration, these analyses focus on overall HNSCC patients and do not accurately predict m 1 A methylation modification patterns in individual tumor patients. Therefore, we developed a scoring system named m 1 A score to assess the pattern of m 1 A modification in individual HNSCC patients based on the identified m 1 A-related genes. We examined the correlation between immune cell infiltration and m 1 A scores to characterize the m 1 A signatures better. The m 1 Ascore significantly positively correlated with the infiltration level of Activated B cell, Activated CD4 T-cell, Activated CD8 T-cell, Activated dendritic cell, CD56bright tural killer cell, Eosinophil cell, Gamma delta T cell, Immature B cell, Immature dendritic cell, MDSC, Macrophage, Mast cell, Monocyte, Natural killer T cell, Natural killer cell, Plasmacytoid dendritic cell, Regulatory T cell, T follicular helper cell, Type 1 T helper cell, and Type.17 T helper cell. m 1 Ascore was significantly negatively associated with the infiltration level of CD56dim tural killer cell and Neutrophil (Figure6D). Considering the complexity of the m 1 A modification assessment, we visualized the workflow of m 1 A score construction using Sankey diagrams (Figure6E). Kruskal-Wallis test found significantly distinct m 1 A scores between gene clusters. Gene-Cluster A indicated the lowest median m 1 Ascore, while geneCluster B had the highest median m 1 Ascore, which demonstrated high m 1 Ascore might be strongly related to immune activation-associated signatures, while low m 1 Ascore may be associated with stromal activation linked signature (Figure6F). More significantly, m 1 Acluster A had significantly higher m 1 A scores than m 1 Acluster B ( Figure 6G). The above results revealed that high m 1 A scores were significantly associated with immune cell activation, and low m 1 A scores were associated with stroma activation. The m 1 A score provides a better evaluation of the m 1 A modification pattern of individual tumors and further assesses tumor TME cell infiltration characteristics to identify false and true TME immune cell infiltration.

Survival analysis and correlation analysis between m 1 A score and clinical features
We also examined the value of m 1 Ascore in forecasting OS in HNSCC patients. The R package "Survminer" was used to determine the m 1 A cut-off value of -2.071, thus dividing the 769 HNSCC patients into high and low m 1 Ascore subgroups (Table S9). Patients with a high m 1 Ascore indicated a conspicuous survival benefit (P＜0.001, Figure 7A). The 3-year survival rate was higher than that of patients with a low m 1 Ascore. Univariate Cox regression analysis showed that age, M stage, and m 1 Ascore are associated with the prognosis (all P＜0.05, Figure 7B). Multivariate Cox analysis also identified that m 1 Ascore is a stable and independent prognostic marker for evaluating patient outcome (P＜0.05, HR: 0.961 (0.936-0.987), Figure7C). In addition, we explored the relationship between m 1 Ascore and the clinical characteristics of HNSCC patients. People in the different ages, male, T3-4 patients, advanced stage (stage III＆IV), N0, N1-3, and M0 patients, with high m 1 Ascore had the better OS than the patient with low m 1 Ascore (all P＜0.05, Figure 7D-K). The above results indicated that m 1 Ascore could also be utilized to assess certain clinical features of HNSCC patients, such as clinical stage, etc.

Relationship between m 1 A score and tumor somatic mutation
There is growing evidence of the relevance of TMB to immunotherapy. High TMB was associated with sustained clinical response to anti-PD-1/PD-L1 immunotherapy [54]. Thus, we analyzed the discrepancies in the distribution of somatic mutation between high and low m 1 Ascore groups in the TCGA-HNSCC cohort through the R package "maftools". TMB mutation frequency was lower in the high-m 1 Ascore group than in the lowm 1 Ascore group ( Figure 8A, B). The m 1 Ascore and TMB demonstrated a negative correlation in different gene clusters ( Figure 8C). The TMB quantification analysis demonstrated that the low m 1 Ascore tumors were related to a higher TMB ( Figure 8D). Survival analysis indicated that patients with higher TMB have a prominently lower OS than those with lower TMB patients ( Figure 8E). We found that HNSCC patients with merged low TMB and high m 1 Ascore had the best OS, while HNSCC patients with a combination of high TMB and low m 1 Ascore had the worst OS ( Figure 8F). Thus, the above finding indirectly indicated that tumor m 1 A modification patterns might be an essential element that mediated clinical reaction to anti-PD-1/PD-L1 immunotherapy. Moreover, the value of the m 1 Ascore in forecasting the outcome of immunotherapy has been indirectly confirmed. These results will also provide new insights to investigate the mechanisms of m 1 A methylation modifications in tumor somatic mutations, the formation of TME, and the role of ICIs therapy.

The role of the m 1 A score in predicting the effect of immunotherapy
Immune checkpoint inhibitors such as anti-PD-L1, PD-1, and CTLA-4 immunotherapy have certainly made considerable breakthroughs in cancer treatment. Novel immune checkpoints with potential therapeutic significance were also identified, including LAG3, CD80, CD86, TNFRSF9, and HAVCR2. Because of the strong correlation between m 1 A score and immune response, we further investigated the relationship between m 1 A score and immune checkpoint gene expression. The expression levels of PD-1, PD-L1, CTLA-4, LAG3, CD80, CD86, TNFRSF9, and HAVCR2, in the high m 1 A score subgroup were significantly higher than that in the low m 1 A score subgroup, respectively ( Figure 8G-N). Besides, Immunophenoscore (IPS) was used to evaluate response to ICIs in HNSCC patients. IPS with CTLA_ negative_ PD1_ negative was significantly increased in the low m 1 A score group ( Figure 8O). IPS with CTLA_ positive_ PD1_ positive was significantly increased in the high m 1 A score group ( Figure 8P). The detailed IPS data of HNSCC patients showed in Table S9. These results indicated that anti-CTLA-4＆anti-PD1 treatment in the high m 1 A score subgroup is better than that in the low m 1 A score subgroup. Furthermore, we successfully validated the role of the m 1 A score in predicting the response to ICB in GSE78220 cohort (melanoma). In the anti-PD-1 cohort, the low m 1 Ascore group has presented a significantly prolonged OS ( Figure 8Q). In summary, our study presents robust evidence that m 1 A methylation modification patterns significantly correlate with tumor immunophenotype and immunotherapy response and that the developed m 1 A score contributes to the prediction of response to immunotherapy. It can further predict the prognosis of HNSCC patients.

Discussion
There has been report that m 1 A modification regulators were correlated with overall cancer survival [55]. Since most studies have centered on a single TME cell type or a single m 1 A regulator, the full signature of TME permeation regulation by the combined action of multiple m 1 A regulators has not been fully characterized [56]. Therefore, understanding the effect of distinct m 1 A modification profiles in the cellular infiltration of TME will contribute to a deeper comprehension of the interplay of m 1 A RNA methylation in the antitumor immune response and facilitate more effective and precise immunotherapeutic approaches. In this present research, for the first time, we have systematically explored the characteristics of TME cell infiltration in HNSCC mediated by the integrative effects of multiple m 1 A regulators and corresponding modification patterns. According to the TCGA-HNSCC, GSE47443, and GSE65858 cohorts, we identified two diverse m 1 A methylation modification patterns with different immunophenotypes and are associated with specific anti-cancer immunity. In addition, we developed a scoring system called m 1 Ascore to assess the m 1 A modification pattern of individual HNSCC patients.
We first performed a comprehensive analysis somatic mutation of m 1 A modulators in HNSCC and found that 22 of 506 patients undergo genetic mutations. YTHDF3 indicated the highest mutation frequency, followed closely by YTHDF1, ALKBH1, ALKBH3, and YTHDF3. Further analysis of the nine m 1 A regulators revealed the general existence of CNV mutations. Six m 1 A regulators were mainly focused on the amplification in copy number. Alterations of m 1 A regulators may lead to dysregulation of m 1 A-associated genes, thus demonstrating their involvement in tumorigenesis. Yeon et al. reported that TRMT6 often perturbed in cancer [57]. Further analysis demonstrated that m 1 A regulator expression with CNV amplification has a higher proportion of active cell cycle pathway, such as YTHDF1, YTHDC1, TRMT61A, TRMT6, and ALKBH1. Chi et al. found that cell cycle pathway was enriched in the group with high expression of YTHDF1 [58]. TRMT6 was also found to be enriched in the cell cycle pathway [59]. While CNV-loss m 1 A regulators were significantly associated with the active RTK and TSC/mTOR pathway, such as YTHDF3 and YTHDF2. In addition, the relationship between YTHDC1, TRMT61A, ALKBH1, YTHDF3, YTHDF2 and the cancer-related pathway has not been reported previously and we are reporting it for the first time. Interestingly, the expression of most of the m 1 A regulators in HNSCC cancer tissues was significantly higher than that in the corresponding normal tissues. Previous reports suggest that dysregulation of m 1 A regulators is associated with gastrointestinal tumorigenesis and progression of hepatocellular carcinoma [14,60,61]. Our research identified that the expression profiles between m 1 A "erasers" and "writers" are unbalanced. In theory, these unbalanced expression profiles might lead to abnormal m 1 A modification patterns, which result in the development of HNSCC. Besides, ALKBH1, ALKBH3, TRMT6, and YTHDF2 had a detrimental effect on OS for HNSCC patients, which was consistent with previous studies [62,63]. Li et al. reported that m 1 A regulators are associated with the regulating of carcinogenic pathways and OS [55]. m 1 A regulators also associated with HCC patient prognosis [61]. However, little is known about the correlation between m 1 A regulators and TME infiltration. Therefore, we found a strong correlation between the four identified m 1 A regulators, associated with the OS of HNSCC, and the current TME immune infiltrating cells, based on the ssGSEA algorithm for immune infiltrating cells, using the Spearman method. The above results indicate a potential role of m 1 A regulators in the TME and progression of HNSCC.
More significantly, these m 1 A regulators are interrelated and form a tight network of interactions. These discoveries have motivated us to construct a comprehensive clustering analysis rather than analyzing the role of individual m 1 A regulators. Therefore, a systematic analysis of all m 1 A regulators can provide a more comprehensive view of m 1 A modification patterns in TME. Previous studies on hepatocellular carcinoma and colon cancers have demonstrated that m 1 A regulator-based classification models have good typing ability, have significant differences in survival between phenotypes, and are involved in regulating distinct pathways. Diverse phenotypes have significantly distinct characteristics of TME cell infiltration [14,61]. Therefore, in this present research, based on eight m 1 A regulators, we determine two different m 1 A methylation modification modes that have different OS, immunophenotypes and are associated with diverse anti-cancer immunity. We found that m 1 Acluster A demonstrated prominent survival advantage, while m 1 Acluster B had the poor OS. Further research found that high-risk m 1  . Our analyses revealed that the m 1 A regulators were significantly highly expressed in m 1 Acluster B, and m 1 Acluster B had the worse prognosis between the two m 1 Aclusters, which indirectly indicated that our classification was very accurate.
Furthermore, we also discovered m 1 Acluster B was distinctly enriched in RNA methylation and carcinogenic activation biological process or pathways, such as RNA methylation, DNA templated transcription termination, regulation of translational initiation, ubiquitin mediated proteolysis, nucleotide excision repair, cell adhesion molecules cams, et al. RNA methylation plays a critical role in tumorigenesis [66]. According to the ssGSEA, two m 1 A modification modes had remarkably different TME cell infiltration features. m 1 ACluster A was categorized as immunoinflammatory phenotype, featured through adaptive immune activation and immune cell infiltration. m 1 Acluster B was classified as an immuno-desert phenotype, characterized through immunosuppression. The immunoinflammatory phenotype, tumors are highly inflamed with CD8+ T cells ('infiltrated' or 'hot' tumors), show massive immune cell infiltration in TME [27,72]. Increased CD8+ T-cell infiltration after relapse is a good prognostic indicator for HNSCC [73]. The m 1 A cluster A has massive immune cell infiltration in TME and the better prognosis after thoroughly investigating the features of TME cell infiltration induced through distinct m 1 A modification patterns, which is consistent with previous research in hepatocellular carcinoma and colon cancer [14,61]. In combination with each cluster's TME cell infiltration features, this identified the credibility of our immunophenotypic categorization of the distinct m 1 A modification patterns.
In addition, a total of 131 DEGs were identified as the m 1 A-associated gene signature from distinct m 1 A modification patterns. Similar to the m 1 A modification clustering results, two genomic clusters according to the m 1 A signature gene were confirmed, significantly correlated with distinct clinicopathologic characteristics, survival outcomes, and TME landscapes. GeneCluster A was related to worse OS and geneCluster B was associated with the better OS. Meanwhile, the higher expression of TRMT6, TRMT61A, YTHDC1, YTHDF1, YTHDF2, YTHDF3, ALKBH1, ALKBH3 was observed in the gene-Cluster A. This again demonstrates the importance of m 1 A modification in shaping diverse TME landscapes. Therefore, a comprehensive assessment of m 1 A modification modalities will enhance our comprehension of TME cell infiltration characteristics. Because of the individual heterogeneity of m 1 A modifications, there is an urgent need to quantify the m 1 A modification patterns of individual tumors. For this purpose, we developed a protocol named m 1 Ascore to assess the m 1 A modification pattern of individual HNSCC patients. The m 1 A modification pattern characterized by the immune-desert phenotype exhibited lower m 1 Ascores, whereas the immune-inflammatory phenotype exhibited higher m 1 Ascores. We also found m 1 Ascore significantly positively correlated with the infiltration level of many immune cells. Patients with a high m 1 Ascore indicated a remarkable survival benefit. These findings provide further evidence that HNSCC patients with high levels of immune cell infiltration have a higher OS. In HNSCC, increased infiltration of CD8+ T cells was the only immune cell type associated with improved survival, regardless of tumor location, stage, and treatment [74,75]. The m 1 Ascore might also be utilized to assess certain clinical features of HNSCC patients, such as clinical stage et al. Further analysis also confirmed that m 1 Ascore is an independent prognostic marker for HNSCC. These results demonstrate that the m 1 Ascore is a credible and stable tool for synthetically assessing individual tumor m 1 Ascore for further recognition of TME infiltration patterns, and our m 1 Ascore indicates a predictive benefit specific immunotherapy for HNSCC.
There is growing evidence that TMB is a predictor of solid cancer patients receiving anti-PD-1/CTLA-4 combination immunotherapy, and high TMB status (≥10 mutations per megabase (MB) genome) indicated a durable clinical reaction to anti-PD-1/PD-L1 immunotherapy, such as pembrolizumab (Keytruda) and nivolumab (Opdivo), and has more prolonged survival in cancers such as melanoma and lung cancer [54,76,77]. Amazingly, unlike other previously reported tumors, patients with HNSCC with high TMB expression had a poorer prognosis. Our research revealed a negative correlation between m 1 Ascore and TMB in different gene clusters. The low m 1 Ascore subgroup indicated more extensive TMB than the high m 1 Ascore subgroup. Patients in the high TMB subgroup had a significantly lower OS than those in the low TMB subgroup. Meanwhile, HNSCC patients with a combination of low TMB and high m 1 Ascore had the best OS. Zhang et al. also reported that HNSCC patients with high TMB have worse prognosis than those with low TMB, and TMB might affect CD4+T cell and B cell infiltration [78]. Braun et al. has reported that tumors are heavily infiltrated with immune CD8 T cells in melanoma and some other cancers, creating a so-called inflammatory or "hot" environment within the tumor that responds better to PD-1 blockade. HNSCC tumors are highly immuno-infiltrative, but have an overall immunosuppressive TME profile [79].
Considering the close association between the m 1 A score and immune response, we further investigated the connection of the m 1 A score with the expression of immune checkpoints. We found that the expression level of PD-1, PD-L1, CTLA-4, LAG3, CD80, CD86, TNFRSF9, and HAVCR2 are significantly higher in the high m 1 A score group. Since the high m 1 A score group has a better prognosis, it also means that the high m 1 A score group can benefit from immunotherapy, which further confirms that immunotherapy can prolong the overall survival of HNSCC patients. These findings also reveal new therapeutic targets (LAG3, CD80, CD86, TNFRSF9, and HAVCR2) beyond PD-1, PD-L1, and CTLA-4 and new combinatorial approaches to further understand immune checkpoint biology in the future. We also found that the effect of anti-CTLA-4＆anti-PD1 treatment in the high m 1 A score subgroup is better than that in the low m 1 A score subgroup. Wang et al. reported >70% of HNSCC lesions responded to intratumoral anti-CTLA-4 [80]. Checkpoint inhibitors, such as anti-PD-1 and anti-PD-L1 antibodies, have been shown to significantly improve overall survival and disease-free survival of HNSCC after the failure of platinum-based chemotherapy [21,81].
Although we reviewed the literature screening, nine recognized m 1 A RNA methylation regulators, comprehensively assessed the m 1 A modification patterns in HNSCC, and systematically linked these modification patterns to the permeability properties of TME cells. However, there are also few deficiencies to this study. First, a novel array of identified regulators needs to be integrated into the model to optimize the accuracy of m 1 A modifications to the model in the future. Second, m 1 A modification patterns and m 1 A scores were verified using a retrospective dataset; therefore, a prospective cohort study of HNSCC patients receiving immunotherapy is necessary to validate our results in future work.
In conclusion, m 1 Ascore can be used to comprehensively evaluate the m 1 A methylation modification patterns of individual HNSCC patients and their corresponding TME cell infiltration characteristics, further clarify the immunophenotype of the tumor and guide more rational clinical treatment. We have also identified that the m 1 A score can be applied to evaluate HNSCC patients' clinicopathological features, immune cell infiltration, and tumor mutation burden. Similarly, m 1 Ascore can be considered as an independent prognostic factor to predict the survival of HNSCC patients. We can also predict the clinical response to anti-PD-1/PD-L1or CTLA-4 immunotherapy by m 1 Ascore. What's more, this research has generated some innovative insights into cancer immunotherapy targeted m 1 A modulators or m 1 A phenotype-associated genes to alter m 1 A modification modes, further inverting poor TME cell infiltration signature, conversion of "hot tumors" into "cold tumors", might facilitate to exploring the improvement of new drug, combination approaches or new immunotherapeutic agents in the future. Our study results offered innovative insights into improving the clinical reaction of HNSCC patients to immunotherapy, identified distinct tumor immune phenotypes, and facilitated personalized cancer immunotherapy in the future.

Conclusions
This study confirms the widespread regulatory mechanism of m 1 A methylation modification on the TME. The diversities in m 1 A modification patterns are a non-negligible element contributing to the heterogeneity and sophistication of the individual TME. A comprehensive evaluation of m 1 A modification patterns in individual HNSCC tumors will facilitate in-depth comprehending of the TME, immune cell infiltration features and guide more efficient immunotherapeutic approaches.  Table S1. The baseline information of TCGA-HNSCC, GSE47443, and GSE65858 cohorts. Table S2. The detailed genomes that mark each TME-infiltrating cell type. Table S3. The correlation analysis between ALKBH1 expression and immune cell. Table S4. The correlation analysis between ALKBH3 expression and immune cell. Table S5. The correlation analysis between TRMT6 expression and immune cell. Table S6. The correlation analysis between YTHDF2 expression and immune cell.