Identification of epigenetically altered genes and potential gene targets in melanoma using bioinformatic methods.

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 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.

Introduction

Melanoma, which is one of the most notoriously aggressive forms of skin cancer, arises from melanocytes (pigment-producing cells in the skin).1 Melanoma in young and middle-aged adults is characterized by increasing incidence and mortality rates worldwide.2 Early detection and excision are the most important factors in successful treatment of malignant melanoma.2 Therefore, it is essential to develop new and effective management strategies.

Malignant melanoma results from the interplay of host, environmental, and genetic factors.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.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.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.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).6 Recent data have reviewed that aberrant promoter hypermethylation at cyclin-dependent kinase inhibitor 2A (CDKN2A) has been found in metastatic cutaneous melanoma samples.7 However, genetic understanding of melanoma remains incomplete.

In the current study, we downloaded and reanalyzed the microarray data of 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 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.8 The two platforms used for gene expression profiling and methylation profiling are GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array and 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 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	GSM790543	GSM790553
	GSM790544	GSM790554
	GSM790545	GSM790555
	GSM790548	GSM790556
	GSM790549
	GSM790550
	GSM790552
DNA methylation profiles	GSM790528	GSM790538
	GSM790529	GSM790539
	GSM790530	GSM790540
	GSM790533	GSM790541
	GSM790534	GSM790542
	GSM790535
	GSM790537

Data preprocessing and screening of DEGs and DMPs

The original data that had already been log2 transformed and loess normalized9 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.10 An unpaired t-test11 in the limma package10 was used to identify DEGs and DMPs between melanoma samples and normal melanocyte samples. DEGs and DMPs were determined using a threshold of |log2 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.12 Kyoto Encyclopedia of Genes and Genomes (KEGG) is an integrated database resource for integration and biological interpretation of large-scale molecular data sets.13 In the present study, in order to analyze the DEGs at the functional level, GO12 and KEGG pathway13 enrichment analyses were carried out on all the DEGs using the Database for Annotation, Visualization, and Integrated Discovery (DAVID).14 Overrepresented GO terms12 in biological process (BP), molecular function (MF), and cellular component (CC) categories as well as pathways were identified based on a hypergeometric distribution algorithm.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 (http://genome.ucsc.edu) database provides convenient and timely access to high-quality genome sequence and annotations.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; Table S1) were identified in melanoma samples compared with melanocyte samples. Additionally, 14,094 DMPs (10,450 upregulated and 3,644 downregulated; 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 (Figure 1). The 167 genes had 493 DMPs of which 345 (69.98%) were upregulated.

Figure 1

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 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

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 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

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 (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 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

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

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.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.18 Moreover, DNA hypermethylation has been demonstrated to lead to abnormal silencing of several TSGs in most types of carcinoma.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.20 The findings of Friocourt et al21 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.22 In addition, Marzese et al23 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.24 Studies have indicated large quantities of glucocorticoid receptors in malignant melanoma cells.25,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.27 Luijsterburg et al28 had demonstrated that DDB2 promoted chromatin decondensation at UV-induced DNA damage. Additionally, Itoh et al29 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.30,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.32 Moreover, E2Fs regulate expression of genes involved in DNA repair.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 al34 had showed that MBP could be a substrate of mitogen-activated protein kinase (MAPK). Besides, Fecher et al35 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.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.37 Moreover, such DNA methylation-based biomarkers of disease can be also used for risk stratification and may have utility in personalized medicine.38 In addition, both genetic and epigenetic abnormalities can be used in monitoring the effect of chemopreventive drugs.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.