M6A-mediated molecular patterns and tumor microenvironment infiltration characterization in nasopharyngeal carcinoma

ABSTRACT N6-methyladenosine (m6A) is the most predominant RNA epigenetic regulation in eukaryotic cells. Numerous evidence revealed that m6A modification exerts a crucial role in the regulation of tumor microenvironment (TME) cell infiltration in several tumors. Nevertheless, the potential role and mechanism of m6A modification in nasopharyngeal carcinoma (NPC) remains unknown. mRNA expression data and clinical information from GSE102349, and GSE53819 datasets obtained from Gene Expression Omnibus (GEO) was used for differential gene expression and subsequent analysis. Consensus clustering was used to identify m6A-related molecular patterns of 88 NPC samples based on prognostic m6A regulators using Univariate Cox analysis. The TME cell-infiltrating characteristics of each m6A-related subclass were explored using single-sample gene set enrichment (ssGSEA) algorithm and CIBERSORT algotithm. DEGs between two m6A-related subclasses were screened using edgeR package. The prognostic signature and predicated nomogram were constructed based on the m6A-related DEGs. The cell infiltration and expression of prognostic signature in NPC was determined using immunohistochemistry (IHC) analysis. Chi-square test was used to analysis the significance of difference of the categorical variables. And survival analysis was performed using Kaplan–Meier plots and log-rank tests. The NPC samples were divided into two m6A-related subclasses. The TME cell-infiltrating characteristics analyses indicated that cluster 1 is characterized by immune-related and metabolism pathways activation, better response to anit-PD1 and anti-CTLA4 treatment and chemotherapy. And cluster 2 is characterized by stromal activation, low expression of HLA family and immune checkpoints, and a worse response to anti-PD1 and anti-CTLA4 treatment and chemotherapy. Furthermore, we identified 1558 DEGs between two m6A-related subclasses and constructed prognostic signatures to predicate the progression-free survival (PFS) for NPC patients. Compared to non-tumor samples, REEP2, TMSB15A, DSEL, and ID4 were upregulated in NPC samples. High expression of REEP2 and TMSB15A showed poor survival in NPC patients. The interaction between REEP2, TMSB15A, DSEL, ID4, and m6A regulators was detected. Our finding indicated that m6A modification plays an important role in the regulation of TME heterogeneity and complexity.


Introduction
Nasopharyngeal carcinoma (NPC) is a subtype of head and neck tumor, which is an epithelial carcinoma that originates from nasopharyngeal mucosal lining and is different from other head and neck tumors. 1 Based on the morphologic characteristics, NPC is divided into epithelial carcinoma (EC), sarcomatoid carcinoma (SC), mixed sarcomatoidepithelial carcinoma (MSEC), and squamous cell carcinoma (SCC). 2 NPC is a rare tumor comparison with other tumors worldwide, which is 129,079 new cases and 73,987 deaths in 2018. 34][5] The geographical distribution of NPC is also extremely unbalanced in China, higher in South China than in other regions. 6Even with advanced regimens, such as liquid biopsies, minimally invasive surgery, radiotherapy, chemotherapy, and immunotherapy widely used for NPC treatment, [7][8][9] approximately 10-20% of NPC patients suffer recurrence after primary treatment. 10,11Also, the asymptomatic nature of NPC hinders the early diagnose of NPC.NPC is a heterogeneous tumor with complexity contributed to the risk factors including genomic variations, Epstein-Barr virus (EBV) infection, smoking, diet, gender, environmental, and other clinical characteristics. 12,13Therefore, finding the novel marker is of important in the early diagnose and treatment of NPC.
Tumor microenvironment (TME) consists of cancer cells, variation immune cells, stromal cells, bone marrow-derived cells, and they secreted cytokine, chemokines, and growth factors, 14 and is characterized by hypoxia, high oxidation, acidity, malnutrition, and inflammation because of the heterogenous of cells in TME. 15,16TME associates to hall makers of cancer, 17 radio-/chemoresistance, 18,19 immune evasion and immunotherapeutic response, 20 tumor progression and recurrence, 21 and prognosis. 22When it refers to NPC, NPC has the most severe stromal infiltration compared to other solid tumors, and the pathological status might alter cellular composition in the NPC microenvironment. 23Liu et al. 24 found that in EBV-induced NPC, the cell infiltration of CD8 + T cells cells and NK cells were limited.In contrast, regulatory T cells (Treg), M2 macrophage, and B cells can promote tumor cell proliferation by inhibiting the activity of CD8+ T cells and promote metastasis, while tumor endothelial cells (TEC) and cancer-associated fibroblasts (CAF) were increased.TME exert a important role in tumor growth and development and its composition needs to be investigated in more depth.
In the present study, we systematically collected mRNA and corresponding clinical information of NPC patients to identify the m6A-related molecular patterns in NPC explored the TME cell infiltrating characteristics, and predicated the therapeutic response for each subclass.Moreover, we also constructed an m6A-related prognostic risk model and identify the prognostic signature, then developed a predicated nomogram to predicate progression-free survival (PFS) for NPC patients.Our finding might provide the potential biomarkers for NPC diagnosis and treatment, and supply the evidence for managing the therapeutic strategy.

Data downloading and processing
mRNA expression data of NPC patients of two datasets, GSE102349, and GSE53819 datasets, were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/).The mRNA expression data of the GSE102349 dataset (Submission date: Aug 08, 2017, Last update date: May 15, 2019, data retrieval: Jan 08, 2022) 42 was generated using Illumina HiSeq 2000, this dataset includes mRNA expression data and clinical information from 113 NPC patients, and 88 cases with progression-free survival (PFS) involved in this study.And the mRNA expression data of the GSE53819 dataset (Submission date: Jan 04, 2014, Last update date: Aug 01, 2019, data retrieval: Jan 08, 2022) 43 was generated using Agilent -014,850 Whole Human Genome Microarray 4 × 44K G4112F, and this dataset contained 18 NPC primary tumor tissues and 18 non-cancerous nasopharyngeal tissues.

Identification and validation of the molecular subclusters
The unsupervised class was a powerful technique in cancer research, and the consensus clustering methods were used for estimating the number of unsupervised classes in a dataset. 44he consensus clustering was performed using the consensusClusterPlus package in R to conduct the classification in the GSE102349 dataset (88 NPC patients) according to their expression profiles of prognosis-related m6A regulators.The k value was determined with consistent cumulative distribution function (CDF) reached an approximate maximum, that represented the optimal number of clusters.T-Distributed Stochastic Neighbor Embedding (t-SNE) is a non-linear dimensionality reduction method with minimum structural information loss and is widely used for bioinformatics. 45In this study, t-SNE was performed to validate the subtype assignments using mRNA expression data of GSE102349 and GSE53819 datasets.Moreover, a Submap matrix was conducted to detect the similarity of molecular classes from GSE102349 and GSE53819 datasets. 46In addition, differential expression of the 10 m6A regulators between subclusters was detected and visualized using pheatmap package in R.

Gene set variation analysis (GSVA) for function enrichment
GSVA is a Gene set enrichment (GSE) method that estimates variation pathway and biological processes activity over a sample population in a non-parametric and unsupervised method. 47Here, we conducted the GSVA using GSVA package in R to investigate the various biological function among distinct m6A subclasses.The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene set (c2.cp.kegg.v7.4.symbols.gmt) and 50 human cancer gene sets (h.all.v7.4.symbols.gmt)were downloaded from Molecular Signatures Database (MSigDB, data retrieval: Jan 08, 2022, https://www.gseamsigdb.org/gsea/msigdb/index.jsp), Biological process gene set (c5.go.bp.v7.4.symbols.gmt)was downloaded from Gene Ontology (GO) database (data retrieval: Jan 08, 2022, http:// geneontology.org/),and pathway gene set (c2.cp.reactome.v7.4.symbols.gmt)was obtained from Reactome pathway database (data retrieval: Jan 08, 2022, https://reactome.org/).Besides, we also obtained the gene sets that related to CD8 T effector, DNA damage repair, Nucleotide excision repair from MSigDB, and a previous study. 48Above gene sets were performed GSVA to identify the enriched biological processes and pathways.Finally, the various biological processes and pathways between distinct subclusters were determined by Wilcox rank-sum or Kruskal-Wallis tests.Adjusted P-value <.05 was considered statistically significant.

Estimate tumor immune cell infiltration
We used the ESTIMATE algorithm to estimate the composition of the tumor microenvironment (TME). 49The ratio of the immune cells and stromal cells in the TME was quantized with the immune score and stromal score.The differences between distinct subclasses were detected by Wilcox rank-sum or Kruskal-Wallis tests.In addition, we used the single-sample gene set enrichment (ssGSEA) algorithm 50 to estimate TME immune cell infiltration according to the marker genes of 24 immune cell types. 51And the relative abundance of each cell infiltration was quantized by conducting ssGSEA using GSVA package in R. the variation between m6A-related subclusters were examined by Wilcox rank-sum or Kruskal-Wallis tests.Also, gene expression data from GSE102349 was used to assess the relative abundance of 22 immune infiltration in two subclasses.Also, the the difference in immune cell infiltration was tested using the Wilcoxon test, and the difference was considered significant when p < .05.

Evaluation of the response for immunotherapy of each subclass
Human leukocyte antigen (HLA) genes and immune checkpoints exert a key role in response for immunotherapy. in this study, we detected the expression of HLA genes (HLA-A/B/C/ E/F/G/DRB1/DQB1/DQA1/DPB1/DRA/DRB5/DPA1/DMA/ DMB/DQA2/DOA/DQB2/DOB) and immune checkpoints (IDO1, PD-L1 (CD274), PD-L2 (PDCD1LG2), TIM-3 (HAVCR2), PD-1 (PDCD1), CTLA-4, CD80, CD86) between m6A-related subclasses.Furthermore, tumor immune dysfunction and exclusion (TIDE) and the SubMap algorithm were used to predict the immune response to checkpoint inhibitors PD-1 and CTLA4 in distinct m6A-related subclasses.TIDE is a computational method to predicate the immune checkpoint blocked (ICB) response based on two primary mechanisms of tumor immune evasion, including the induction of T cell dysfunction in tumors with high cytotoxic T lymphocytes (CTL) levels and the prevention of T cell infiltration in tumor with low CTL level. 52The TIDE score reflects the efficiency of the immunotherapeutic outcome, high TIDE score indicates a poor treatment response rate. 53oreover, the Submap mapping was used to compare the response for anti-PD-1 and anti-CTLA4 treatment between m6A-related subclusters.P-value <.05 and Bonferroni corrected P-value <.05 considered statistically significant.

Evaluation of the response for chemotherapy of each subclass
We also investigated the chemotherapy response rate using pRRophetic package in R. According to the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/), 53pRRophetic algorithm was conducted to identify the antitumor drugs from 138 common antitumor drugs using Ridge's regression with the half-maximum inhibitory concentration (IC50) of each patient.The low IC50 value indicated a high response for chemotherapy.

Characterization of m6A-based subclasses in NPC
The differentially expressed genes (DEGs) among NPC subclasses were identified using edgeR package with parameters of |log2 (fold change, FC) |>1 and FDR < 0.05.ggplot2 and pheatmap R packages were used to visualize the DEGs among NPC subclasses.

Biofunction enrichment of DEGs of m6A-based subclass
GO and KEGG pathway enrichment was conducted using clusterProfiler package in R with P-value <.05 or q-value <0.05.GO annotation included biological process (BP), molecular function (MF), and cellular components (CC) terms.

Construction of the m6A-related gene signature and risk model
88 samples in the GSE102394 dataset were divided into a training set and test set at a 7:3 ratio.The expression data of DEGs among m6A-related subclasses and their corresponding clinical information were used for univariate Cox and Kaplan-Meier analyses with P-value <.05.Then, the m6Arelated genes related to prognosis were identified using LASSO regression analysis using glmet package in R. The penalty parameter (λ) of the model was confirmed by tenfold cross-validation.To estimate the prognostic value of the above prognostic signature, the risk score of the training set, test set, and whole set was calculated according to the formula:

Development of a nomogram for prognosis predicating
A nomogram predicted model was constructed based on the risk factors using rms package in R. and the prediction accuracy of the nomogram was estimated using the calibration curve.Moreover, we also investigated the correlation between prognostic factors and m6A regulators using Pearson correlation analysis.And P-value < 5 and correlation coefficient > 0.3 considered significant correlation.

Gene Set Enrichment Analysis (GSEA)
Single gene GSEA was conducted in accordance with KEGG (c2.cp.kegg.v2022.1.Hs.entrez.gmt)gene sets in GSEA database (http://www.gsea-msigdb.org/gsea/msigdb,data retrieval: Nov 10, 2023) using clusterProfilerj and org.Hs.eg.db package in R script.The sample was classified as high-and low-gene expressed group referred to the expression of REEP2, TMSB15A, DSEL, and ID4, and spearman correlation was calculated between the expression of m6A-related prognostic signature and other genes in the dataset.The pathway was selected when padjust < 0.05 and |NES| > 1, and TOP 10 pathways for each gene was selected according to padjust.

Clinical samples
A total of 255 primary NPC samples and adjacent normal nasopharyngeal samples were collected from the First Affiliated Hospital of Kunming Medical University.None of the patients had received chemotherapy or radiotherapy before surgery.All patients involved in this study were known the purposes of this study and wrote the informed consent.Our study was approved by the Ethics Committee of the First Affiliated Hospital of Kunming Medical University.The clinical characteristics of these patients were summarized in Table 1.All samples were quickly frozen in liquid nitrogen and used for subsequent experiments.

Statistical analysis
In this study, all statistical analyses and visualized were performed using R software version 3.4.4and GraphPad Prism 8 software.The continuous variables were shown as mean ± standard deviation (SD), Chi-square test was used to analysis the significance of difference of the categorical variables.And survival analysis was performed using Kaplan-Meier plots and log-rank tests.p-value <.05 was considered statistical significance.

Construction of two m6A-related subclasses in NPC
The systematic analyses of this study have been illustrated in Figure 1a.We obtained the mRNA expression files and corresponding clinical information of 88 NPC patients from the GEO database.Twenty-one m6A regulators were identified in this study from the above dataset.According to the results of Univariate Cox, we found 10 regulators (IGF2B1, ALKBH5, YTHDF2, ELAVL1, WTAP, LRPPRC, HNRNPA2B1, CBLL1, RBM15B, YTHDF1) significantly associated with PFS (Figure 1b, Table 2).Therefore, we clustered NPC patients into different subclasses according to the expression of the 10 m6A regulators who with prognostic values.As shown in   3. T-SNE was performed to validate the above clustering by reducing the dimension of the feature, resulting in that the distinct assignment similar to consensus clustering (Figure 1d).Besides, we performed the clustering analysis in another independent dataset (GSE53819), the results indicated that there are also two distinct subclasses of NPC in the GSE53819 dataset (Figure S1B).And the t-SNE results are concordant with consensus clustering results (Figure 1e).Moreover, a SubMap analysis was conducted to validate the correlation between the subclasses from different datasets, the results indicated that C1 and C2 subclasses in the GSE102349 dataset were significantly correlated with the subclasses in the GSE53819 dataset (Figure 1f, Figure S1C).These results suggested that there were two distinct subclusters of NPC with the different m6A regulator patterns.We also found 10 m6A regulators upregulated in cluster 2 compared with cluster 1 (Figure 1g, Table S1).Kaplan-Meier PFS indicated the clinical relevance of m6A regulators of two m6A-based subclasses in NPC, suggesting that NPC patients belong to cluster 2 with poor prognosis than the patients in cluster 1 (Figure 1h).

Characteristics analysis of TME cell infiltration in two m6A-related subclasses
GSVA was performed to investigate variation pathways among distinct m6A subclasses.We firstly explored the pathways enriched in each subclass, as shown in Figure 2a and Table S2, there were distinct pathways enriched in two m6A-related subclasses.The special genes of cluster 1 were primarily enriched in 46 pathways, most of those pathways were immune-related pathways and some metabolism pathways.The significant immune-related pathways which enriched in immune-related pathways included the B cell receptor signaling pathway, T cell receptor signaling pathway, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, and Toll-like receptor signaling pathway, antigen processing, and presentation.And the metabolism pathways included fatty acid metabolism, betaalanine metabolism, tyrosine metabolism cytochrome P450, porphyrin, and chlorophyll metabolism.However, the special genes of cluster 2 enriched in 14 pathways, including the stromal and tumorigenic pathways, such as ECM receptor interaction, Wnt signaling pathway, adheres junction, Hedgehog signaling pathway.These results indicated that stromal activation in cluster 2 might inhibit the antitumor effects of immune cells and increase immune escape.Whereas immune-related pathways activation in cluster 1 might imply the sensitivity to antitumor therapy.Moreover, we also identify 301 pathways enriched in cluster 1 and 211 pathways enriched in cluster 2 based on Reactome database (Table S3), that accordance to above results  which were the special genes of cluster 1 primarily enriched in immune-related pathways and the special genes of cluster 2 mostly enriched in stromal-related pathways.Additionally, GO enrichment analysis was performed using GVSA, resulting in the special genes of cluster 1 enriched in immune-related biological processes, such as antigen processing and presentation endogenous lipid antigen via MHC class IB, B cell differentiation, T cell activation, immune response, positive regulation of macrophage chemotaxis, positive regulation of T help 2 cell differentiation, etc. Inversely, the special genes of cluster 2 enriched in negative regulation of B cell differentiation, stem cell fate commitment, etc. (Figure 2b, Table S4).
TME consists of tumor cells, diverse immune and stromal cells, which have been reported to the significance in predicting patient outcomes and treatment efficacy. 54Therefore, we investigated the important biological pathways in TME, including epithelialmesenchymal transition (EMT), angiogenesis, Wnt signaling pathway, and CD8 T effector, DNA damage repair, and nucleotide excision related pathways from previous research. 48And the results showed that the special genes of cluster 1 enriched in CD8 T effector, TLR signaling, TCR signaling, immune checkpoint, INF-γ, antigen processing machinery, CLR signaling, and EMT1 pathways.However, the special gene of cluster 2 enriched in Homologous recombination, Wnt signaling, DNA damage repair pathways (Figure 2e, Table S6).These results suggested that cluster 1 might elucidate high immunoreactivity, but cluster 2 might imply high immunosuppression in TME.

Immune infiltration characterization in two m6A-related subclasses
Based on previous analyses, we further explored the immune landscape of m6A-related subclasses in NPC.We obtained the ImmuneScore, StromalScore, and ESTIMATEScore of each sample, resulting in ImmuneScore and StromalScore being higher in cluster 1 than cluster 2 (Figure 3a-c, Table S7).We further used ssGSEA to identify the distinction of the infiltrating immune cells in two m6A-related subclasses.Twenty-four immune cell types were divided into two subclusters and shown in Figure 3d, e and Table S8, we observed that B cells, CD8 T cells, cytotoxic cells, dendritic cells (DCs), eosinophils, immature DCs (iDCs), mast cells, NKCD56dim cells, plasmacytoid DCs (pDCs), T cells, T helper cells, central memory T (Tcm) cells, effector memory T (Tem) cells, follicular helper T (TFH) cells, T help (Th) 1 cell, Th17 cells, Treg cells increased cluster 1 compared to cluster 2, whereas natural killer (NK) cells and Tgd cells elevated in cluster 2 than cluster 1.These results showed that there were significant differences in immune cells between two m6A-related subclasses, which suggested that m6A modification is associated with the immune landscape in NPC.
Also, the differences in the immune status of m6A-related subclasses in NPC was assessed by the CIBERSORT algorithm.The stacked abundance of 22 immune cells was showed in Figure 3f.As result, B cells memory and T cells CD8 were more enriched in cluster 1 than cluster 2, while Macrophage M0 and T cells CD4 memory resting were less enriched in cluster 1 than cluster 2 (Figure 3g, Table S8).The infiltration of the co-exist immune cells from two immune infiltration analysis, i.e.B cells memory (CD20 and CD27), 55 T cells CD8 (CD4 and CD 8), Macrophage M0 (CD68) and T cells CD4 memory resting (CD4 and CD45RO) 56 was detected by IHC in tumor and corresponding para-tumor tissues.The positive rate of CD4, CD45RO, and CD68 was higher, while the positive rate of CD20, CD27 and CD8 was lower in tumor tissues than those in corresponding non tumor tissues, which was coincident with the result of immune infiltration analysis (Figure S3).

Correlation of the two m6A-related subclasses and the immunotherapy responses
According to the previous results, we speculated that cluster 1 was clustered as immune activation phenotype, and cluster 2 was clustered as immune suppression phenotype.Thus, we further test the expression of HLA family and immune checkpoints between two m6A-related subclasses, and the results indicated that HLA-A/B/C/E/F/DRB1/DPB1/DRB5/DPA1/ DMA/DMB high expression in cluster 1 than cluster 2 (Figure 4a, Figure S4, Table S9), as well as the expression of CD80, CD86, CTLA4, HAVCR2, IDO1, LAG3, PDCD1, TIGIT, TNFTSF9 upregulated in cluster 1 compared to cluster 2 (Figure 4b, Table S10).We further detected the TIDE score and found a low TIDE score in cluster 1, but a high TIDE score in cluster 2 (Figure 4c, Table S11).We also compare the response to immunotherapy between two clusters by Submap, the results indicated that the patients in cluster 1 showed a better response to anit-PD1 and anti-CTLA4 treatment (Figure 4d).These findings suggested that the patients in cluster 1 may respond better to immunotherapy treatment.

Correlation of the two m6A-related subclasses and chemotherapy response
Since m6A molecular patterns of NPC are associated with the response of treatment and clinical outcome, the relationship between m6A-related subclasses and chemotherapeutic efficacy was explored.We performed Ridge's regression to predicate the sensitivity of each sample of two subclasses to 139 common chemotherapeutic drugs.The IC50 values were calculated and demonstrated that the patients in cluster 1 showed a better response to 17 chemotherapeutic drugs, including BMS.708163, ATRA, Temsirolimus, Methotrexate, BMS.536924, AKT inhibitor III, RDEA119, PF.02341066, DMOG, CI.1040, Bryostatin.1, AICAR, PD.0332991, Lenalidomide, NVP.TAE684, AZD6244, and LFM.A13 (Figure S3), especially BMS.708163, ATRA, Temsirolimus, Methotrexate (Figure 5).These findings indicated that m6A modification affected the chemotherapeutic outcome in NPC.

Identification of the DEGs and their function enrichment in two m6A-related subclasses
To further explore the potential regulatory mechanism of m6A modification in NPC, we identified 1558 DEGs (928 upregulated and 630 downregulated) between two m6A-related subclasses (Figure 6a, b, Table S12).And GO analysis indicated that the biological processes of those DEGs mainly involved in biological behaviors, such as immune cell T cell activation, regulation of immune effector process, mononuclear cell, lymphocyte, leukocyte, lymphocyte, mononuclear cell proliferation, leukocyte cell, or regulation of leukocyte cell-cell adhesion, positive regulation of immune effector process (Figure 6c, Table S13).And the CC terms enrichment included the external side of the plasma membrane, presynapse, motile cilium, cytoplasmic region, plasm membrane-bounded cell projection cytoplasm, the anchored component of membrane, axoneme, ciliary plasm, motile cilium, axonemal dynein complex (Figure 6c, Table S13).Besides, the molecular function of those DEGs involved in receptor-ligand, signaling receptor activator, metal ion transmembrane transporter, gate channel, purinergic nucleotide receptor, nucleotide receptor, and G-protein-coupled purinergic nucleotide receptor activity, carbohydrate and glycosaminoglycan binding, and ATPdependent microtubule motor activity, minus-end-directed (Figure 6c, Table S13).Moreover, KEGG pathway enrichment analysis indicated that those DEGs significantly enriched in 16 pathways, most of them were immune-related pathways, such as hematopoietic cell lineage, cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway, complement and coagulation cascades, the intestinal immune network for IgA production, natural killer cell-mediated cytotoxicity, and other pathways, such as cell adhesion molecules, calcium signaling pathway, Metabolism of xenobiotics by cytochrome P450 (Figure 6d, Table S14).These results suggested that m6A modification affected tumor progression of NPC mainly through modulating immune-related pathways.

Identification of m6A-related prognostic signature in NPC
We constructed the prognostic signature using univariate Cox and LASSO regression analyses to determine the prognostic value of m6A-related DEGs in the training set from the GSE102349 dataset.As shown in Table 4, 160 DEGs have been found associated with PFS.LASSO regression analysis was then performed to reduce the number of DEGs (Figure 7a, b).And REEP2, TMSB15A, DSEL, and ID4 were screed out and used to construct the risk model (Figure 7c).Then, the risk score of each sample in the training set, test set, and the whole set was calculated (Table S15-17), and all samples of each set were divided into a high-risk score and low-risk score groups, and  the high-risk score associated to high mortality rate both in training, test, and whole sets (Figure 7d-f).Then, Kaplan-Meier curves were plotted with log-rank test and ROC curves to detect the sensitivity and specificity of the risk score to predicate PFS.Kaplan-Meier curves predicated the poor survival rate of patients with the high-risk score than patients with low-risk score (Figure 7g-i).Moreover, ROC curves showed high sensitivity and specificity of a risk score to predicate PFS of NPC patients, the AUC was 0.988, 0.879, 0.784 at 1, 2, 3-year in the training set, the AUC was 0.829, 0.802, 0.802 at 1, 2, 3-year in the test set, the AUC was 0.915, 0.861, 0.772 at 1, 2, 3-year in whole set (Figure 7j-l).These results indicated and encouraged the sensitivity and reliability of the risk model for prognostic predicting.We also collected the NPC samples to validate the expression of m6A-related prognostic signature (DSEL, ID4, REEP2, and TMSB15A) in NPC.As shown in Figure 8a, b, DSEL, ID4, REEP2, and TMSB15A were upregulated in the NPC samples compared with the nontumor samples.However, the expression of DSEL, ID4, REEP2, and TMSB15A not related the tumor progression (Figure 8c).In addition, the Kaplan-Meier curves indicated that high expression of ID4 showed poor survival than low expression of ID4 in NPC (Figure 8d).These data indicated that DSEL, ID4, REEP2, and TMSB15A could be used as the biomarkers for NPC.

Development of a predicated model for NPC
To further validate the predictive accuracy prognostic signature, we developed a nomogram to predicate the survival probability at 1, 2, 3 years.As shown in Figure 9a, a high total score was calculated to predicate the poor survival probability at 1, 2, 3 years.And calibration curves were plotted to evaluate the accuracy and specificity of the prognostic  nomogram, and the plots showed that the prediction lines of 1, 2, 3-year survival closed to the observation lines (Figure 9b-e), suggesting that the excellent accuracy of the prognostic nomogram.In addition, Pearson's correlation analysis was performed to test the relationship between prognostic signature and m6A regulators, resulting in ID4 positively related to 10 m6A regulators, REEP2 positively related to 9 m6A regulators excepted LRPPRC, TMSB15A positively related to ALKBH5, CBLL1, ELAVL1, RBM15B, YTHDF1, YTHDF2, and DSEL positively related to IGF2BP1 (Figure 9f).These results indicated that m6A modification of prognostic signature might affect the clinical outcome of NPC patients.

GSEA and TFs prediction
To further investigate the signaling pathways and potential biological mechanisms for hub genes in NPC, GSEA was used.When threshold was set as padjust < 0.05 and |NES| >     10d, Table 8, Table S25-26).Also, Transcription factors were predicted, a total of 21 TFs and 26 cooperative relationships between m6A-related prognostic signature and TFs was shown in Figure 10e, i.e. 4 for REEP2, 10 for ID4, 3 for DSEL, 8 for TMSB15A (Table S18).

Discussion
In this study, NPC patients were clustered as the two m6A-related molecular patterns based on 21 m6A regulators expression.And the characteristics of each cluster analysis indicated that cluster 1 associated immune activation and cluster 2 was associated with stromal activation.And the special genes of cluster 1 are involved in the immune-related pathways and metabolism pathways, and the special genes of cluster 2 are involved in the stromal-related pathways and tumorigenic pathways.The immune-related pathways involved in cluster 1 included the B cell receptor signaling pathway, T cell receptor signaling pathway, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, and Toll-like receptor signaling pathway, antigen processing, and presentation.B and T cells are the key players of the adaptive immune system and are recognized and activated by the B cell receptor (BCR) and T cell receptor (TCR). 57BCR acts as a key role at multiple checkpoints of B cell biology, 58 and the BCR signaling pathway exerts a crucial role in normal B cell development and adaptive immunity, targeting BCR signaling may have anticancer activity beyond B cell malignancies. 595][66][67] Those immune activation-related pathways enriched in cluster 1, elucidate high immunoreactivity in cluster 1.Nevertheless, the special genes of cluster 2 enriched in stromal-related pathways, such as ECM receptor interaction, Wnt signaling pathway, adheres junction, Hedgehog signaling pathway.ECM is a noncellular component of TME stromal that forms as a scaffold in the tumor and responds to promoting tumor aggression. 68The Wnt signaling pathway is a high conserved signaling pathway that involves embryonic, organ development, and cancer progression. 691][72] Hedgehog (Hh) signaling pathway also acts as an essential role in maintaining CSC self-renewal and tumor progression, targeted treatment of Hh signaling pathway emerges the novel therapeutic strategy in a wide range of solid cancers. 73,74herefore, we speculated that cluster 2 of NPC patients might associate with antitumor immunosuppress.Numerous pieces of evidence have highlighted that m6A modification exerts an indispensable role in regulating the immune system by maintaining the naïve and pluripotent status of immune cells. 75Crosstalk between m6A modification and the immune system has the great significance to discover novel pathogenic mechanisms and to conduct promising therapeutic targets for tumors. 76M6A modification involvement in immune modulation in tumors includes regulating innate and adaptive immune cells 77 and anti-tumor immunity. 78herefore, we further investigated the immune landscape of m6A-related classes.And the results are consistent with the previous finding, B cells memory and T cells CD8 were more enriched in cluster 1 than cluster 2, while Macrophage M0 and T cells CD4 memory resting were less enriched in cluster 1 than cluster 2. ESTIMATE algorithm estimated ImmuneScore, StromalScore, and ESTIMATEScore of each sample, and results indicated that high ImmuneScore and StromalScore of cluster 1 than cluster 2. As known to us, the immune system can be recognized as an innate and adaptive immune system, that contributes to recognize and remove foreign pathogens and tumors. 79Adaptive immunity mainly includes T and B lymphocytes, which harbor plenty of TCRs and BCRs, that respond to different antigens in the pathological process. 80As expected, the immunotherapy response exploring indicated that the patients in the cluster exhibited a better response to anti-PD1 and anti-CTLA4 treatment.Additionally, we also found that the patients in cluster 1 were more sensitive to BMS.708163, ATRA, temsirolimus, and methotrexate.BMS.708163 is a gamma-secretase inhibitor and exerts antitumor effects in lung cancer by reversing EGFR inhibitor resistance, 81 but it currently has not been reported for NPC treatment.ATRA is an active metabolite of vitamin A and used as an effective antitumor in acute promyelocytic leukemia, 82  And emerges the antitumor effects in other tumors by a combination of ATRA with other agents, including chemotherapy, epigenetic modifiers, and others. 83ATRA also has been demonstrated that exert the antitumor effects on NPC. 84Temsirolimus is an mTOR inhibitor widely used in renal cell carcinoma and achieves good efficacy in the therapeutic course, 85 it can rarely use in NPC treatment.Methotrexate (MTX) is widely used for antitumor therapy in a variety of childhood and adult cancers, 86 it has been successfully treated NPC using an effective delivery method. 87he DEGs between two m6A-related subclasses were identified and the mechanism analysis indicated that DEGs mainly enriched immune-related pathways, suggesting that m6A modification regulated tumor progression by modulating immune-related pathways.The m6A-related prognostic signature trained in the GEO dataset was included REEP2, TMSB15A, DSEL, and ID4.The prognostic risk model and predicated nomogram for NPC were construed to estimate the prognosis of NPC.

Conclusion
Taken together, our study demonstrated the significance of m6A modification in NPC by regulating TME.M6A modification patterns caused the heterogeneity and complexity of individual TME, resulting in different responses to therapy.Our findings contributed to our understanding of the pathogenic mechanisms of TME cell infiltrating and provided promising therapeutic targets for NPC.Although our research revealed the significance for understating and conducting treatment strategy of NPC, limitation of this study is small size of the NPC samples to validate the biomarkers, in the subsequent analyses, we will collect the larger NPC samples with corresponding clinical information to demonstrate our results.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Figure
Figure1cand FigureS1A, all patients in the GSE102349 dataset were divided into two distinct subclasses based on the molecular pattern of m6A regulators by consensus clustering, including cluster 1 (64 patients) and cluster 2 (24 patients).The clinical characteristics of 88 patients from the GSE102349 dataset were shown in Table3.T-SNE was performed to validate the above clustering by reducing the dimension of the feature, resulting in that the distinct assignment similar to consensus clustering (Figure1d).Besides, we performed the clustering analysis in another independent dataset (GSE53819), the results indicated

Figure 1 .
Figure 1.Construction of two m6A-related subclasses in NPC A. The workflow of this study.B. Forest plot showing the results of the Univariate Cox regression analysis between m6A regulators and PFS.C. Consensus matrices of the GSE102349 dataset for k = 2. D-E.T-SNE analyses show the stratification into two m6A-related subclasses of GSE102349 and GSE53819 datasets.F. Submap matrix analysis validated the similarity of molecular classes from GSE102349 and GSE102349 datasets.G.The expression of 10 prognostic-related m6A regulators (IGF2B1, ALKBH5, YTHDF2, ELAVL1, WTAP, LRPPRC, HNRNPA2B1, CBLL1, RBM15B, YTHDF1) between two m6Arelated subclasses.H. Kaplan-meier survival plot showing the PFS between two m6A-related subclasses.

Figure 2 .
Figure 2. Characteristics analysis of tme cell infiltration in two m6a-related subclasses a. heatmap showing gsva enrichment analysis of the various pathways between two m6a-related subclasses.b. heatmap showing gsva enrichment analysis of the biological processes between two m6a-related subclasses.C. heatmap showing gsva enrichment analysis of 50 hallmark cancer-related pathways between two m6a-related subclasses.D. heatmap showing gsva enrichment analysis of oncogenic and TME cell-infiltrating related pathways between two m6a-related subclasses.E.boxplot showing the distinction of oncogenic and TME cell-infiltrating related pathways between two m6a-related subclasses.

Figure 3 .
Figure 3.Immune infiltration characterization in two m6A-related subclasses a-c.Violin plot indicating stromal score, immune score and estimate score between two m6A-related subclasses.d.Heatmap showing the TME cell infiltrating between two m6A-related subclasses.e. Boxplot showing distinct TME cell infiltrating between two m6A-related subclasses.f.The staked abundance of immune cells between two m6A-related subclasses.g.Differences in immune cells between two m6A-related subclasses.

Figure 5 .
Figure 5. Correlation of the two m6A-related subclasses and chemotherapy response a-d.Violin plot showing IC50 value of BMS.708163, ATRA, temsirolimus, and methotrexate for each patient from two m6A-related subclasses.

Figure 6 .
Figure 6.Identification of the DEGs and their function enrichment in two m6A-related subclasses a. Volcano plot showing the DEGs between two m6A-related subclasses.b.Heatmap showing the DEGs between two m6A-related subclasses.c.Bubble plot indicating the GO enrichment of DEGs between two m6A-related subclasses, including top 10 biological processes (BP), top 10 cellular components (CC), and top 10 molecular functions (MF) terms.d.Bubble plot indicating the KEGG pathways enrichment of DEGs between two m6A-related subclasses.

Figure 7 .
Figure 7. Identification of m6A-related prognostic signature in NPC A. Least absolute shrinkage and optimal LASSO coefficients of the 160 DEGs.B. Determination of the number of factors by the LASSO regression analysis.C. The sheet of showing the LASSO coefficients of four genes (REEP2, TMSB15A, DSEL, and ID4).D-F.The differences between high-risk score and low-risk score groups both in the training set, test set, and whole set.Upper: the distribution of risk core.Middle: the survival status of patients both in the high-risk score and low-risk score groups.Bottom: heatmap showing expression of m6A-related prognostic signature both in the high-risk score and low-risk score groups.G-I.Kaplan-Meier survival plot showing the PFS between high-risk score and low-risk score groups in the training set, test set, and the whole set.J-L.ROC curves indicate the 1-, 2-, 3-years survival of NPC patients both in the training set, test set, and whole set.

Figure 9 .
Figure 9. Development of a predicated model for NPC A. A predicating nomogram predicting PFS of NPC patients.B-E.Calibration plots assess the accuracy and specificity of the prognostic nomogram.F.Heatmap showing the correlation between prognostic signature and 10 m6A regulators.

Table 1 .
Clinical characteristics of the NPC patients.

Table 2 .
Identification of 21 m6A regulators with prognostic significance in the GSE102349 dataset based on Univariate Cox analysis.

Table 3 .
Clinical characteristics of NPC patients from GSE102349 dataset in distinct subclasses.

Table 4 .
Univariate Cox regression analysis of m6A-related DEGs associated to survival in GSE102349 dataset.