Gene expression data and clinical information for melanoma patient samples were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and TCGA database (https://portal.gdc.cancer.gov/), which are publicly and freely available. We screened the melanoma dataset and eventually included a total of 1,020 patients in the study for subsequent analysis, including the GSE19234, GSE22154, GSE50509, GSE59455, GSE65904, GSE22153, GSE54437, and TCGA-SCKM datasets. For data pre-processing, we downloaded the “CEL” files from the GEO database, relying on the “affy” and “simpleaffy” R packages for background correction and normalization, while the RNA sequencing data from the TCGA database is downloaded in normalized FPKM format and then converted to transcripts per kilobase million (TPM) format. With reference to previous experience, the “ComBat” method of the “sva” R package was used to reduce the batch effect between different datasets, which was mainly a problem for datasets from the GEO database (Dai et al., 2018). Both somatic mutation data and copy number variation data of TCGA-SCKM were curated from the UCSC Xena database and Davoli et al. (Davoli et al., 2017). The copy number variation of 23 m⁶A regulators in human chromosomes was mapped by the “RCircos” R package. For non-synonymous mutations such as frameshift mutations, inflammatory mutations, missense mutations, nonsense mutations, and splice site mutations, numbers represent the tumor mutational load (TML). Supplementary Table S1 presents the clinical information of the samples from the meta-GEO and TCGA-SKCM databases.

A literature review of m⁶A methylation modifications revealed that there are now 23 recognized m⁶A regulators, which constitute the modification pattern of m⁶A methylation (Zhao et al., 2017; Chen et al., 2019a; He et al., 2019; Zaccara et al., 2019). Specifically, eight writers include CBLL1, KIAA1429, METTL14, METTL3, RBM15, RBM15B, WTAP, and ZC3H13; two erasers include ALKBH5 and FTO; 13 readers include ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, LRPPRC, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3; and 23 of them shared the key task of m⁶A methylation modification. Using non-negative matrix decomposition (NMF), we performed a clustering analysis of the 23 m⁶A regulators based on their expression, which could identify different types of m⁶A modification patterns. The expression matrix A of the 23 m⁶A regulators was first split into non-negative matrices W and H, as A ≈ WH, and then the matrix A was subjected to repeated factorization, and finally the output data was summarized, which gave the clustering results of the melanoma samples. It is crucial to consider factors such as covariance, dispersion, and silhouette coefficient to determine the optimal number of clustering groups. The “NMF” R package to perform the clustering analysis used the “brunet” and “200 nruns” algorithms.

The Hallmarker gene set (Subramanian et al., 2005) and Mariathasan et al. (Mariathasan et al., 2018) constructed gene set were used as well-defined biometric backgrounds for gene set variation analysis (GSVA) with “GSVA” R package (Hänzelmann et al., 2013), which was designed to explore the variation in biological processes across different m⁶A modification patterns. In the gene ontology (GO) analysis, we annotated the functions of 23 m⁶A regulators under three entries of cellular component (CC), molecular function (MF), and biological process (BP), which was done using the “clusterProfiler” R package. For GSVA and GO analysis, the cut-off value was set to a false discovery rate (FDR) < 0.01.

In quantifying the relative abundance of 28 immune cell types curated by Charoentong et al. in the tumor microenvironment, we refer to recent studies using the single sample gene enrichment analysis (ssGSEA) method, which well marks the specific functional gene panels of each immune cell type (Charoentong et al., 2017; Jia et al., 2018). As in the previous study (Chen et al., 2020a), we expressed the relative abundance of various immune cell types in the form of enrichment scores, and they were normalized to a uniform distribution from 0 to 1. In terms of biosimilarity, infiltrating immune cells were evaluated and acted upon using multidimensional scaling (MDS) and Gaussian fitting models, and moreover, the deconvolution approach CIBERSORT (Newman et al., 2019) (http://cibersort.stanford.edu/) was then used to estimate the abundance of 22 different subpopulations of leukocytes, which have melanoma gene expression profiles.

T cell-inflamed gene expression profile (GEP) is a superior predictor of response to anti-PD-1 regimens, which contained IFN-γ-responsive genes related to antigen presentation, cytotoxic activity, and adaptive immune resistance (Ayers et al., 2017). The T cell-inflamed scores were calculated and weighted by averaging of the included genes for the IFN-γ (6-gene) and expanded immune (18-gene) signatures. In modeling different types of tumor immune evasion mechanisms, we drew on the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm proposed by Jiang et al. (Jiang et al., 2018). This algorithm integrates the dysfunction of tumor-infiltrating toxic T lymphocytes (CTLs) and rejection of CTLs by immunosuppressive factors. The higher TIDE score implies greater chance of immune escape of tumor cells and represents a possible poor outcome of treatment with ICIs. The method of Estimation of Stromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) (Yoshihara et al., 2013) was adopted to calculate the immune score of tumors, and this algorithm can be better based on transcriptional profiles to estimate the cellularity of the tumor and the purity of the tumor. The level of infiltrating immune and stromal cells is the basis for tumor purity, which is predicted by the immune score of the tumor. In detail, a high immune score of a tumor is an indication of a high infiltration of immune cells in the tumor tissue, or a low tumor purity.

The MutSigCV algorithm was used to identify significantly mutated genes (SMGs) (Lawrence et al., 2013; Chen et al., 2019b), which takes into account the specific background mutation rate in the mutation context before evaluating the significant enrichment of non-resting somatic mutations in a gene. We considered q < 0.1 as statistically significant, and these genes needed to be certified in the Cancer Cell Line Encyclopedia of Humans (CCLE) (Ghandi et al., 2019) to be defined as SMGs (Chen et al., 2020c) (Supplementary Table S2). The “maftools” R package (Mayakonda et al., 2018) was used to characterize genes in the TCGA-SKCM cohort that underwent m⁶A modification, the mutation details of SMGs, and the capture of mutational features in the genomic data. The ExtractSignatures function based on Bayesian variation non-negative matrix decomposition was used for model construction; specifically, using this function, we split the mutation portrait matrix into two non-negative matrices and noted as “signature” and “contribution,” where “signature” represents the mutation process and “contribution” represents the corresponding mutation activities (Chong et al., 2021a). Better still, the SignatureEnrichment function allows determining the optimal number of extracted mutation features and assigning them appropriately to each sample. For comparison and annotation, using the Catalogue of Somatic Mutations in Cancer (COSMIC) (Kandoth et al., 2013) as a reference, we performed a cosine similarity analysis on the extracted melanoma mutation portraits.

Patients were classified into three clusters of m⁶A modification patterns using a consensus clustering algorithm, and then the “limma” R package (Ritchie et al., 2015) was used to find differentially expressed genes between groups. Voom normalized data were then subjected to “lmFit” and “eBayes” function algorithms, which in turn allowed the calculation of specific data for differential expression. In this process, we set adjusted p-values <0.001 as statistically significant differences.

Based on principal component analysis (PCA), we constructed an m⁶A score system to quantify the level of m⁶A modifications in specific patients. According to DEGs, they are the intersecting parts of different m⁶A clusters, and we analyzed the prognostic impact of each gene on melanoma patients with the help of univariate Cox regression models. Deeper feature selection was performed for genes that significantly affect the prognosis of melanoma patients, and this process was computed by the recursive feature elimination (RFE) method of random forest and the 10-fold cross-validation method included in the “caret” R package. Further, we obtained the gene expression profiles based on the above steps, and the principal components 1 and 2 obtained from PCA analysis were the basis of our feature score. The specific formula for this score system is referred to a previous study (Zhang et al., 2020; Chong et al., 2021b), m⁶Sig score = ∑(PC1i + PC2i).

The gene expression profiles of patients treated with ICI were retrieved in publicly available databases, focusing on matching with clinical information. Ultimately, we included metastatic melanoma treated with PD-1 (nivolumab or pembrolizumab) or PD-1 combined with CTLA-4 (ipilimumab) (Liu et al., 2019), and metastatic urothelial carcinoma (mUC) treated with atezolizumab (anti-PD-L1 mcAb) (Mariathasan et al., 2018) in this study. The gene expression profiles of the samples were converted in TPM format.

All statistical analyses in the study were performed with R 3.6.1. Student’s t-test was performed for quantitative data conforming to a normal distribution, and Wilcoxon rank sum test was performed for non-normally distributed data. When more than two sets of analyses were performed, the nonparametric test was the Kruskal-Wallis test, while the parametric test was the analysis of variance (Hazra and Gogtay, 2016). The Fisher exact test was used for the calculation of contingency rates. Kaplan-Meier survival analysis and Cox regression analysis were performed using the “Survminer” package, and the m⁶Sig score subgroup stratum was “survival” package with the surv-cutpoint function completed. “timeROC” package completed the evaluation of the m⁶Sig score model, which plotted the corresponding subject operating characteristic curve (ROC) and calculated the area under the curve (AUC). In analyzing the relationship between patient’s clinical characteristics and the m⁶Sig score system, multivariate regression models were used to adjust for confounding factors in this. p < 0.05 was considered as statistical significance, and the Benjamini-Hochberg method was used to perform multiple hypothesis testing for false discovery rate (FDR) (Love et al., 2014).

In our study, we explored the possible physiological roles of 23 m⁶A methylation-regulated genes in melanoma, including the “writers” CBLL1, KIAA1429, METTL14, METTL3, RBM15, RBM15B, WTAP, and ZC3H13; the “readers” ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, LRPPRC, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3; and the “erasers” ALKBH5 and FTO. These m⁶A regulators not only recognize, remove, and add m⁶A modification sites but also, as revealed by GO enrichment analysis and Metascape analysis, can actually alter biological processes, such as regulating mRNA stability, RNA modifications, and RNA metabolism (Figure 1A). Among 467 melanoma patient samples with genomic sequencing, 133 (28.48%) had somatic mutations in m⁶A regulators, which mainly included missense mutations, nonsense mutations, and code-shifting mutations (Figure 1B). KIAA1429 had the highest mutation frequency, followed closely by IGF2BP1, and the next in the gradient were YTHDC1, LRPPRC, YTHDC2, ZC3H13, YTHDF1, and IGF2BP3. Interestingly, IGF2BP1, ZC3H13, and YTHDF1 had only missense mutations in the relatively high mutation frequencies. Analysis of the co-mutation profiles of the 23 m⁶A regulators revealed significant co-mutations between FMR1 and IGF2BP1, IGF2BP2 and IGF2BP1, ZC3H13 and LRPPRC, YTHDF3 and RBM15, ALKBH5 and METTL3, and ALKBH5 and FMR1 (Figure 1C). When performing CNV mutation analysis, we concluded that there was widespread CNV amplification in IGF2BP1, YTHDF1, and KIAA1429, while CNV deletion was more widespread in WTAP and RBM15 (Figure 1D). Comparing primary melanoma and metastases, we found that ALKBH5, ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1/2/3, KIAA1429, LRPPRC, RBM15, YTHDC1/2, YTHDF1/3, and ZC3H13 were significantly upregulated in metastases, while RBM15B and METTL3 were significantly upregulated in primary melanoma (Figure 1E). The expression of m⁶A regulators with aberrant CNV amplification was also upregulated in metastases compared to primary melanoma (ALKBH5, FMR1, HNRNPA2B1, IGF2BP1/2/3, KIAA1429, YTHDC1, YTHDF1/3), and conversely, m⁶A regulators with aberrant CNV deletion were also downregulated (METTL3, RBM15B, YTHDC2), which are obtained by combining Figure 1D. The m⁶A regulator network mapped in Figure 1F showed the interaction relationships between 23 molecules that are interconnected and influence each other, which further modulates the prognosis of melanoma patients. This implies that there is a complex and well-organized crossover network between the regulators of writers, readers, and erasers, and this network allows the m⁶A modification pattern to further refine and take effect, influencing the development and metastasis of melanoma. Using the Spearman correlation test, we found that there is a mutual regulatory relationship between these m⁶A regulators. Interestingly, ALKBH5 was negatively correlated with most of the m⁶A regulators, while FMR1 and HNRNPA2B1 were positively correlated with most of them (Supplementary Figure S1A). We further analyzed the association between tumor purity and 23 m⁶A modification regulators (Figure 1G) and found that most of the m⁶A regulators were positively correlated with tumor cell purity, whereas WTAP has a negative association, suggesting that WTAP was enriched in non-tumor cell components. Forest plots with Cox regression model were employed to speculate the relationship between m⁶A regulators and the prognosis of melanoma patients. We found that samples with high expression of WTAP, FMR1, and METTL14 were associated with improved overall survival, while an opposite tendency was observed in RBM15B and ELAVL1 (Supplementary Figure S1B). Taken together, we integrated the genomic and transcriptomic landscapes of m⁶A regulators in melanoma, and noticed the changes in the expression levels and genetic variation of m⁶A regulators driving the development and progression of melanoma.

The m⁶A Methylation Modification Pattern Consisting of 23 m⁶A Regulators Is Associated with Prognosis in Melanoma Patients.

Further, we stratified melanoma samples into three m⁶A modification patterns according to the expression of m⁶A regulators, a process based on consensus clustering analysis of the NMF algorithm (Supplementary Figure S2A,C). We named the three clusters as m⁶A-C1, m⁶A-C2, and m⁶A-C3, respectively (Figure 2A). The samples of m⁶A-C2 cluster had significantly different from the other two groups with regard to ELAVL1, RBM15B, YTHDF1/2/3, IGF2BP1/2/3, WTAP, METTL3, ZC3H13, RBM15, HNRNPA2B1, CBLL1, and LRPPRC. Besides, YTHDF1, IGF2BP3, METTL3, ZC3H13, and LRPPRC were significantly upregulated in the m⁶A-C1 subtype, WTAP and RBM15 were significantly upregulated in the m⁶A-C3 subtype, while ELAVL1, IGF2BP3, ZC3H13, and LRPPRC were significantly decreased in the m⁶A-C3 subtype. This conclusion was validated in the meta-GEO cohort consisting of five datasets, which include GSE19234, GSE22154, GSE50509, GSE59455, and GSE65904 (Supplementary Figure S2B). Patients in the m⁶A-C3 cluster have a significant survival advantage than other clusters in both TCGA and meta-GEO cohort (log-rank test, TCGA: p < 0.0001, Figure 2B; meta-GEO: p = 0.0015, Figure 2D). This model remained significant after multivariate Cox proportional risk regression analysis adjusted for clinicopathological factors of age, gender, and stage (TCGA: HR = 0.34 (0.21–0.51), p < 0.001; meta-GEO: HR = 0.51 (0.35–0.76), p < 0.001; Figures 2C,E).

We performed GSVA analysis against on Hallmarker gene set in an attempt to discover differences in the biological behavior of the three m⁶A modification clusters. As shown in Figure 3A, m⁶A-C1 cluster was associated with cell proliferation and differentiation and glucose transport, including oxidative phosphorylation, PI3K/AKT/mTOR signaling, DNA repair, and glycolysis. m⁶A-C2 cluster is distinguished by cancer and immune surveillance, involving epithelial mesenchymal transition, TGF-β signal, TNF-α signaling via NF-κB, and IL2/STAT5 signaling. As for m⁶A-C3, it was significantly enriched in signaling pathways related to inflammation and innate immune response, such as interferon-γ response, interferon-α response, allograft rejection, IL6/JAK/STAT3 signaling, and inflammatory response. The GSVA results further corroborate that these three m⁶A methylation modification clusters are directly related to different molecular mechanism, and m⁶A-C3 was strongly associated with antitumor immunity. In addition, we further evaluated the immune enrichment level of m⁶A methylation modification clusters using the ImmuneScore model constructed by ESTIMATE algorithm. The results showed significant differences in different clusters in both TCGA-SKCM cohort and the meta-GEO cohort (Supplementary Figure S3A,B). There is a coherence between the immune activation and survival time, which cluster with higher ImmuneScore having greater survival benefit for patients, like m⁶A-C3 in the TCGA-SKCM cohort and meta-GEO cohort. Thorsson et al. (Thorsson et al., 2018) divided the tumor immune landscape into six immune subtypes, represented with immune infiltration and stromal activation. Consistent with our findings, the m⁶A-C1 cluster is more inclined to the “Proliferation” and “Wound Healing” subtypes, m⁶A-C2 is highly expressed in “TGF-β Response,” and m⁶A-C3 is mainly dominated by “Lymphocyte Infiltration Signature Score,” “Macrophage Regulation,” and “IFN-γ Response” subtypes (Figure 3B). In addition, we performed a comparative analysis of immune checkpoint-related key genes (IDO1, CD274, TIM-3, PDCD1, CTLA-4, LAG3, and PDCD1LG2) among the three clusters. The results indicated the expression levels of seven key genes differ significantly between the three clusters, and the highest expression was all observed in the m⁶A-C3 cluster. TCGA-SKCM has established molecular typing based on the genomic landscape and transcriptomic profile. Although the proportion of m⁶A modification clusters among different mutational-based molecular subtype (BRAF-Mut, RAS-Mut, NF1-Mut, and Triple Negative) were not significant (Figure 3D), an obvious difference was found in transcriptomic-based subtype across three m⁶A modification clusters (Figure 3E). Samples with TCGA-Immune subtype account for 93.1% of the m⁶A-C3 cluster, followed by 64.8% in m⁶A-C2 subtype. However, the m⁶A-C1 were dominated by TCGA-Keratin (55.7%) and TCGA-MIFT-low (32.8%) subtype, which demonstrated the desert-related immune phenotype in m⁶A-C1. It is feasible to determine immune cell type abundance and expression from bulk tissues with digital cytometry (Newman et al., 2019). We also compared the immune cell infiltration level among the three m⁶A modification clusters in Figure 3F. The m⁶A-C1 cluster showed more myeloid-derived suppressor cell (MDSC), regulatory T cells, and T helper cell infiltration, while the m⁶A-C3 cluster exhibited infiltration of most types of T cells, natural killer cells, and dendritic cells. This suggested that the lower survival risk of melanoma patients with m⁶A-C3 clusters may be due to effective activation of the pre-existing immunity to inhibit tumor growth and malignant progression. Likewise, evaluation on immune cell abundance by CIBERSORT algorithm also corroborate our conclusions (Supplementary Figure S3C). In addition, the association between each m⁶A regulator and immune cell infiltration was also explored. As expressed in Supplementary Figure S3D, upregulation of WTAP and ALKBH5 was positively correlated with enhanced immune infiltration, while high expression of LRPPRC, METTL3, YTHDF1/3, and ZC3H13 was mostly associated with immunosuppression.

Since RNA N6-methyladenosine (m⁶A) modification plays an important role in post-transcriptional regulation, we further examine the potential impact on gene expression change of each m⁶A modification cluster in melanoma. To clarify these queries, we employed the Bayesian-based method to identify differentially expressed genes (DEGs) that are differentially regulated across the three m⁶A methylation modification clusters. As illustrated in the Venn diagram of Figure 4A, there are a total of 636 DEGs which may play the crucial role in distinguishing the three m⁶A modification clusters (Supplementary Table S3). Based on these 636 representative DEGs of m⁶A signature, we further stratified the melanoma samples into three well typed transcriptomic phenotypes (denoted as m⁶Sig-SI, m⁶Sig-SII, and m⁶Sig-SIII) by unsupervised consensus clustering analysis (Supplementary Figure S4A). We also compared the m6A clusters and m6A signature-derived subtype, and found a significant association among these two-stratification method (Supplementary Table S4, adjusted χ² test, p < 0.0001). Patients in m⁶Sig-SII were proved to be associated to better prognosis, while m⁶Sig-SI had a worse outcome (p < 0.001, log-rank test; Figure 4B). PD-L1 and ImmmuneScore were also highly expressed in m⁶Sig-SII subtype than the other subtypes (p < 0.0001, Kruskal Wallis test; Figures 4C,D). Mariathasan et al. (Mariathasan et al., 2018) summarized and formed a suit of gene set for assessing the activation of immune and stroma signaling pathway, whereby we adopted to evaluate the m⁶Sig signature. The m⁶Sig-SI subgroup was enriched in cell proliferation and DNA damage repair related pathways, m⁶Sig-SIII was characterized by CD8 T effector and antigen processing related pathway, whereas m⁶Sig-SII was focused in epithelial mesenchymal transition (EMT) related signaling pathway (Figure 4E). The expression level of 23 m⁶A regulators in three gene subgroups was also compared and shown in Supplementary Figure S4B. We observed significant differences of m⁶A regulator expression in the three m⁶A gene-signature subgroups, which was consistent with the expected results of the m⁶A methylation modification clusters.

The results of the previous parts of the study can be concluded that m⁶A methylation modification has a strong association with prognosis and immune regulation in melanoma patients. In order to be able to accurately predict the m⁶A methylation modification cluster of tumors in individual melanoma patients, we developed the m⁶Sig score system to quantify the m⁶A modification cluster based on the identified m⁶A-related signature genes. In Figure 5A, the Sankey diagram uncovered the workflow of the m⁶Sig score system in melanoma. It can be found that m⁶A-C3 was linked to a higher m⁶Sig score, and lower for keratin and MITF-low subtypes. We sought to evaluate the significance of the m⁶Sig score system in predicting the prognosis of patients with melanoma through survival analysis. As shown in Figures 5B,C, the m⁶Sig score system clearly distinguished patient with different prognosis in both the TCGA-SKCM cohort and meta-GEO cohort (patients with high m⁶Sig scores had a better prognosis). We performed the Kruskal Wallis test and showed that the m⁶Sig score could be clearly distinguished in the immune and keratin subtypes (Figure 5D). Encouragingly, the m⁶Sig score system can be extremely well distinguished among the previous m⁶A clusters and m⁶Sig clusters (Supplementary Figure S5A,B), which are clusters constructed based on m⁶A methylation modification clusters and DEGs, respectively. Compared with clinicopathological staging, the m⁶Sig score system can better evaluate the prognosis of melanoma patients (AUCs for stage and m⁶Sig score are 0.613 and 0.681, respectively; Supplementary Figure S5C). In detail, multivariable cox regression model revealed that melanoma patients with low m⁶Sig scores had a worse survival outcome in TCGA (Supplementary Figure S5D) and meta-GEO cohort (Supplementary Figure S5E). This score system was also explored by ImmuneScore, and patients with high m⁶Sig scores had a higher ImmuneScore (p < 0.0001, Supplementary Figure S5F,G). Among the six immune subtypes, the m⁶Sig score was also differentially distributed. Samples with high scores were mostly clustered in the “Macrophage Regulation,” “Lymphocyte Infiltration Signature Score,” and “IFN-γ Response” subtypes (Supplementary Figure S5H), which was similar to the m⁶A-C3 cluster. Heatmap of correlation matrix demonstrated that the m⁶Sig score was markedly positively correlated with the immune activation process and negatively correlated with cell cycle and DNA damage repair (Supplementary Figure S5I). To better evaluate the effectiveness of our m⁶Sig score system in predicting the prognosis of melanoma patients, we introduced two independent cohorts (GSE22153 and GSE54437) to perform a survival analysis, and the results showed that patients with high m⁶Sig score had a better prognosis (Supplementary Figure S6A,B). Furthermore, a survival analysis after combining all patients involved in this study revealed that a high m⁶Sig score continued to indicate a survival benefit in melanoma patients (p < 0.0001, Supplementary Figure S6C). The TCGA-SKCM cohort also demonstrated that melanoma patients with high m⁶Sig score had prolonged disease-free survival (DFS, p = 0.0064, Supplementary Figure S6D). In addition, m⁶Sig score also negatively correlated with somatic copy number alternation (SCNA) level (r = −0.49, p < 0.0001, Figure 5E), which is a significant predictor of immunotherapy resistance in melanoma. The PD-L1 expression levels were also positively correlated with the m⁶Sig score (r = 0.75, p < 0.0001, Figure 5F), suggesting that melanoma patients may also benefit from the m⁶Sig score system for PD1/PD-L1 treatment regimens. More deeply, we performed significant mutation gene (SMG) analysis of melanoma samples based on m⁶Sig score, and the waterfall plot of mutation landscape noted that BRAF (54%/48%), SIRPB1(11%/5%), and KNSTRN (7%/2%) had higher somatic mutation rates in the high-score group, although BRAF was not statistically significant (Figure 5G). These data assist us to more comprehensively understand the m⁶Sig score system mapping to genomic variants, predicting that m⁶A methylation modification is closely linked to somatic mutations in melanoma patients.

The m⁶Sig Score System Can Be a Better Predictor of the Effectiveness of Immunotherapy in Cancer

Cancer treatment regimens based on immune checkpoint inhibitors have provided a landmark innovation in the treatment of malignancies, mostly in melanoma. In addition to TML and PD-L1, TIDE and T cell-inflamed GEP have been recommended to predict immune response in recent years (Chen et al., 2019c; Chen et al., 2019d). We compared the established m⁶Sig score system with the T-cell inflamed gene expression profile (GEP) score and found that melanoma patients with high m⁶Sig score had elevated T-cell inflamed GEP score in both the TCGA-SKCM cohort and the Meta-GEO cohort (p < 0.0001, Figures 6A,B). In contrast, TIDE showed increased levels in patients with low m⁶Sig scores, implying that greater chance of tumor immune escape and resistance in low m⁶Sig scores subgroup (p < 0.0001, Figures 6C,D). These results further demonstrate that m⁶A modification clusters play a critical role in the immune response of tumors, thereby affecting the immune microenvironment of tumors.

Aforementioned data point to a strong association between m⁶A modification and immune response, we next investigated whether the m⁶Sig score could predict patients’ response to ICI treatment in independent immunotherapy cohorts. Patients with high m⁶Sig score exhibited significantly longer survival time (log-rank test, p = 0.0082, Figure 6E) and markedly clinical response to PD-1/CTLA-4 treatment in melanoma (response rate, high vs low m⁶Sig score subgroup, 52.3 vs 31.2%, Figure 6F). This result was also identified in an anti-PD-L1 metastatic uroepithelial cancer cohort (Mariathasan et al., 2018), in which patients with high m⁶Sig scores significantly benefited from PD-L1 immunotherapy (log-rank test, p = 0.0005, Figure 6G; response rate: high vs low m6Sig score subgroup, 32.7 vs 17.8%, Figure 6H). Furthermore, we found that patients with PD-1 immune response also had a higher m⁶Sig score (p = 0.0012, Figure 6I). Metastatic uroepithelial carcinoma patients with immune inflamed phenotype had a higher m⁶Sig score than immune excluded and desert phenotype (Figure 6J). A significant elevation of PD-L1 was identified in high m⁶Sig score subgroup (p < 0.0001, Figure 6K). Therefore, we divided the overall population into four subgroups according to the TMEsig-score and PD-L1 distribution, including TMEsig-score-H + PD-L1-H, TMEsig-score-H + PD-L1-L, TMEsig-score-L + PD-L1-H, and TMEsig-score-L + PD-L1-L. The TMEsig-score-H + PD-L1-H subgroup exhibited the best prognostic outcome compared with the other three subgroups (log-rank test, p < 0.0001, Figure 6L). Based on the results of the above analysis, our established m⁶Sig score system enables the prediction of responsiveness and prognosis to cancer immunotherapy.