The gene expression data and clinicopathological data of 457 melanoma patients and 233 normal skin tissues were obtained from the UCSC Xena platform (http://xena.ucsc.edu/). Gene mutation maps of 457 melanoma patients were downloaded from TCGA (https://portal.gdc.cancer.gov/). Gene expression data and the clinical data of GSE65904 (n = 210), GSE59455 (n = 122), GSE54467 (n = 78), GSE91061 (n = 51), GSE19234 (n = 44), GSE78220 (n = 26), GSE22154 (n = 22) and GSE100797 (n = 21) were obtained from gene expression omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). Other gene expression data and the clinical data were obtained from previous studies, including a melanoma cohort from Nathanson et al. (n = 24), an anti-CTLA4 melanoma cohort (n = 39), an anti-PD1 melanoma cohort (n = 121), and IMvigor210 cohort (n = 298). Finally, 233 normal samples and 1,513 tumor samples were included in the follow-up analyses. The sequencing platform and sample size of each cohort are summarized in Table 1.

A total of 21 regulators related to m6A modification were included in the follow-up analyses. These 21 m6A regulators included eight writers (METTL3, RBM15, METTL14, RBM15B, CBLL1, WTAP, KIAA1429, ZC3H13), eleven readers (YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, IGF2BP1, FMR1, HNRNPA2B1, ELAVL1, HNRNPC, LRPPRC) and two erasers (ALKBH5, FTO) (Zhang and Wu, 2020a).

The limma package was used to compare the difference in gene expression between normal samples and tumor samples. The mutation maps of the TCGA tumor samples were obtained, and the mutation waterfall map of the 19 mutated m6A regulators was drawn with the GenVisR package. The landscape of the CNV variation in the chromosome of these m6A regulators was plotted via R package OmicCircos. We acquired seven cutaneous melanoma driver genes from the previous study, included BRAF, NRAS, CDKN2A, NF1, MECOM, COL5A1, DACH1, in which mutation frequencies were more than 10% (Bailey et al., 2018). Then we performed a permutation test to examine the p-value of the expression of m6A regulators between the driver genes-mutant and -wild samples.

Unsupervised clustering analysis of melanoma patients was performed using the ConsensusClusterPlus package in the R software. We used the GSVA enrichment analysis to explore the differences in pathways between the two m6A clusters. We downloaded the gene sets of “c2.cp.kegg.v7.4” from the MSigDB database for the GSVA enrichment analysis (Subramanian et al., 2005). FDR < 0.05 was considered statistically significant. We combined eight independent melanoma cohorts to validate the unsupervised clustering model, including GSE78220, GSE91061, GSE100797, a melanoma cohort from Nathanson et al. study, GSE65904, GSE19234, GSE22154, GSE59455.

We used the estimate R package to estimate the immune infiltration in TCGA melanoma patients. The 28 immune-related gene sets were obtained from the study of Charoentong et al. They contained both innate immune cells, such as eosinophils, mast cells, macrophages, and so on, and adaptive immune cells, such as CD4⁺ T cell, Gamma delta T cell, and so on (Charoentong et al., 2017a). The single-sample gene set enrichment analysis (ssGSEA) of 28 immune-related gene sets was used to quantify the abundances of 28 immune cells using GSVA package in melanoma patients. The estimate package was used to evaluate the immune scores, stromal scores, and tumor purity of the tumor patients. In addition, we utilized CIBERSORT, MCPcounter, TIMER algorithms to assess the abundances of immune infiltration between the two m6A-clusters. These three analyses were performed via the CIBERSORT function, MCPcounter, and IOBR R packages, respectively. Differential gene expression analysis between the two m6A clusters was performed via the limma package of the R language.

We selected differential expression genes between two m6A clusters (n = 849, p < 0.05 and | log fold change | >0.5). We performed GO enrichment analysis for the differential expression genes via the clusterProfiler package. Then we screened genes present in both the GEO and our total of differentially expressed genes via univariate Cox regression analyses (n = 31, p < 0.0001). We used multivariate Cox regression analyses to narrow the range of target genes. The Akaike information criterion (AIC) was used to measure the prognostic prediction ability of the Cox proportional hazard regression model by calculating for different variables. Here, lower values of AIC indicate a better model fit. Finally, 12 genes were identified to establish a risk signature (Gaulke et al., 2019). The Benjamini and Hochberg method was applied for multiple-testing adjustments, and the p-value remained significant. The patients were divided into the high-risk and low-risk groups using the surv_cutpoint function (Chu et al., 2020). The integrated discrimination improvement (IDI) and the net reclassification improvement (NRI) were used to evaluate the performance of the signature (survIDINRI package). To test the predictive efficiency of this signature on immunotherapy, we performed the validation analysis through integrating five melanoma immunotherapy cohorts, including an anti-PD1 melanoma cohort, an anti-CTLA4 melanoma cohort, GSE91061, GSE78220 and GSE100797 (ICB-therapy-combined melanoma, n = 258). Data were corrected for batch effects via the combat function from the sva R package. To predict the survival rates, we constructed a nomogram model based on risk score, gender, stage, and age in 457 TCGA melanoma patients. We obtained the weights of the immune-related genes (four categories) from the previous study, including major histocompatibility complex (MHC)-related molecules, suppressor cells, effector cells, and checkpoints or immunomodulators. The weighted averaged Z-score was calculated based on the expression. Based on the sum of the weighted average Z-scores of four categories, the IPS was calculated (ranging from 0 to 10) (Charoentong et al., 2017b).

R 3.6.2 was used for data analysis and graph drawing. The detailed methods, software, algorithms, and sources of each cohort are summarized in Table 2. The expression correlations between 21 m6A regulators were analyzed by Spearman correlation coefficient analysis. Univariate Cox analysis was used to screen for genes that were significantly associated with prognosis. Then, a risk signature was constructed via the multivariable Cox regression model based on the stepwise regression with the minimum Akaike information criteria (AIC). A Wilcoxon rank-sum test was used to analyze the differences in continuous variables between the two groups. The log-rank test was used to compare the difference in survival rates between the two groups. Differences were established by chi-square tests and Fisher’s exact test. p < 0.05 was considered to be statistically significant.

A schematic of the process is presented in Figure 1. The landscape of the expression correlation of these m6A regulators was shown in Figure 2A. We found most m6A regulators presented a significant positive correlation in expression, except for ALKBH5. To better understand the role of m6A regulators in the pathogenesis and progression of melanoma, we analyzed the differences in the expression of 21 m6A regulators between normal and tumor samples via the R package limma. Among the 21 m6A regulators, 11 were significantly more highly expressed in the tumor tissues than in the normal tissues (RBM15B, RBM15, ALKBH5, YTHDF1-3, IGF2BP1, CBLL1, ELAVL1, LRPPRC, and ZC3H13; p < 0.05, Figure 2B). Then, we analyzed the somatic mutations of 21 m6A regulators in TCGA melanoma; 19 m6A regulators displayed mutations, with the exception of only two genes (METTL14 and ALKBH5) in melanoma tissues (Figure 2C). Among the 19 m6A regulators showing somatic mutations, IGF2BP1, KIAA1429 and YTHDC1 displayed the highest mutation frequency (7.49, 7.06, and 4.28%, respectively) (Figure 2D). To evaluate the potential signaling pathways controlling the mutation patterns of these three genes, we performed GO enrichment analysis in 67 melanoma cell lines using the mutation and RNA-seq expression profiles from the Cancer Cell Line Encyclopedia (CCLE) dataset (Ghandi et al., 2019). We found that the IGF2BP1 and KIAA1429 mutant cell lines exhibited a significant association with the functions related to DNA repair like melanin biosynthetic process and DNA strand elongation involved in DNA replication, as well as DNA replication and Mismatch repair. While Melanin metabolism pathway, including melanin biosynthetic process, melanosome organization, and Melanogenesis, were remarkably rich in the YTHDC1 mutant cell lines (Supplementary Figure S1).

The somatic copy number alteration (CNV) analysis revealed that most m6A regulators displayed deletion in copy number, except for CBLL1, HNRNPA2B1, IGF2BP1, KIAA1429, LRPPRC, and YTHDF1-3 (Figure 2E). The position of the CNV variation in the chromosome of these m6A regulators is shown in Figure 2F. The high expression of eight m6A regulators (CBLL1, FMR1, HNRNPA2B1, HNRNPC, LRPPRC, RBM15, RBM15B, ZC3H13) was enriched in driver genes-mutated samples. The low expression of one m6A regulator (ALKBH5) was enriched in driver genes-mutated samples (permutation test, p < 0.005, Figure 3). These results showed the high heterogeneities of m6A regulators and indicated that the dysregulation of m6A regulators had a critical role in the pathogenesis and progression of melanoma.

To comprehensively analyze the role of multiple m6A regulators in melanoma, we performed unsupervised clustering of the tumor samples using the expression values of 21 m6A regulators via the R package ConsensusClusterPlus. The consensus clustering algorithm was used to determine the most stable consensus clusters. According to the consensus matrix, when k = 2, the number of patients was evenly distributed in each cluster, without any cluster containing an exceptionally high or low number of patients. The clusters showed a low correlation. When K = 2, the CDF curve is flat, revealing the correct K (Wilkerson and Hayes, 2010; Senbabaoglu et al., 2014; Hu et al., 2021). The maximum Calinski-Harabasz (CH) index indicated that K = 2 was the optimal number of clusters (Supplementary Figure S2). Principal component analysis indicated a significant difference in transcriptome between the two clusters (Figure 4A). Two optimum clusters of m6A were obtained, named m6A-clusterA and m6A-clusterB, respectively (Figure 4B). Heatmap demonstrates the clinicopathological manifestations between the two m6A-clusters (Figure 4C). The association between the m6A clusters and the clinical indexes was analyzed, and we found that the low expression of FMR1, METTL3 and YTHDC2 was associated with the m6A-clusterA (Supplementary Table S1). To further understand the differences in biological processes between the two clusters, we employed GSVA enrichment analysis. The results indicated that the m6A-clusterA exhibited a significant association with the pathways related to immunological responses, such as antigen processing and presentation, natural killer cell-mediated cytotoxicity, chemokine signaling pathway, complement and coagulation cascades, and so on (Figure 4D). To verify the findings in the TCGA melanoma dataset, we performed external validation using a combined melanoma cohort to validate the unsupervised clustering model (n = 520). The result is consistent with the model based on TCGA, and K = 2 was the optimal value (Supplementary Figure S3).

The correlation between the m6A RNA methylation regulators and immune infiltration was further evaluated. To assess the different immunological responses between the two m6A clusters, we measured the immune score, tumor purity, stromal score, mRNAsi, and infiltrating abundances of 28 immune cells using the Wilcoxon rank-sum test. The mRNA expression-based stemness index (mRNAsi) was obtained from the study of Malta et al. It was calculated with an innovative one-class logistic regression machine-learning algorithm (OCLR) based on mRNA expression, giving values between 0 and 1 (Malta TM. et al., 2018). The results showed that the m6A-clusterA had a higher immune score (Figure 5A) and stromal score (Figure 5B) but lower tumor purity (Figure 5C) and mRNAsi (Figure 5D) compared with the m6A-clusterB (Wilcoxon rank-sum test, p < 0.05). The infiltrating abundances of 28 immune cells indicated that most immune cells were enriched in the m6A-clusterA compared with the m6A-clusterB except for the Activated CD4 T cell, Effector memory CD4 T cell, Memory B cell, Type 2 T helper cell, Eosinophil, Gamma delta T cell, Immature dendritic cell, Natural killer cell (Wilcoxon test, p < 0.05, Figure 5E). Similar to previous conclusions, the m6A-clusterA was characterized by a higher immune score and immune cell infiltration in the validation cohort (Supplementary Figure S4).

In addition, we also utilized CIBERSORT, MCPcounter, and TIMER algorithms to evaluate the abundances of tumor immune infiltration. The heatmap was used to show the landscape of immune infiltration based on these five algorithms (Figure 5F). These results showed that multiple m6A regulators-mediated modification patterns might involve regulating tumor immune infiltration of melanoma.

To predict the prognosis and guide individualized treatment, we used the differentially expressed genes with prognostic significance between the two m6A clusters to construct an m6A-related signature. We selected differential expression genes between two m6A clusters (n = 849, p < 0.05 and | log fold change | >0.5) (Supplementary Table S2). Then we performed GO enrichment analysis for these differential expression genes. These differential expression genes were enriched in m6A-related pathways, such as RNA splicing, regulation of mRNA processing, and so on (Supplementary Figure S5). So we considered these differential expression genes as m6A phenotype-related genes, which can be used to construct the m6A-related signature (Zhang et al., 2020b). Then we screened differential genes with prognosis significance via univariate Cox regression analyses (n = 31, p < 0.0001) (Supplementary Table S3). Moreover, 12 genes were used to construct the optimized risk signature with minimum Akaike information criterion (AIC) value via multivariate Cox regression analysis (IL6ST, MBNL1, NXT2, EIF2A, CSGALNACT1, C11orf58, CD14, SPI1, NCCRP1, BOK, CD74, PAEP) (Sun et al., 2020). The risk scores were derived from the following equations: Risk scores = ∑ (coefi × Xi), wherein coefi is the coefficient and xi is the relative expression value, transformed by Z-score, of the 12 genes (Supplementary Table S4).

The patients in the cohort were divided into high- and low-risk groups using the surv_cutpoint function. The Kaplan Meier survival curves showed significant differences in overall survival between the two risk groups. The above results were independently validated in three melanoma cohorts from the gene expression omnibus database (GEO) (log-rank test, p < 0.05, Figures 6A–D). We plotted receiver operating characteristic curves (ROC curves) using timeROC in the R package. We evaluated the prediction accuracy of the risk signature according to the area under these ROC curves (AUC); high AUC values indicate good sensitivity and specificity (Supplementary Figure S6). The association between risk signature and survival remained statistically significant after age, gender, and stage were taken into account [TCGA, HR: 2.79 (95% CI: 2.01–3.85), p < 0.001, Figure 6A; GSE65904, HR: 2.20 (95% CI: 1.43–3.36), p < 0.001, Supplementary Figure S7A; GSE54467, HR:2.70 (95% CI: 1.41–5.20), p = 0.003, Supplementary Figure S7B; GSE78220, HR: 10.49 (95% CI: 1.30–84.15), p = 0.027, Supplementary Figure S7C; TCGA + GEO: HR: 2.41 (95% CI: 1.95–2.97), p < 0.001, Figure 6E]. We calculated the integrated discrimination improvement (IDI) and the net reclassification improvement (NRI) for the comparation of the old model (AJCC stage only, without risk score) and combined model (combined stage and risk score). The IDI for 1- and 3-years were 5.5% (p < 0.001) and 16.5% (p < 0.001), respectively. Consistently, the NRI for 1- and 3-years were 30.2% (p = 0.006), 40.1% (p < 0.001), respectively, (Supplementary Figures S8A,B). From the ROC curves, the AUCs of combined model for the survival prediction in 1- and 3-years were 0.721 and 0.745, respectively, higher than the old model (Supplementary Figures S8C,D). In addition, we developed a nomogram based on risk score, age, gender, and stage for survival rate prediction (Figure 6F). The C-index of the nomogram was 0.726 (95% CI, 0.685–0.767), markedly higher than AJCC stage (0.621, 95% CI: 0.580–0.662). A calibration curve at 3 or 5 years showed high consistency between predicted survival probability and actual survival proportion (Figure 6G).

The low-risk group had higher PDL1 (Figure 7A), PD1 (Figure 7B) and CTLA4 expression (Figure 7C), and immunophenoscore (IPS) (Figures 7E,F) compared with the high-risk group, indicating that the low-risk group may exhibit a better response to immunotherapy. The results were validated in GEO-combined cohorts (Supplementary Figure S9). The correlation between the 28 immune infiltration cells and the risk scores was examined using Spearman correlation analysis (Figure 7G).

Notably, nearly all the writers and readers were significantly upregulated in the low-risk group, indicating the potential correlation of m6A levels and immunotherapy response in melanoma (Figure 7H).

Furthermore, the GO enrichment analysis revealed that the low-risk group was enriched in pathways associated with immune fully activation, including MHC class II protein complex, regulation of dendritic cell apoptotic process, and T cell receptor complex (Figure 7I). In comparison, the high-risk group was prominently related to Melanin metabolism pathways such as melanin biosynthetic process, tyrosine metabolic process, and melanocyte differentiation (Figure 7J). In addition, the risk scores of the BRAF wild-type group were higher than that of the mutation group (Figure 7D).

To further test the predictive efficiency of m6A-related signature in immunotherapy cohorts, we analyzed the proportion of patients with response to immune checkpoint blockade therapy in low- and high-risk groups. The immunotherapy cohorts included a combined melanoma cohort (ICB-therapy-combined melanoma) and urothelial cancer patients treated with anti-PDL1 (ICB-therapy-UC). As regards ICB-therapy-combined melanoma, the patients in the low-risk group had a significantly improved survival rate (log-rank test, p < 0.05, Figure 8A) and a higher response rate to ICB therapy (Fisher’s exact test, p < 0.05, Figure 8B) compared with the high-risk group. Consistent with ICB-therapy-combined melanoma, the patients in the low-risk group exhibited a significantly prolonged survival rate (log-rank test, p < 0.05, Figure 8C) and a higher response rate to immune-checkpoint blockade therapy (Fisher’s exact test, p < 0.05, Figure 8D) compared with the high-risk group in ICB-therapy-UC.