Increased expression of YTHDF1 and HNRNPA2B1 as potent biomarkers for melanoma: a systematic analysis

The incidence and mortality of melanoma is increasing around the world. To deeply explain the mechanism insight into it, we conducted a systematic analysis to examine the levels of regulatory genes of the common RNA epigenetic modification-N6-methyladenosine (m6A) in patients with melanoma compared by the healthy. We analyzed the expression of m6A Eraser, Writer, and Reader genes based on publicly available datasets on Oncomine and validated the results with a gene expression omnibus dataset. Hub genes were identified with Cytohubba and the frequency of copy number alterations was analyzed with the cBioPortal tool. The results revealed the up-regulation of YTHDF1 and HNRNPA2B1 in melanoma. Combining the two genes improved the efficacy in diagnosing melanoma by about 10% compared to each gene alone. Hub genes identified with four analysis methods were compared and the overlapping genes were selected. These genes were enriched in several gene ontology terms. Genes related to p53-signaling consisted of CDK2, CDK1, RRM2, CCNB1, and CHEK1. All five genes were positively correlated with either YTHDF1 or HNRNPA2B1, suggesting that both genes may affect m6A modification by the five genes, further up-regulating their expression and facilitate their roles in inhibiting p53 to suppress tumorigenesis. We also observed major mutations in YTHDF1 and HNRNPA2B1 that led to their amplification in melanoma. Significant differences were observed in the clinical characteristics of patients with altered and unaltered m6A regulatory genes such as tumor stage and treatment response. We, for the first time, identified a combination of m6A regulatory genes to diagnose melanoma. We also analyzed m6A-related genes more comprehensively based on systematic complete data. We found that YTHDF1 and HNRNPA2B1 were altered in melanoma and might influence the development of the disease through signaling pathways such as p53.


Background
Melanoma is one of the fastest developing malignancies with strong aggressive ability, however no proper curative treatments exist at present [1,2]. It is a type of epithelial malignant tumor originating from melanocytes with obviously increasing incidence and mortality around the world [1]. Based on insight into the molecular mechanism underlying the disease, serum lactate dehydrogenase is considered as a biomarker to predict the development of melanoma [3].
With the advancement of sequencing technology, related mutations such as B-Raf proto-oncogene (BRAF) and signaling pathways such as the mitogen-activated protein kinase (MAPK) pathway have been found. These discoveries have led to the emergence of targeted drug treatment for melanoma such as inhibitors of BRAF, mitogen-activated protein kinase kinase 7 (MAP2K7, also known as MEK), and mitogen-activated protein kinase 1 (MAPK1, also known as ERK) [2,4,5]. However, the actual pathological mechanism in melanoma remains unknown.
N6-methyladenosine (m 6 A) is a type of RNA epigenetic modification [6]. Originally reported only in mRNA, experts have found m 6 A in other types of RNA and involved in various bioprocesses with the development of high-throughput sequencing technology [7]. Similar to DNA methylation modification, m 6 A is also regulated by methyltransferase and demethylase [7,8]. Previous reports showed that m 6 A methylation mainly interacted with three classes of proteins: firstly, m 6 A methyltransferase whose coding genes were called Writers and included methyltransferase-like 3 (METTL3), methyltransferase-like 14 (METTL14), and WT1 associated protein (WTAP); secondly, m 6 A demethylase whose coding genes were called Erasers and included alpha-ketoglutarate dependent dioxygenase (FTO) and alkB homolog 5 (ALKBH5); thirdly, proteins that could attach to the m 6 A methylation site in RNA and further played a specific role in physio-pathological processes with coding genes called Readers, which included YTH N6-methyladenosine RNA-binding protein 1 (YTHDF1), YTH N6-methyladenosine RNA-binding protein 2 (YTHDF2), E74-like ETS transcription factor 3 (ELF3), heterogeneous nuclear ribonucleoprotein C (HNRNPC), and heterogeneous nuclear ribonucleoprotein A2/B1 (HNRNPA2B1) [6,9]. Recent studies have shown that m 6 A methylation is closely related to tumorigenesis and tumor development [10]. Vu et al. [11] reported that CD34+ hematopoietic stem cells with lower expression of METTL3 had higher levels of phosphorylated protein kinase B and that METTL3 promoted normal stem cells into acute myeloid leukemia cells. Zhou et al. [12] reported increased expression of FTO in tissue lesions of cervical squamous cell carcinoma (CSCC) and observed that patients with high FTO expression exhibited chemotherapy tolerance. However, few studies have investigated the role of m 6 A modification in the development or pathological process of melanoma.
In this research, our goal was to examine the expression of m 6 A regulatory genes in melanoma based on open biological data and to identify the related pathways. Firstly, we used the online Oncomine tool to determine whether the m 6 A-related genes presented above were altered in melanoma. We then extracted the RNA-seq data from the gene expression omnibus (GEO) website (http://www.ncbi.nlm.nih.gov/geo) to validate the results from Oncomine and investigate the altered genes. Next, we used statistical methods to evaluate the hub genes involved in RNA processing and calculate the correlation between the hub genes and the selected target genes. Finally, we analyzed all the critical genes and identified the signaling pathways and mutation hot spots for m 6 A regulatory genes to uncover new therapy targets and determine the m 6 A-related biological mechanism in melanoma. A summary of this study is shown in Fig. 1.

Ethics statement
All clinical data, copy number alterations (CNAs), mutations, RNA-seq data, and other data were abstracted from Oncomine, cBioPortal (http://cbiop ortal .org), or the GEO platform, which are open to the public. All analyses were performed based on the data of previously published studies. Therefore, no ethical approval or patient consent was required for this study.

Overview of the expression of 10 genes in Oncomine
The Oncomine (https ://www.oncom ine.org) tool was used to determine the expression of m 6 Table S1. We used "Differential analysis" online tool in Oncomine portal to systematically compare the DNA or transcript levels of the normal and patients with melanoma. The whole data in the five studies above were corrected and then median ranked analysis was performed, results were shown in Fig. 2.

Validation of the expression of all genes by RNA-seq data from GEO
We obtained the RNA-seq data from GEO (23 normal vs 57 patients with melanoma). Patient information is shown in Additional file 1: Table S2. R package affyPLM and RColorBrewer were used to do the quality control. Robust Multi-Array Average method was used in data pre-processing, then the files from the normal and patients' samples were merged. Package impute in R was used to normalize and correct the files and package limma was to determine the different expression files. Heat map and volcano graphs were drew by gplots package in R software. The expression files of target genes were abstracted to do further analysis.

Functional enrichment analysis of the normalized RNA expression files from GEO
Gene ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes analysis were performed by ClueGO [13,14]. Biological Process, Cellular Component, Molecular Function, Immune System Process ontologies and KEGG pathways were selected. GO tree interval was from 3 to 8, the minimal number of genes in GO term or pathway was 3. GO term/pathway network connectivity (Kappa Score) was 0.4, the statistical option was Enrichment/Depletion (Two-sided hypergeometric test) and pV correction was Bonferroni step down. The String tool (https ://strin g-db.org/) was used to draw the general network. Cytoscape (CytoHubba) was used to predict and explore the critical genes, nodes' score was calculated in this software and the scores of maximal clique centrality (MCC), maximum neighborhood component (MNC), edge percolated component (EPC), and degree methods of each genes were obtained to determine the key genes in this study.

Analysis of mutation diagram for each gene with cBioPortal
We used cBioPortal to analyze the frequency of mutations and CNAs in patients with melanoma. A total of 653 samples were included in this research. Detailed information for the patients is shown in Table 1. Oncoprint tool in cBioPortal was used to scan mutation frequency of each gene in every sample, samples with mutation were presented in Fig. 5. There were three studies (Snyder, MEJM; TCGA; Van Allen, Science) included in this research. Cancer Type tool was used to show the mutation types of every gene or total genes in each study. Mutations tool in cBioPortal was used to showed the mutation sites for each gene in melanoma, and the protein post-translational modification (PTM) site plugin in mutations tool was used to predict the occurrence of phosphorylation, acetylation, ubiquitination, methylation, and O-linked glycosylation, primarily between 0 and 370aa.

Statistics
For the comparison of two different groups, the t test or Mann-Whitney test was used according to the data distribution. The Chi squared test was performed to compare the discrete data and the Kruskal-Wallis test was used to compare the expression of more than two groups that did not conform to the conditions of the parametric test. The relationship between two continuous variables was evaluated with the Pearson correlation coefficient. Receiver operating characteristic (ROC) curves and scatter diagrams were drawn using GraphPad Prism 6.0. IBM SPSS Statistics version 21.0 was used to calculate all statistical parameters. The significance level was P < 0.05.

Systematic comparison of m 6 A regulatory genes in melanoma and normal populations
We used the Oncomine database to compare the expression of m 6 A-related genes in patients with melanoma and healthy individuals reported in five studies. The inclusion criteria were as follows: patients involved in the studies had a definite diagnosis or histopathological diagnosis of melanoma; the samples were from patients confirmed to have clinical data; the data were collected from tissue samples not cell lines; gene expression profiles were available. We excluded the datasets if : they were without clinical information; the samples' type was cell line or primary cell culture; their gene expression files were not provided on the Oncomine portal; they did not have samples of both patients with melanoma and the healthy. Samples were included in this study regard less of the age, race or sex. As shown in Fig. 2, YTHDF1 and HNRNPA2B1 were significantly over-expressed in patients with melanoma while ELF3 was down-regulated. The detailed P-values are shown in Table 2. Other genes such as HNRNPC were also altered to some degree. We further validated our findings using normalized RNA-seq data from GEO.

Validation of the expression of target genes
Original RNA-seq files for 53 patients with melanoma and 27 healthy individuals were obtained from GEO. The gene expression for all participants is shown in a cluster heat map in Fig. 3a and the up-regulated (red) and down-regulated (green) genes are visually presented in volcano plots (Fig. 3b). We examined the 10 genes of interest in this dataset and found that only HNRNPA2B1, YTHDF1, ELF3, and YTHDF2 were significantly altered in patients with melanoma (Table 3, Fig. 3c, Additional file 2: Figure S1). Integrating the results from GEO and Oncomine, the expression of HNRNPA2B1, YTHDF1, and ELF3 showed the same patterns and were thus persuasive. We assessed the specificity and sensitivity of these three genes and observed that HNRNPA2B1 and YTHDF1 had higher ROC areas than YTHDF2, whereas ELF3 had a lower value (Table 3). Therefore, we selected HNRNPA2B1 and YTHDF1 for further research and combined the expression of the two genes to evaluate the combined specificity and sensitivity. The equation for the combination used in this study was X = logit (P) = ln (P/1 − P) = − 7.777 + 3.93*YTHDF1 + 1.024*HNRN-PA2B1. Next, the predicted probability of melanoma was calculated with the following equation: predicted probability (P) = e x /(1 + e x ) [15,16]. The results indicated the efficiency of combining YTHDF1 and HNRNPA2B1 compared to the results for each gene individually, with an improvement of almost 10% (Fig. 3c, Table 3). The expression and ROC curves for ELF3, and YTHDF2 are presented in Additional file 2: Figure S1.

m 6 A-associated mechanism in melanoma
We extracted the GO terms for all genes using ClueGO and selected 611 genes involved in RNA processing or the epigenetic activation processes (Fig. 4a). After filtering the unconnected genes from the network using the String tool (Fig. 4b), we chose 586 nodes to determine the important genes using CytoHubba in Cytoscape. The top 50 genes according to the MCC, MNC, EPC and degree were extracted and a Venn diagram was created using Venny tool online. Thirty-one genes overlapped in the results from the four types of analysis (Fig. 4c). The network and MCC score heat nodes are illustrated in Fig. 4d, which shows that cyclin dependent kinase 2 (CDK2), CDK1, cyclin-B1 (CCNB1), and checkpoint kinase 1 (CHEK1) were altered. The detailed ranks for MCC and the other analysis methods are shown in Table 4. Further we performed the correlation analysis of all 31 genes and YTHDF1 or HNRNPA2B1 by calculating their Pearson correlation coefficients (Table 5). Only 17 genes were ether positively or negatively correlated with YTHDF1 or HNRNPA2B1 (P < 0.05) ( Table 5).
Next, we created a network map to visualize the connections between these genes and further explored the pathways behind them (Fig. 4g). YTHDF1 and HNRN-PA2B1 interacted with targeted genes involved in bioprocesses such as the p53 signaling pathway and positive regulation of DNA replication, which might provide insight into the molecular mechanism underlying melanoma (Fig. 4g).

Table 2 P-values of different genes in five studies selected in Oncomine
The

Mutation and CNAs of m 6 A-related genes in patients with melanoma
We analyzed the data for melanoma patients with CNA and mutation profiles on cBioPortal. All mutated samples are shown in Fig. 5a. Integration of the mutation profiles with the clinical information of the patients revealed statistical differences between patients with alterations and those without alterations in terms of the tumor stage, durable clinical benefits, neo-antigen load, serum lactate dehydrogenase, and treatment response ( Table 1). The mutation diagram shows the top three altered genes were YTHDF1 (8%), ELF3 (7%), and HNRNPA2B1 (6%), with a high frequency of amplification. Three studies were included in our systematic analysis. The amplification of genes was dominant in Snyder's study, while mutation and amplification were the major modifications in the other two studies (Fig. 5b). For YTHDF1 and HNRNPA2B1, the predominant alteration was amplification, followed by mutation (Fig. 5b). Other genes are shown in Additional file 2: Figure S2a.

Discussion
Melanoma is a tumor that consists of malignant melanocytes, which are pigment-producing cells of neuroectodermal origin that can be found throughout the body (including in the iris, rectum, and skin) [2]. With the advancement of science and technology coupled with increased recognition of the disease, experts' interpretation of the pathogenesis of melanoma has changed from the previous explanation of simple molecular interactions to epigenetic regulation [17,18]. As one of the most common chemical modifications in messenger RNA (mRNA), m 6 A has received substantial attention in epigenetics recently [10,18]. Last year, Yang et al. [19] found that FTO, an m 6 A demethylase, could promote tumorigenesis and anti-PD-1 resistance in melanoma, highlighting the role of FTO in immunotherapy for the disease. In the same year, Jia et al. [18] observed a decrease in m 6 A levels in ocular melanoma samples, which could indicate poor prognosis and predict tumor progression. However, few studies based on systematic data have been conducted to explore m 6 A-related genes in melanoma and their general mechanisms. In this study, we summarized the expression of 10 Writer, Reader, and Eraser genes of m 6 A for the first time using Oncomine with more than 500 melanoma samples. We found that YTHDF1, ELF3, HNRN-PA2B1, YTHDF2, FTO, and HNRNPC were significantly altered (Fig. 2, Table 2). By integrating the results from Oncomine and GEO, the up-regulation of HNRNPA2B1 and YTHDF1 was identified (Fig. 3, Table 3). We further examined whether the combination of YTHDF1 and HNRNPA2B1 had higher efficacy in diagnosing melanoma by using linear regression and found the efficacy was improved from about 0.760 to 0.857 (P < 0.0001). This combination provides the first optimized method to distinguish patients with and without melanoma based on m 6 A regulatory genes rather than common molecular biomarkers such as serum lactate dehydrogenase.
We also analyzed the GO terms for all genes. For the terms associated with RNA processing, 518 genes were enriched in RNA metabolic processes, with positive regulation of RNA metabolism as the secondary term. The genes for both processes were up-regulated, indicating the occurrence of dynamic and active RNA metabolic processes in melanoma (Fig. 4a). YTHDF1, an m 6 A Reader, was found to augment EIF3C translation by binding to m 6 A-modified EIF3C mRNA and to  further promote tumorigenesis and metastasis in ovarian cancer [20]. Han et al. [21] demonstrated that targeting YTHDF1 with intestinal stem cells in established tumors could result in tumor shrinkage and longer survival. In our research, we found genes involved in the p53 signaling pathway such as CCNB1, CDK1, CHEK1, RRM2, and CDK2 were positively correlated with either YTHDF1 or HNRNPA2B1 (Fig. 4f ). This suggested that YTHDF1 or HNRNPA2B1 may interact with the related genes above and further influence the p53 signaling pathway, resulting in the development of melanoma (Fig. 4g). CCNB1 silencing activates the p53 signaling pathway, further inhibiting the proliferation of cells and promoting cell senescence in pancreatic cancer [22]. CDK1/2 inhibitors activate p53 in tumors, especially the inhibitors of CDK2, which maintains the balance of S-phase regulatory proteins and thus coordinates subsequent p53-independent G2/M checkpoint activation [23,24]. Hala et al. confirmed the existence of CHEK1/p53 association in human colorectal cancer in vivo and demonstrated that tumors lacking p53 had higher levels of CHEK1 [25]. The RRM1/RRM2B enzyme is capable of retaining activity in hypoxia and the replication pressure it induces is one of the factors proposed to promote selection stress that results in the loss of key ingredients for the DNA damage response, including p53 [26][27][28][29]. Therefore, YTHDF or HNRNPA2B1 might up-regulate CCNB1, CDK1, CHEK1, RRM2, or CDK2 (Fig. 4f ), all of which could inhibit the role of p53 in suppressing tumors and promote the development of melanoma. GO terms were found related to the activation of the cell cycle, DNA integrity checkpoints, histone phosphorylation, DNA damage response, signal transduction by p53 class mediation resulting in cell cycle arrest, regulation of transcription involved in the G1/S transition in the mitotic cell cycle, etc. (Fig. 4g). These findings suggested that the gene-interaction chain pathway might affect the proliferation as well as the normal physiological and metabolic processes of melanocytes, eventually resulting in the development of melanoma. All the genes in the identified pathways were also positively related to YTHDF1 or HNRNPA2B1 (Fig. 4f ), suggesting the expression of both genes might influence m 6 A modifications in mRNA genes and potentially result in melanoma heterogeneity at the cellular level. As shown in Fig. 4f-g, STAT1 and CD86 were negatively regulated by YTHDF1. Previous studies demonstrated that the loss of STAT1 activation or expression was found in malignant cells derived from various tumors, while the loss of CD86 expression resulted in a decrease in tumor-infiltrating T lymphocytes in diffuse B-cell large-cell lymphoma [30][31][32][33][34][35]. YTHDF1 might negatively regulate STAT1 and further facilitated melanoma development, but the roles of STAT1 and CD86 in the general network were not large enough to mitigate the trend of disease development. These genes together promoted the formation of melanoma, as shown in Fig. 4g. To fully explore the role of YTHDF1 or HNRNPA2B1, we assessed their mutations using the cBioPortal Tool (Fig. 5) and found the two genes were mainly amplified in melanoma. Major mutations occurred in the non-coding domain and affected the post-transcription of related proteins by influencing modifications to the proteins such as phosphorylation (Fig. 5c). In addition, the alterations of m 6 A-related genes influenced the patients' clinical features such as the tumor size, tumor stage, or treatment response (Table 3), indicating their important roles in the progression of melanoma.
However, this research also had some limitations. Firstly, the heterogeneity among the five studies included in Oncomine portal could not be ignored. The heterogeneity might result from patients' races, treatment or other factors that could have an influence on the expression levels of patients' transcript files. Secondly, there were few studies including gene qualification documents from both the normal and melanoma sites, therefore the number of samples used to examine the results got from Oncomine was limited. Thirdly,   we compared the levels of ten m 6 A regulatory genes between melanoma and the normal groups. Those ten genes included Writers as METTL14, Erasers as FTO, and Readers as YTHDF1, HNRNPA2B1, which were common genes that regulated m 6 A modification in RNA, while with the development new regulatory genes might emerge. Fourthly, the further experiments that explored the mechanism insight into the way m 6 A regulatory genes affected the development of melanoma needed to be performed in the future.

Conclusion
We systematically examined m 6 A-related genes in melanoma for the first time based on large datasets on publicly available databases. Several sources were used to confirm our results. Statistical analysis revealed the significance of YTHDF1 and HNRNPA2B1 and the combination of the two genes may provide a better approach to diagnose melanoma. Moreover, YTHDF1 or HNRN-PA2B1 might interact with RNA processing-related genes such as CDK1/CDK2 and further promote the formation and progression of melanoma. This study may provide new targets for the treatment of melanoma and the evidence to explain the pathophysiology of the disease.