Abstract This study aimed to analyze epigenetically and genetically altered genes in melanoma to get a better understanding of the molecular circuitry of melanoma and identify potential gene targets for the treatment of melanoma. The microarray data of [39]GSE31879, including mRNA expression profiles (seven melanoma and four melanocyte samples) and DNA methylation profiles (seven melanoma and five melanocyte samples), were downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) and differentially methylated positions (DMPs) were screened using the linear models for microarray data (limma) package in melanoma compared with melanocyte samples. Gene ontology (GO) and pathway enrichment analysis of the DEGs were carried out using the Database for Annotation, Visualization, and Integrated Discovery. Moreover, differentially methylated genes (DMGs) were identified, and a transcriptional regulatory network was constructed using the University of California Santa Cruz genome browser database. A total of 1,215 DEGs (199 upregulated and 1,016 downregulated) and 14,094 DMPs (10,450 upregulated and 3,644 downregulated) were identified in melanoma compared with melanocyte samples. Additionally, the upregulated and downregulated DEGs were significantly associated with different GO terms and pathways, such as pigment cell differentiation, biosynthesis, and metabolism. Furthermore, the transcriptional regulatory network showed that DMGs such as Aristaless-related homeobox (ARX), damage-specific DNA binding protein 2 (DDB2), and myelin basic protein (MBP) had higher node degrees. Our results showed that several methylated genes (ARX, DDB2, and MBP) may be involved in melanoma progression. Keywords: melanoma, DNA methylation, differentially expressed genes, gene ontology, pathway enrichment analysis, transcriptional regulatory network Introduction Melanoma, which is one of the most notoriously aggressive forms of skin cancer, arises from melanocytes (pigment-producing cells in the skin).[40]^1 Melanoma in young and middle-aged adults is characterized by increasing incidence and mortality rates worldwide.[41]^2 Early detection and excision are the most important factors in successful treatment of malignant melanoma.[42]^2 Therefore, it is essential to develop new and effective management strategies. Malignant melanoma results from the interplay of host, environmental, and genetic factors.[43]^3 Since the development of high-throughput genomic analysis, there has been significant progress in understanding of the genetic and molecular characteristics of melanomas, with the identification of signaling genes and pathways critical to their initiation and progression.[44]^1 For instance, mutations of v-raf murine sarcoma viral oncogene homolog B (BRAF) and neuroblastoma RAS viral (v-ras) oncogene homolog (NRAS) genes have been discovered with high frequency in nevi and cutaneous melanomas.[45]^4 Additionally, DNA methylation often occurs at the C5 position of cytosine in CpG dinucleotide. Methylation of the promoter regions can inhibit gene expression either by directly blocking the binding of transcriptional activators or silencing gene expression.[46]^5 Aberrant DNA hypermethylation is a frequent event found in melanoma and represents an essential mechanism utilized by shutting off different tumor suppressor genes (TSGs).[47]^6 Recent data have reviewed that aberrant promoter hypermethylation at cyclin-dependent kinase inhibitor 2A (CDKN2A) has been found in metastatic cutaneous melanoma samples.[48]^7 However, genetic understanding of melanoma remains incomplete. In the current study, we downloaded and reanalyzed the microarray data of [49]GSE31879 from a public database. The differentially expressed genes (DEGs), differentially methylated positions (DMPs), and differentially methylated genes (DMGs) between melanoma samples and human melanocyte samples were identified. Gene ontology (GO) and pathway enrichment analyses of the DEGs were carried out. Furthermore, the transcriptional regulatory network was constructed based on the interactions between the transcription factors (TFs) and the overlapping genes of DEGs and DMGs. We sought to get an increased understanding of epigenetic mechanisms in melanoma development and to identify potentially critical gene targets for melanoma. Further understanding of the underlying molecular mechanisms in melanoma may represent a novel and valuable therapeutic approach for melanoma. Materials and methods Microarray data Microarray data of [50]GSE31879, containing gene expression profiles along with methylation profiles, were downloaded from the Gene Expression Omnibus (GEO) database in the National Center of Biotechnology Information (NCBI), which was deposited by Tommasi et al on 6 September 2011.[51]^8 The two platforms used for gene expression profiling and methylation profiling are [52]GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array and [53]GPL13118 NimbleGen Human DNA Methylation 3×720 K CpG Island Plus RefSeq Promoter Array (100718_HG18_CpG_Refseq_Prom_MeDIP). A total of 23 samples as shown in [54]Table 1 were selected for the follow-up analysis, including seven primary human melanoma samples as well as five normal human melanocyte samples of DNA methylation profiles and seven primary human melanoma samples as well as four normal human melanocyte samples of mRNA expression profiles. Table 1. Sample data information Chip type Melanoma specimens Melanocyte samples mRNA expression profiles [55]GSM790543 [56]GSM790553 [57]GSM790544 [58]GSM790554 [59]GSM790545 [60]GSM790555 [61]GSM790548 [62]GSM790556 [63]GSM790549 [64]GSM790550 [65]GSM790552 DNA methylation profiles [66]GSM790528 [67]GSM790538 [68]GSM790529 [69]GSM790539 [70]GSM790530 [71]GSM790540 [72]GSM790533 [73]GSM790541 [74]GSM790534 [75]GSM790542 [76]GSM790535 [77]GSM790537 [78]Open in a new tab Data preprocessing and screening of DEGs and DMPs The original data that had already been log[2] transformed and loess normalized[79]^9 were downloaded. For the gene expression profiles, the probe sets were annotated according to the annotation information. The expression values of the probe sets were averaged if multiple probes were mapped to one gene. The gene expression matrix was obtained. For the DNA methylation profiles, the corresponding genes were identified according to the promoter information on the methylated probes. Linear models for microarray data (limma) software package provides an integrated solution for assessing differential expression for microarray studies.[80]^10 An unpaired t-test[81]^11 in the limma package[82]^10 was used to identify DEGs and DMPs between melanoma samples and normal melanocyte samples. DEGs and DMPs were determined using a threshold of |log[2] fold change (FC)| >1 and P-value <0.05. Functional and pathway enrichment analysis of DEGs The GO project is a community-based bioinformatics resource that uses a set of structured, controlled vocabularies and classifications to represent biological knowledge of genes, gene products, and sequences.[83]^12 Kyoto Encyclopedia of Genes and Genomes (KEGG) is an integrated database resource for integration and biological interpretation of large-scale molecular data sets.[84]^13 In the present study, in order to analyze the DEGs at the functional level, GO[85]^12 and KEGG pathway[86]^13 enrichment analyses were carried out on all the DEGs using the Database for Annotation, Visualization, and Integrated Discovery (DAVID).[87]^14 Overrepresented GO terms[88]^12 in biological process (BP), molecular function (MF), and cellular component (CC) categories as well as pathways were identified based on a hypergeometric distribution algorithm.[89]^15 The P-value <0.05 was chosen as the threshold. Identification of TF–mRNA interactive pairs The genes where the DMPs were located were named DMGs. According to the identified DEGs, DMPs, and DMGs, the overlapping genes among DMGs and DEGs were identified. University of California Santa Cruz (UCSC) genome browser ([90]http://genome.ucsc.edu) database provides convenient and timely access to high-quality genome sequence and annotations.[91]^16 In this study, combined with the TF binding site information at the UCSC genome browser database, the TF binding sites for the overlapping genes were detected and a TF–mRNA transcriptional regulatory network was constructed. Results Screening of DEGs and DMPs After analysis of the mRNA expression profiles and DNA methylation profiles, a total of 1,215 DEGs (including 199 upregulated and 1,016 downregulated genes; [92]Table S1) were identified in melanoma samples compared with melanocyte samples. Additionally, 14,094 DMPs (10,450 upregulated and 3,644 downregulated; [93]Table S2) were identified, and these DMPs corresponded to 3,004 DMGs. In addition, 167 overlapping genes among DEGs and DMGs were found of which 142 genes (85.03%) were downregulated ([94]Figure 1). The 167 genes had 493 DMPs of which 345 (69.98%) were upregulated. Figure 1. Figure 1 [95]Open in a new tab Venn diagram showing the number of DEGs (including up- and downregulated genes) and DMGs and the number of overlapping genes between DEGs and DMGs. Abbreviations: DEG, differentially expressed gene; DMG, differentially methylated gene. GO and pathway enrichment analysis of the DEGs The results of pathway enrichment analysis showed that the DEGs were enriched in different pathways. Top 10 (ranked by P-value) enriched pathways are presented in [96]Table 2. The upregulated genes were mainly enriched in sulfur metabolism, complement and coagulation cascades, and homologous recombination. Additionally, the downregulated genes were mainly enriched in the pathways associated with biosynthesis and metabolism, such as steroid biosynthesis, metabolic pathways, and sphingolipid metabolism. Table 2. Enriched KEGG pathway of DEGs Expression changes KEGG ID KEGG pathway Gene counts P-value Downregulation 4,142 Lysosome 29 1.241E–10 1,100 Metabolic pathways 112 3.651E–08 100 Steroid biosynthesis 8 8.371E–06 280 Valine, leucine, and isoleucine degradation 11 5.37E–05 600 Sphingolipid metabolism 9 0.0006024 62 Fatty acid elongation in mitochondria 4 0.0008343 563 GPI-anchor biosynthesis 6 0.0035269 4,962 Vasopressin-regulated water reabsorption 8 0.0050073 4,966 Collecting duct acid secretion 6 0.0053069 5,110 Vibrio cholerae infection 9 0.0054328 Upregulation 920 Sulfur metabolism 2 0.0035376 4,610 Complement and coagulation cascades 3 0.0120513 3,440 Homologous recombination 2 0.0160474 4,672 Intestinal immune network for IgA production 2 0.0438785 [97]Open in a new tab Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene; GPI, glycosylphosphatidylinositol. On the other hand, the DEGs were assigned to different GO terms. Top 10 GO terms are shown in [98]Table 3. GO analysis showed that the upregulated DEGs were mainly associated with the regulation of compound metabolic processes, such as regulation of nitrogen compound metabolic process and nucleobase-containing compound metabolic process. Besides, the downregulated DEGs were mainly linked to the metabolic process such as small molecule metabolic process, organonitrogen compound metabolic process and the oxidation-reduction process. Moreover, the downregulated DEGs were also found to be associated with GO terms such as pigment cell differentiation and developmental pigmentation. Table 3. Top 10 enriched GO terms of DEGs Expression changes GO ID Term Gene counts P-value Downregulation GO:0044281 Small molecule metabolic process 265 5.551E–16 GO:0044710 Single-organism metabolic process 382 1.31E–14 GO:0044711 Single-organism biosynthetic process 137 1.327E–13 GO:0006629 Lipid metabolic process 122 2.809E–13 GO:0019752 Carboxylic acid metabolic process 96 1.674E–10 GO:1901564 Organonitrogen compound metabolic process 193 2.056E–10 GO:0043436 Oxoacid metabolic process 102 7.596E–10 GO:0044255 Cellular lipid metabolic process 89 8.928E–10 GO:0055114 Oxidation–reduction process 97 1.11E–09 GO:0006082 Organic acid metabolic process 102 1.787E–09 Upregulation GO:0051171 Regulation of nitrogen compound metabolic process 53 1.608E–05 GO:0031323 Regulation of cellular metabolic process 62 2.143E–05 GO:0080090 Regulation of primary metabolic process 62 2.595E–05 GO:0006139 Nucleobase-containing compound metabolic process 67 2.932E–05 GO:0019219 Regulation of nucleobase-containing compound metabolic process 51 4.263E–05 GO:0031326 Regulation of cellular biosynthetic process 48 4.618E–05 GO:1901360 Organic cyclic compound metabolic process 69 6.055E–05 GO:0034641 Cellular nitrogen compound metabolic process 69 6.161E–05 GO:0009889 Regulation of biosynthetic process 48 6.311E–05 GO:0046483 Heterocycle metabolic process 67 7.705E–05 [99]Open in a new tab Abbreviations: GO, gene ontology; DEG, differentially expressed gene. Transcriptional regulatory network analysis According to the information of TF binding sites in the UCSC database, a total of 29 DMPs were identified to be covered by TF binding sites. The interactions between the TFs and the DMGs that corresponded to the 29 DMPs were represented in the transcriptional regulatory network ([100]Figure 2). The network was composed of 54 interactive pairs, 19 gene nodes, and 43 TF nodes. The nodes with higher degrees are shown in [101]Table 4. The results showed that genes such as Aristaless-related homeobox (ARX; degree =7), damage-specific DNA binding protein 2 (DDB2; degree =6), potassium voltage-gated channel, Isk-related family, member 4 (KCNE4; degree =5), and myelin basic protein (MBP; degree =5) could be regulated by more TFs. Moreover, of all the TFs in the regulatory network, CCAAT/enhancer binding protein (C/EBP), CCAAT/enhancer binding protein alpha (C/EBPalpha; degree =3), and glucocorticoid receptor-alpha (GR-alpha; degree =3) could regulate more genes. In particular, DDB2 was mainly regulated by the members of the E2F family of TFs. Figure 2. [102]Figure 2 [103]Open in a new tab Transcriptional regulatory network. Notes: The quadrate nodes indicate TF. The rounded nodes represent DMGs. Abbreviations: TF, transcription factor; DMG, differentially methylated gene. Table 4. Nodes with a higher degree in the regulatory network Gene Degree TF Degree ARX 7 C/EBPalpha 3 DDB2 6 GR-alpha 3 KCNE4 5 ATF-2 2 MBP 5 C/EBPbeta 2 HR 4 LMO2 2 SLC7A5 4 POU2F1 2 TBC1D16 4 PPAR-gamma1 2 BDNF 3 PPAR-gamma2 2 PCDH17 3 ZIC2 2 WNT5B 3 SCARB1 2 [104]Open in a new tab Abbreviation: TF, transcription factor. Discussion Melanoma arises from the accumulation of different gene alterations, and it is crucial to characterize the genetic changes during the progression of melanoma.[105]^17 In this study, integrated analysis was performed on the mRNA expression profiles and methylation profiles in human melanoma samples. Using the microarray platforms, we identified 1,215 DEGs (199 up- and 1,016 downregulated) and 14,094 DMPs (10,450 up- and 3,644 downregulated) in melanoma samples compared with melanocyte samples. The upregulated and downregulated DEGs were significantly associated with different GO terms and pathways, such as pigment cell differentiation, biosynthesis, and metabolism. Moreover, the transcriptional regulatory network showed that DMGs such as ARX, DDB2, and MBP could be regulated by more TFs, such as GR-alpha and E2F family. According to the literature, DNA methylation is an epigenetic process that can heritably result in gene expression changes without altering the DNA sequence.[106]^18 Moreover, DNA hypermethylation has been demonstrated to lead to abnormal silencing of several TSGs in most types of carcinoma.[107]^19 Consistent with this notion, our study showed that 14,094 DMPs (10,450 upregulated and 3,644 downregulated) were identified, and these DMPs corresponded to 3,004 DMGs. Thus, we concluded that most methylated positions in melanoma had increased methylation levels compared with melanocyte samples, suggesting that aberrant DNA methylation may be an essential mechanism of melanoma development and progression. ARX is a homeobox-containing gene expressed during brain development.[108]^20 The findings of Friocourt et al[109]^21 demonstrated that overexpression of ARX increased the length of cell cycle and ARX played an important role in controlling cortical interneuron migration and differentiation. A study had shown that melanoma cells that colonized the brain, which was a common target of metastases for melanoma patients, harbored numerous epigenetically and genetically altered genes.[110]^22 In addition, Marzese et al[111]^23 revealed that a genome-wide demethylation in low CpG density and also an increased methylation level in CpG islands were observed in melanoma brain metastasis (MBM) patients. Moreover, the results of our study demonstrated that ARX was downregulated and ARX was a DMG found in the regulatory network analysis. The network showed that ARX were regulated by the TF GR-alpha. Evidence demonstrates that GR-alpha is ubiquitously expressed and can alter expression patterns of target genes.[112]^24 Studies have indicated large quantities of glucocorticoid receptors in malignant melanoma cells.[113]^25^,[114]^26 Therefore, we suggested that methylation of ARX may play an essential role in melanoma progression and in melanoma metastasis to the brain. In addition, DDB2 encodes a protein that is necessary for the repair of ultraviolet (UV)-damaged DNA or nucleotide excision repair.[115]^27 Luijsterburg et al[116]^28 had demonstrated that DDB2 promoted chromatin decondensation at UV-induced DNA damage. Additionally, Itoh et al[117]^29 revealed that DDB2 was a haploinsufficient tumor suppressor and could control spontaneous germ cell apoptosis and tumor incidence. Besides, studies have demonstrated that UV can cause DNA damage via photosensitized reactions that lead to the production of oxygen radical species and UV exposure has been implicated as a potential role in the pathogenesis of melanoma.[118]^30^,[119]^31 In keeping with the previous studies, our study showed that DDB2 was downregulated in human melanoma samples and DDB2 was a DMG found in the regulatory network analysis. Furthermore, DDB2 could be regulated by members of the E2F family of TFs. An active E2F element is revealed in the DDB2 promoter.[120]^32 Moreover, E2Fs regulate expression of genes involved in DNA repair.[121]^33 Taken together, we suggested that DDB2 methylation may play an important role in melanoma progression via affecting repair of DNA damage and cell fate. Furthermore, our study showed that MBP was another DMG in melanoma samples. MBP, the second most abundant protein in myelin of the central nervous system, is in charge of adhesion of cytosolic surfaces of the multilayered compact myelin. The study of Pearson et al[122]^34 had showed that MBP could be a substrate of mitogen-activated protein kinase (MAPK). Besides, Fecher et al[123]^35 had revealed that MAPK pathway was one of the several potentially targetable pathways in melanoma. As shown in the transcriptional regulatory network in this study, MBP could be regulated by TF GR-alpha. The MAPK pathway mediates behavioral effects of glucocorticoids.[124]^36 In line with the previous studies, we suggested that MBP may be essential in melanoma progression via involvement in the MAPK pathway. In recent years, DNA methylation alterations as an epigenetic abnormality have come into focus and it is indicated that methylation of a subset of genes may be useful biomarkers as tools for screening and early detection of cancer.[125]^37 Moreover, such DNA methylation-based biomarkers of disease can be also used for risk stratification and may have utility in personalized medicine.[126]^38 In addition, both genetic and epigenetic abnormalities can be used in monitoring the effect of chemopreventive drugs.[127]^38 Therefore, the findings of this study may represent a novel and valuable therapeutic approach for melanoma. Conclusion Our study identified several methylated genes (ARX, DDB2, and MBP) that may be involved in the mechanism of melanoma progression. Aberrant DNA methylation was an essential mechanism of melanoma development and progression. Methylation of ARX may play an essential role in melanoma progression and in melanoma metastasis to the brain. Additionally, DDB2 methylation may play an important role in melanoma progression via affecting repair of DNA damage and cell fate. Moreover, MBP may be essential in melanoma progression via involvement in the MAPK pathway. However, further experiments and investigations are needed to verify our results, exploring useful novel DNA methylation predictors and effective treatment strategies. Footnotes Disclosure The authors report no conflicts of interest in this work. References