Identification and Verification of Molecular Subtypes with Enhanced Immune Infiltration Based on m6A Regulators in Cutaneous Melanoma

Background As the most aggressive type of skin cancer, cutaneous melanoma (CM) is experiencing a rapidly rising mortality in recent years. Exploring potential prognostic biomarkers or mechanisms of disease progression therefore has a great significance for CM. The purpose of this study was to identify genetic markers and prognostic performance of N6-methyladenosine (m6A) regulators in CM. Method Gene expression profiles, copy number variation (CNV), and single nucleotide polymorphism (SNP) data of patients were obtained from The Cancer Genome Atlas (TCGA) database. Results Genomic variation and association analysis of gene expressions revealed a high degree of genomic variation in the presence of m6A-regulated genes. m6A patients with high-frequency genomic variants in the regulatory gene tended to develop a worse prognosis (p < 0.01). Unsupervised cluster analysis of the expression profiles of m6A-regulated genes identified three clinically distinct molecular subtypes, including degradation-enhanced subgroup and immune-enhanced subgroup, with significant prognostic differences (p = 0.046). A novel prognostic signature, which was established according to m6A-related characteristic genes identified through genome-wide expression spectrum, could effectively identify samples with poor prognosis and enhanced immune infiltration, and the effectiveness was also verified in the dataset of the chip. Conclusion We identified genetic changes in the m6A regulatory gene in CM and related survival outcomes. The findings of this study provide new insights into the epigenetic understanding of m6A in CM.


Introduction
As the most aggressive type of skin cancer, cutaneous melanoma (CM) resulted in 287,723 new cases and 60,712 deaths worldwide in 2018 [1]. Although it accounts for only 4% of all skin cancer cases, it is far more dangerous than other skin cancers. The treatment and control of CM still remain less effective, though great efforts have been devoted to the improvements in managing advanced regional and metastatic CM [2]. UV exposure is known to be closely associated with the development of melanoma [3]. The incidence of CM shows regional variations such as ethnic skin phenotypes and differences in sun exposure [4]. A study demonstrated that the mortality of melanoma patients is rapidly increasing in recent years [5] and also pointed out the management of melanoma patients remains a complex and evolving subject [6]. Thus, exploring potential prognostic biomarkers or the mechanism of disease development is highly necessary for CM.
RNA methylation has become a common phenomenon and a key regulator of transcript expressions [7]. N6-Methyladenosine (m6A), which is methylated at the N6 position of adenosine, is considered to be the most common conserved internal transcriptional modifications of messenger RNAs (mRNAs), microRNA (miRNAs), and long noncoding RNAs (lncRNAs) [8]. The deposition of m6A is encoded by a methyltransferase complex involving three homologous factors, including methyltransferases (termed as "writers"), demethylases (termed as "erasers"), and recognition from m6A binding proteins (termed as "readers"). The m6A dysregulation consists of multiple cellular processes even in human cancers. It has been demonstrated that m6A mRNA demethylase FTO could regulate melanoma tumorigenicity and react to anti-PD-1 blockade [9]. The m6A methyltransferase METTL3 promotes osteosarcoma progression by regulating the m6A level of LEF1 [10]. Furthermore, METTL3 (an oncogene) maintains SOX2 expression via the m6A-IGF2BP2-dependent mechanism in colorectal carcinoma cells and therefore could serve as a potential biomarker for cancer prognosis prediction [11].
So far, however, little is known about the relationship between m6A-related genes and CM. The Cancer Genome Atlas (TCGA) database allows an easy access to human cancer data of gene expressions, copy number variation (CNVs), and single nucleotide polymorphism (SNPs), which all play important roles in the development and progression of human cancers [12]. However, the CNV and SNP of m6Arelated genes remain unknown to us.
This retrospective study developed a novel gene signature based on m6A regulators in CM with data acquired from TCGA database and also analyzed the performance of the signature. By analyzing genomic variations and gene expression associations, we found that these m6A regulatory genes showed high genomic variations; interestingly, patients with m6A regulatory genes of high-frequency genomic variation tended to develop a worse prognosis (p < 0:01). Three clinical subtypes with different molecular characteristics and significant prognostic differences were identified by performing unsupervised clustering analysis on the expression profiles of m6A regulatory genes. m6A-related characteristic genes were determined based on genome-wide expression profiles. A new prognostic signature was then built for identifying subgroups with poor prognosis and enhanced immune infiltration, and its performance was validated in the training set and validation set. In summary, we determined genetic changes in the m6A regulatory genes in CM and patient survival outcomes. The findings of this study provide new insights into the epigenetic understanding of m6A in CM.

Materials and Methods
2.1. Data Resource and Processing. All CM clinical data, copy number variation (CNV), single nucleotide polymorphism (SNP), and mRNA expression data were retrieved from TCGA-Assembler of the TCGA website [13] and downloaded in May 2019. For transcriptome data, we obtained 472 samples and downloaded data of read counts. Data was normalized by R package DESeq. For SNP data, we obtained a total of 469 samples, and the downloaded data was processed by MuTect [14]. For CNV data, there were 940 samples acquired by R package RTCGA. Here, the "deletion" was defined as segment_mean < −0:2, and the "amplification" was defined as the segment_mean > 0:2. For clinical information data, initially, there were a total of 389 clinical samples; after integrating the data and excluding samples with a survival time shorter than 90 days, a total of 250 CM samples were finally decided for further analysis.
The expression spectrum dataset GSE65904 [15] with 214 melanoma samples from Illumina HumanHT-12 V4.0 Expression BeadChip platform was downloaded from the Gene Expression Omnibus (GEO) database. The samples with a survival time shorter than 90 days were selected, resulting in a final of 189 samples. Multiple probes corresponding to a gene were retained and shown as the median of the gene expression level, while probes corresponding to multiple genes were eliminated. Finally, the expression profile of genes was obtained. The work flow chart is shown in Figure 1.

Univariate Cox Survival Analysis.
Here, the prognostic performance of m6A regulatory genes was examined. Based on the expression profile of m6A regulatory genes and clinical follow-up data of the samples, each m6A regulatory gene was evaluated using R software package survival, and the forest map was plotted applying R software package forestplot.

Clustering Analysis.
For a better classification of the patients, an unsupervised clustering method was employed in hierarchical clustering analysis of the expression profile of CNV-related m6A regulatory genes. The clinical subtypes were selected according to the minimum intragroup variance and the maximum intergroup variance of the method.

Differential Expression Gene
Analysis. Gene differences between each clinical subgroup were analyzed using R software package DESeq2, with |log 2ðfold changeÞ | >2 and FDR < 0:05 serving as the thresholds.
2.5. Analysis of Immune Cells and Immune Infiltration. The ssGSEA (single-sample gene set enrichment analysis) algorithm was applied to quantify the relative abundance of cell invasion in each tumor microenvironment (TME) of a CM patient. The gene sets labeled as TME-infiltrating immune cell type were obtained from Charoentong et al.'s study [16], which investigated a variety of human immune cell subtypes, including activated CD8 T cells, activated dendritic cells, macrophages, natural killer T cells, and regulatory T cells. The enrichment fraction calculated by ssGSEA was used to represent the relative abundance of TMEinfiltrating cells in the samples. Patient's infiltration score was assessed by the R package Estimate.

KEGG Pathway Enrichment Analysis of Different Clinical
Subgroups. ssGSEA KEGG pathway analysis was performed on each sample using R software package GSVA [17] to calculate the difference of enrichment score of each pathway in different clinical subgroups. Moreover, the relationship between enrichment scores of each sample pathway and clinical subgroup was analyzed by performing Pearson correlation. p < 0:05 was the threshold value in determining the KEGG pathway that was the most relevant to the clinical subgroup. BioMed Research International model was established to calculate the risk score of each sample according to the model. Samples scored above 75% were assigned into the high-risk group with immunorejection phenotypes, while those scored below 75% were in the low-risk group.
2.8. Statistical Analysis. All statistical data were analyzed by SPSS 23 (IBM, Chicago, USA) and R language. The association between CNV and SNP of m6A regulatory genes and clinicopathological characteristics was examined by the chisquared test. The association between three CM key genes and CNV and SNP of m6A regulatory genes was analyzed by the chi-squared test. The Kaplan-Meier curve and logrank test were applied to evaluate the prognostic performance of the alterations in m6A regulatory genes. All statistical results with a p value ≤ 0.05 were considered to be significant.

Results
3.1. The m6A Regulatory Genes Showed a High Frequency of Genomic Variation in CM Patients. A total of 17 m6A regulatory genes were recruited in this study. 136 out of the 242 SNP of m6A regulatory genes with sequencing data appeared as independent samples (Table 1). Among them, the "writer" gene KIAA1429 in 44 samples and the "reader" gene IGF2BP1 in 41 samples demonstrated a higher SNP, with all the "reader" genes showing a higher SNP frequency than that of the "writer" or "eraser" genes; noticeably, the frequency of SNP of the "eraser" gene was the lowest (Figure 2(a)). We also observed that in all 469 CM samples with CNV data, the m6A regulatory genes had a high frequency of CNV ( Figure 2(b)). For example, the "writer" gene WTAP showed the highest frequency of 57.08% of its CNV events, and the frequency of the "reader" gene IGF2BP3 was 52.1% (Table 2). These results indicated a high frequency of genomic SNP in the m6A regulatory gene in CM, and these abnormalities may affect gene expression and lead to different clinical outcomes in CM patients.

CNV of m6A Regulatory Genes Was Associated with
Adverse Clinical Outcomes. Considering that changes in CNV can affect gene expression levels through dosage compensation effects, to this end, the effect of m6A-regulated gene changes on mRNA expression was evaluated. In 472 KM samples, mRNA expression levels showed a significant correlation with CNV patterns. For all 17 regulatory genes, except for IGF2BP3 and IGF2BP1, increased CNV of the rest 15 genes was related to a higher mRNA expression and deletion in mRNAs with decreased expression (Figure 3(a)). Furthermore, the prognostic differences between samples of m6A regulatory genes with the high-frequency and lowfrequency CNV abnormalities were analyzed. We observed that patients with a high CNV tended to show worse clinical outcomes (p < 0:001), directly indicating that abnormal copy number of m6A regulatory genes was related to poor clinical outcomes ( Figure 3(b)). Finally, by analyzing the relationship between the expressions of m6A regulatory genes and patients' prognosis, we observed that RBM15, YTHDF1, WTAP, and METTL14 genes were significantly associated with CM prognosis (Figure 3(c)). Specifically, highexpressed YTHDF1 was indicative of a poor prognosis and could therefore be considered a risk factor, at the same time, low-expressed RBM15, WTAP, and METTL14 were also related to a poor prognosis and were all seen as protective factors. In order to further determine the prognostic value of m6A regulators, a risk score model (m6AScore) was established based on a multivariate regression method to calculate the risk score of each patient in TCGA cohort, and the prognostic value of m6AScore was determined by univariate and multivariate Cox survival analyses. The results of univariate analysis demonstrated that m6AScore and a variety of clinical features such as T, N, and stages were strongly correlated with CM patients' prognosis ( Figure 3(d)), while multivariate

BioMed Research International
analysis showed that only m6AScore, N3, N2, and T4 had a significant prognostic correlation (Figure 3(e)). In addition, HR of N3, N2, and T4 was abnormally high. These results indicated that m6AScore has a more stable prognostic performance than T, N, M, and stage or other clinical features.

Different Molecular Subgroups Identified by Unsupervised
Cluster Analysis Based on Expressions of m6A Genes. Given that m6A regulatory genes are associated with CM prognosis, unsupervised hierarchical cluster analysis was conducted based on the expression profiles of 15 CNV-related m6A regulatory genes. We observed that these 15 genes were divided into three main categories (Figure 4(a)) with different expression patterns of m6A genes. Further analysis showed that 12 genes (80%) out of the 15 genes in the three subtypes had significant expression differences (Figure 4(b)). Subsequently, prognostic differences among the three types of samples were examined, and the data demonstrated that m6A Cluster3 was related to a significantly better prognosis than m6A Cluster1 and m6A Cluster2 (Figures 4(c) and 4(d)). These results suggest that the three m6A Clusters have different molecular characteristics, which will lead to different clinical outcomes.
3.4. The Regulatory Characteristics of m6A Cluster. To explore the relationship between the biological behaviors of CM and m6A Cluster, the enrichment scores of each sample in the KEGG pathway were calculated by ssGSEA to further determine the differences of the enrichment scores ( Figure 5(a)). Significant differences in 44 KEGG pathways can be observed; interestingly, 40 pathways showed significant positive correlations with m6A Cluster2, and they were mainly important signaling pathways involved in tumor development, for example, VEGF_signaling_pathway, basal_cell_carcinoma, and drug_metabolism_cytochrome.
In addition, the remaining 4 KEGG pathways, namely, ubiq-uitin_mediated_proteolysis, non_homologous_end_joining, basal_transcription_factors, and RNA_degradation, were significantly upregulated in m6A Cluster3 and are mainly related to the degradation process. 28 immune cell components were calculated for each patient by the ssGSEA, and a significant difference in the distribution of 19 immune cells (67.9%) was detected ( Figure 5(b)). Among these significant immune cells, memory B cell and type 2 T helper cell showed the highest scores in m6A Cluster1, while the remaining 17 immune cells all had the highest scores in m6A Cluster2. These results suggested that m6A Cluster3 is related to immune processes and a better CM prognosis, while m6A Cluster2 is related to a stable immune microenvironment and a worse prognosis. We also compared the relationships between the 3 m6A-related molecular subtypes and the previously reported genomic subtypes (BRAF, RAS, NF1, and triple-WT) [18]). The data showed that m6A Cluster3 had the largest intersection with the RAS subtype, whereas m6A Cluster1 and m6A Cluster2 had the largest intersection with triple_WT ( Figure 5(c)). In our study, the results that m6A Cluster1 and m6A Cluster2 were indicative of a poorer prognosis and higher frequency of m6A copy number are consistent with triple-WT with significantly more copy number fragments.

Identification of m6A Gene Markers and Molecular
Characteristics. To determine the gene expression differences among the m6A cluster samples, we identified 5,533 differentially expressed genes based on the gene expression profile. It has been found that m6A Cluster2 and m6A Cluste1/m6A Cluster3 had the most differentially expressed genes, as compared with m6A Cluster1 (Figure 6(a)). Gene expression profiles of these 5533 genes were extracted for cluster analysis,      (Figure 6(b)). Further analysis demonstrated that Cluster2 was related to a significantly worse prognosis than Cluster1 (Figure 6(c)). Moreover, we also detected the expression differences of 15 m6A regulatory genes with CNV abnormalities in Cluster1 and Cluster2 ( Figure 6(d)), and 12 (80%) out of the 15 genes showed significant expression differences, and among the 12 genes, ALKBH5 and YTHDF1 were significantly highly expressed in Cluster2, while the remaining 10 genes were significantly highly expressed in Cluster1. These results indicated that a class of molecular subgroups (m6A Cluster1/Cluster2) with unfavorable prognosis which resulted from abnormalities of immune microenvironment could be identified by m6A regulatory genes or m6Arelated gene markers.

Establishment and Verification of Prognostic Signature
Based on m6A Regulatory Genes. The subgroup of immunerelated prognosis was determined by performing the principal component analysis based on the expression profiles of m6A-related gene markers to divide the samples into a high-risk group and low-risk group. A careful analysis of the prognostic differences between the two groups showed that the prognosis in the high-risk group was significantly worse than the low-risk group (Figure 7(a)). Moreover, a higher immune score was observed in the high-risk group, as compared with the low-risk group (Figure 7(b)). Further-more, the patients were accordingly group based on the chip platform GSE65904. Here, the data showed that patients' prognosis in the high-risk group sample was greatly more unfavorable than that in the low-risk group (Figure 7(c)). Immunomicroenvironment analysis demonstrated that the samples from the high-risk group had higher immune microenvironment scores (Figure 7(d)). In addition, the Swegene Center for Integrative Biology at Lund University melanoma cohort (GSE22153) was introduced for verification. It has been found that the prognosis of the high-risk group was significantly worse than that of the low-risk group samples (Figure 7(e)). Immunomicroenvironment analysis showed that the high-risk group samples had a higher immune microenvironment score (Figure 7(f)). These results suggested that the expression profiles of m6A-related gene markers could serve as prognostic markers to evaluate the prognosis of CM patients.

Discussion
The present study developed a novel gene signature based on m6A regulators to CM using TCGA database and also assessed its predictive performance. Genetic alterations of the m6A regulatory genes have been found to be closely related to CM patients' survival outcomes. Our findings provide new insights into the epigenetic understanding of m6A in CM.
As the most common internal chemical modification of mRNA, m6A is widely involved in many pathological processes of cancer development. The deposition of m6A is encoded by a methyltransferase complex involving three homologous factors including methyltransferases (such as

12
BioMed Research International METTL3/14, WTAP, RBM15/15B, and KIAA1429, termed as "writers"), demethylases (such as FTO and ALKBH5, termed as "erasers") and recognition from m6A binding pro-teins (such as YTHDF1/2/3, IGF2BP1 and HNRNPA2B1, termed as "readers") [8].       16 BioMed Research International mRNA level of lymphoid enhancer binding factor 1 (LEF1), thereby inhibiting the activity of the Wnt/β-catenin signaling pathway. Moreover, METTL3 promotes osteosarcoma cell progression by regulating the m6A level of LEF1 and activating the Wnt/β-catenin signaling pathway [10]. Previously, upregulated expression of METTL3 has been observed in human hepatocellular carcinoma and multiple solid tumors [19]. Knockdown of METTL3 can reduce hepatocellular carcinoma cell proliferation, migration, and colony formation, and more interestingly, knockout of METTL3 remarkably suppresses tumorigenicity and lung metastasis. In colorectal cancer, METTL3 acts as a tumor suppressor on cell proliferation, migration, and invasion via the p38/ERK pathway [20]. In addition, reducing m6A methylation could activate oncogenic Wnt/PI3K-Akt signaling and promote malignant phenotypes of gastric cancer cells [21]; moreover, METTL3 knockdown could inhibit the level of total RNA methylation of m6A, cell proliferation, and migration of gastric cancer cells [22]. Noticeably, m6A methylation has also been reported in bladder cancer [23], renal cancer [24], acute myeloid leukemia [25], breast cancer [26], and hepatocellular cancer [27]. The regulators of m6A RNA methylation function critically in the malignant progression of glioma and may be use-ful for the development of prognostic stratification and therapeutic strategies [28]. In this study, 15 m6A regulatory genes with abnormal copies of the genome were identified based on TCGA gene expression data and CNV data. Our analysis showed that a higher frequency of the CNV of m6A regulatory genes was related to a worse the prognosis of CM patients. Three clinical subgroups with different molecular characteristics were detected after performing unsupervised cluster analysis. Here, m6aCluster2 was associated with abnormal immune microenvironment and poor prognosis. A risk assessment model was established with the m6A-related gene markers, and subgroups with high immune scores and poor prognosis were identified in both sequencing data and chip data. Our results indicated that m6A regulatory genes play an important role in CM, as they could influence immune-related processes during tumor development. Zhou et al. [29] proved that a low expression of the writer gene METTL3 is related to activations of adipogenesis and mTOR pathways. However, m6A regulators in CM have not been reported before.
Some cancer-related pathways are dysregulated during CM development. In this study, We determined a molecular subgroup with favorable prognosis based on m6A regulatory genes and found that degradation pathways such as ubiquitin  Figure 7: Establishment and validation of risk scores based on m6A-related genetic markers: (a) KM curves for prognostic differences between high-risk and low-risk groups; (b) differences in immune microenvironment scores between high-risk and low-risk groups; (c) KM curves for prognostic differences between high-risk and low-risk groups in the GSE65904 dataset; (d) differences in immune microenvironment scores between high-risk and low-risk groups in the GSE65904 dataset; (e) KM curves for prognostic differences between high-risk and low-risk groups in the GSE22153 dataset; (f) differences in immune microenvironment scores between high-risk and low-risk groups in the GSE22153 dataset.
pathways, mediated_proteolysis, and RNA_degradation were active in the subgroups, and interestingly, these pathways play critical roles in regulating cell growth and proliferation by controlling the abundance of key cyclins [30]. Evidences increasingly showed that abnormal proteolysis of cell cycle regulators significantly promotes tumorigenesis, while enhanced degradation of cell cycle regulators (i.e., protooncoproteins) can inhibit tumor metastasis. The driving force of the cell cycle is the activation of cyclin-dependent kinases (CDKs), whose activity is controlled by the proteolysis of ubiquitin-mediated key regulators (such as cyclins and CDK inhibitors), and disorder of proteolytic system may lead to uncontrolled proliferation, genomic instability, and cancer initiation [31]. Selective RNA degradation has been shown to be able to strongly resist disease development and is therefore believed to have the potential to treat cancer [32]. In antioral tumor therapy, activation of peroxisomal proliferator activator (PPAR) upregulates RNA degradation pathways, leading to reprogramming of tumor metabolism [33]. These results suggest that degradation-enhancing molecular subgroup identified by the m6A regulatory genes could be explored as a novel effective strategy for CM treatment.

Conclusion
In summary, this study first detected the genetic alterations of m6A regulatory genes in CM and identified molecular subsets with enhanced degradation and immune enhancement based on genomic mutagenesis of m6A regulatory genes. A m6A signature with a high predictive accuracy in predicting the prognosis of CM patients was constructed in this study. Our findings provide new insights into the epigenetic understanding of m6A in CM.

Data Availability
The datasets analyzed in this study are available from the GEO repository (https://www.ncbi.nlm.nih.gov/geo/).