Integrative Genomics Identifies Gene Signature Associated with Melanoma Ulceration

Background: Despite the extensive research approaches applied to characterise malignant melanoma, no specific molecular markers are available that are clearly related to the progression of this disease. In this study, our aims were to define a gene expression signature associated with the clinical outcome of melanoma patients and to provide an integrative interpretation of the gene expression-, copy number alterations-, and promoter methylation patterns that contribute to clinically relevant molecular functional alterations.


Introduction
Melanoma is an aggressive, therapy-resistant malignancy of the melanocytes. The incidence of this disease has been increasing worldwide, resulting in a growing public health problem [1,2]. Interactions between molecular and environmental risk factors that promote melanomagenesis are the subject of ongoing research [3]. There are currently no systemic therapies available to significantly extend the life expectancy of patients with advanced melanoma. Therefore, early diagnosis and conventional treatments remain the key to improved survival for all affected individuals [1,4,5]. Recently, several high-throughput genomic and gene expression studies have been designed to reveal the molecular mechanisms involved in melanoma progression [6][7][8][9][10][11][12][13][14].
DNA microarray technology, a high-throughput screening methodology, offers an opportunity to simultaneously analyse the gene expression levels of thousands of genes in the cancer genome. This technology has been extensively used in cancer research to identify tumour subclasses, predict disease outcomes and identify genes associated with drug resistance [15,16]. Previous transcriptome analyses have provided valuable information for the assessment of patient group classifications, such as subgroups of patients who are likely to respond to a particular therapy [17]. Many research groups have used different gene expression platforms to obtain a better understanding of the stepwise tumourigenesis involved in melanoma development. These studies have identified cohorts of genes that facilitate distinction of benign nevi from malignant melanomas [18,19], sub-classification of metastatic melanomas into distinct subgroups [6,[19][20][21][22][23] and prediction of distant metastasis-free survival [20,24]. Most of these research groups have found at least two naturally occurring subgroups using unsupervised clustering. Some researchers have been unable to find any association between subgroups with any clinical variables [19,21]. However, other groups found that the clusters differed in terms of how advanced tumour were, defined by thickness and survival [20], and the transition from the radial growth phase (RGP) to vertical growth phase (VGP) [6,22,23]. Microarray studies have been limited in utility because of the lack of concordance from one study to the next. This limitation suggests that in addition to tumour heterogeneity, inconsistencies in the experimental design, sample preparation and platforms used for analysis have contributed to the diversity of data.
It has recently become apparent that melanoma heterogeneity may be caused by both alterations in gene expression and distinct genomic changes [25][26][27][28][29]. Taking a global view of the available malignant genome and transcriptome has proven to be valuable for understanding the intricate process of neoplastic transformation and for obtaining profound downstream benefits from these discoveries. Despite rapid progress in exploring commonly amplified and deleted regions and epigenetic changes in melanoma, integrated analyses have focused only on melanoma cell lines. Simultaneous studies have been aimed at exploring commonly affected genes that exhibit copy number changes and altered gene expression [30][31][32], but there has been no study conducted on melanoma tissues [33].
In this study, our aim was to define the gene expression signature of previously untreated primary melanomas in patients with a known clinical history using the Affymetrix GeneChip Human Genome U133 Plus 2.0 expression arrays platform. We extended these analyses and correlated the transcriptomic data with copy number alterations defined using the NimbleGen HG18 CGH 4x72K WG Tiling v2.0 array and the CpG methylation of promoters analysed by the Illumina GoldenGate Methylation Assay.
Integrated analysis of different types of data that are characteristic of the same tumour might allow us to identify molecular pathways and genetic aberrations that could influence the biological and clinical behaviour of melanomas and pinpoint progression-related genes, facilitating molecular classification of the disease. These candidate genes could ultimately be utilised as therapeutic targets for treating human malignant melanoma.

Melanoma Samples
Thirty-six primary melanomas were involved in our gene expression studies and 17 out of the 36 primary melanomas was analysed by tiling array CGH and DNA methylation assay. Tumour tissues were obtained from the Department of Dermatology, University of Debrecen, Hungary. All human studies were conducted in accordance with principles outlined in the Declaration of Helsinki and were approved by the Regional and Institutional Ethics Committee of the University of Debrecen Medical and Health Science Center and was conducted according to regulations (Protocol #2836-2008. Written informed consent was obtained from each patient. The clinicopathological data on the 36 primary melanomas are summarised in Table 1 according to the new melanoma TNM staging system [34,35].
High quality total RNA and genomic DNA were isolated from fresh tissues. Detailed protocols for routine procedures and quality controls are available under the Text S1.

Microarray Based Experiments
Gene expression analysis. Gene expression analysis was carried out using Affymetrix GeneChip Human Genome U133 Plus 2.0 expression arrays The microarray data was published in the Array Express Archive repository (a MIAME compliant database) with accession number E-MTAB-946.
Analyses of gene expression data were performed using GeneSpring 7.3 (Agilent Technologies, Santa Clara, CA, USA) and BRB Array Tools 4.3.0 developed by Dr. Richard Simon and BRB Array Tools Development Team. The raw intensity values from each chip were normalised to the 50th percentile of the measurements to reduce chip-wide variations in intensity. Each gene was normalised to the average level of that particular gene to allow comparison of relative changes in gene expression levels between different conditions. Only genes showing detection flag present in at least 50% of the samples were used in further analyses (25,886 probe sets). The log 2 transformed data of all samples were classified according to their expression pattern via hierarchical clustering using average linkage and Pearson's correlation. The genes, that expressed differently among the predefined melanoma groups (Table 1) at the 2-fold change and p#0.001 level after Benjamini and Hochberg multiple testing correction, were assessed by univariate t-test for each gene and were visualized by Volcano plots. Principal component analysis (PCA) of differentially expressed genes was applied to verify the proportion and the significance of the first 3 components covering the total variation.
Kaplan-Meier survival curves and Cox proportional univariate test were performed to whether the gene expression of a particular gene significantly influences survival at the p,0.001 level. To control for patients' age, gender and clinical covariates on survival and to predict a survival prediction we used the Penalized Cox Regression model.

Copy Number Analysis by Array CGH
An array CGH (NimbleGen HG18 CGH 4x72K WG Tiling v2.0) and the calculation of ratio values were performed on genomic DNA samples by Roche NimbleGen (Reykjavik, Iceland). All of the aCGH data was published in the Array Express Archive repository (a MIAME compliant database) with accession number E-MTAB-947.
Statistical analyses were conducted with Nexus Copy Number 5.1. For segmentation, we applied the FASST2 algorithm. To adjust the sensitivity of the segmentation algorithm, we determined a significance threshold of 1.0E-6 and specified 1000 kb as the maximum spacing between adjacent probes. The minimum number of probes per segment was 5 to eliminate small CNAs. To detect gains and losses, we set the thresholds as 60.3 for gains and losses.
Significantly different copy number events between ulcerated versus non-ulcerated melanomas were identified by applying Fisher's exact test, and FDR adjustment for correcting multiple testing was calculated using the Nexus Copy Number 5.1 Comparison feature.

Investigation of the Promoter Methylation Patterns
The quantitative methylation status of the 1505 CpG sites corresponding to 807 cancer-related gene promoters corresponding to 14 primary melanomas was performed with the Illumina GoldenGate DNA Methylation arrays as described previously [36]. Briefly, genomic DNA from tumour tissues was prepared by overnight proteinase K treatment, phenol-chloroform extraction, and ethanol precipitation. Sodium bisulfite modification was performed on 500 ng DNA using the EZ DNA Methylation-Gold Kit (Zymo Research). BeadStudio v3.2 software (Illumina) was used for initial filtering and clustering analysis as described previously [36]. We identified 45 overlapping genes represented by 98 different CpGs between our downregulated gene list (987 genes) and the Illumina Assay. To estimate if the 98 overlapping CpG sites are methylated differently, we performed univariate ttests for each gene by BRB Array Tools. Statistically significant difference was considered when the P value was less than 0.05 below the 0.1 Benjamini and Hochberg FDR.

Integration of Gene Expression, Copy Number Variation and Promoter Methylation Datasets
To study the relationship between DNA copy number gains/ losses and mRNA levels, we exported the median of the replicate probe log 2 ratios and the expression values corresponding to the same genomic region for determining Pearson's correlation. The genelist generated through this analysis was imported to the Nexus 5.1 package. We assessed the genomic locations that exhibited significant copy number losses and concentrations of downregulated genes simultaneously. The statistics-generated P value was based on the number of deregulated transcripts located in significantly deleted genomic regions after multiple testing correction.
To verify the effects of cis-regulatory copy number elements on gene expression and to assess trans-acting copy number alterations, we used the 'lol' R package applying cross-validation as optimizer [37][38][39]. Detailed description of this approach is available under the Text S2.
To assess relation between gene expression and promoter methylation patterns, we applied Pearson's correlation for each of the 98 CpG sites corresponding to the 45 overlapping genes between the datasets.

Pathway Analysis of the Gene Expression Data
Functional characterisation of differentially expressed genes (p,0.05) that showed differences of greater than 2-fold between the melanoma groups (ulcerated melanomas versus melanomas without an ulcerated surface) was performed using Ingenuity Pathways Analysis software (IPA, Ingenuity Systems: www. ingenuity.com) and the Database for Annotation Visualisation and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov) Gene Functional Classification Tool. Datasets containing gene identifiers (Affymetrix ID) and the corresponding expression values were evaluated. Each gene identifier was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base (IPKB). Functional analysis of the genes was based on the collection of data summarised and stored in the IPKB database. Genes identified using the IPA program were categorised according to location, cellular components and their reported biochemical, biological, and molecular functions.

Validation of the Gene Expression Microarray Results by QPCR
The expression status of selected genes was performed using quantitative real-time PCR in all melanoma samples with the ABI PrismH 7900HT Sequence Detection System (Applied Biosystems, Carlsbad, California, USA). Reverse transcription (RT) was carried out on total RNA (600 ng) using the High Capacity cDNA Archive Kit, according to the protocol of the supplier (Applied Biosystems, Carlsbad, California, USA). A TaqMan Low Density Array with pre-designed TaqManH Gene Expression Assays (Applied Biosystems, Carlsbad, California, USA) was used to perform qPCR for 17 genes (Table S1).
Selection of genes covered several genes from Table 2 Table S1 comprising 987 genes that had been proved to be downregulated in ulcerated melanomas by microarray technique.
The house keeping genes b-actin (Hs99999903_m1) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH; Hs99999905_m1) were used as controls for accurate normalisation of the gene expression data [40]. Pearson's correlation was used to evaluate the strengths of the relationships between the PCR and the microarray expression data. The Student's t-test was performed to validate statistically significant differences among the predefined sample groups, and the data were considered significant at p,0.05.

Gene Expression Alterations of Primary Melanomas
Affymetrix oligonucleotide microarrays were used to investigate the gene expression patterns in 36 primary melanomas. The follow-up period of the patients was 5 years. The clinicopathological parameters of the samples are summarised in Table 1.
After the normalisation and quality control steps were performed for the microarray data, 25,886 probe sets were used in further analyses. Univariate t-tests were performed on logtransformed expression values of the probes to determine differentially expressed genes (P value cut off 0.001, greater than 2-fold change), among all the predefined clinical groups ( Table 1). The analysis resulted in 1,080 genes that expressed differentially between ulcerated and non-ulcerated melanomas with FDR false discovery rate of less than 5%. Our results are visualized by Volcano plot (Figure 1 A). The majority of the differentially expressed genes were downregulated (987 genes, Table S1), whereas only 93 genes (Table S2) were significantly upregulated in ulcerated specimens compared to non-ulcerated melanomas. Alternatively, hierarchical clustering of the 1,080 differentially expressed genes distinguished also the primary melanoma subgroups, and melanomas with a tumour surface ulceration were clustered mainly together ( Figure 1B). To measure how ulcerated and non-ulcerated groups can be distinguished according to the 1080 differentially expressed genes, we performed Principal Component Analysis (PCA). Based on the global test of clustering (p,0.001), the first 3 components covered the 59% of the total variance.
Using univariate t-tests and PCA for further verification we were unable to describe differentially expressed genes among other sample groups categorized by histologycal subtype, Breslow thickness, Clark level and metastasis formation.
For survival analysis we established Penalized Cox Regression model that revealed no relationship between patients' survival and gene regulation when we considered patients' age, gender and other clinical cofounders, including survival.

Functional Categorisation of Differentially Expressed Genes
The functional annotation of upregulated genes was performed using the web-accessible programs of the Database for Annotation Visualisation and Integrated Discovery (DAVID). Of the 93 overexpressed genes, 79 could be identified by DAVID, and 23 of these genes belonged to the well-known KEGG pathways (Kyoto Encyclopaedia of Genes and Genomes). These data are summarised in Table 3. The overexpressed genes in the ulcerated melanoma cluster included many genes with reported functional roles in cell cycle regulation and proliferation, such as cell division cycle 25a (CDC25a) and MAPKK12. Out of the 93 genes, the Table 2. Networks and pathways of downregulated genes associated with ulcerated melanoma surface.
According to the functional role of the genes, 26 networks were identified via IPA analysis. These networks are related to hair and skin development and functions, dermatological diseases and conditions, cellular development, cellular growth and proliferation, cardiovascular system development and function, cell morphology, cellular assembly and organisation, and cell-to-cell signalling and interactions (Table 2). For example, the ''Cancer, Cellular Movement, Cellular Growth and Proliferation'' network encompasses EGFR, FGFR2 and other key genes, such as AURKA, a cell cycle checkpoint mediator, and SPP1, LAMA3 and GJA1, which play a role in cell communication.
In the analysed melanomas, the ''Cell adhesion, Cadherin'' network was significantly enriched for genes for which decreased expression is correlated with ulceration (p,0.0001). This network involves genes playing a significant role in the Wnt (WNT16, Figure 1. Gene expression patterns of primary melanomas. A: The Volcano plot displays the relationship between fold changes and statistical significance: the x-axis visualizes the log 2 values of the fold change while the y-axis belongs to the log 10 values of significance. The green line of the xand y-axes indicate greater than 2-fold change between ulcerated and non-ulcerated samples and P value cut off 0.001 revealed by univariate t-tests for each gene. The Benjamini and Hochberg false discovery rate was applied as multiple test correction. B: Unsupervised hierarchical clustering of the 1080 genes which expressed differentially among ulcerated and non ulcerated melanomas. Samples are displayed vertically and genes are displayed horizontally. Annotation track marks ulcerated samples with violet and non-ulcerated with green colour. The colour in each cell of the table represents the median-adjusted expression value of each gene. Yellow indicates increased expression relative to the median, while blue represents decreased expression. C: Principal Component Analysis for distinction of ulcerated (green dots) and non-ulcerated (blue dots) samples based on the 1080 differentially expressed genes. The analysis revealed that according to the first 3 components which covered the 59% of the total variance the 2 groups were significantly different (p,0.001). doi:10.1371/journal.pone.0054958.g001 WNT3, WNT4, WNT5a), beta-catenin, EGFR, FZD10 and p21 signalling pathways ( Figure S1). Interestingly, the members of the Wnt family that have been reported to be overexpressed in the majority of melanomas were overexpressed only in one subset of our melanoma samples. This analysis also revealed that the expression levels of distinct components of the ''Junctional Mechanism Regulation Pathway'' (i.e., ZO-2, OCLN, GJA1, CLDN1, DSP, JUP) were differentially expressed in the analysed tumours ( Figure S2) ZO-2 (TJP2) is a tight junction plaque adapter for the recruitment of cytosolic molecules implicated as integral membrane proteins involved in cell signalling, and Occludin (OCLN) and Claudin 1 (CLNDN1) serve as linkers of the actin cytoskeleton to the plasma membrane and contribute to cytoskeleton organisation, small GTPase-mediated signal transduction, and cellular component organisation. Downregulation of laminin 5, collagen IV and ezrin is related to integrin-mediated cell-matrix adhesion, playing a key role in cell surface structure adhesion, migration and organisation. Compared to melanomas without ulceration, ulcerated melanomas were characterised by downregulation of members of the kallikrein serine protease family (i.e., KLK7, KLK5, KLK11), which are related to angiogenesis and the degradation of extracellular matrix components. The top five identified networks that were affected in tumours with ulcerated surfaces were centred on key genes involved in the p53, Nf-kappaB, WNT/beta-catenin pathways.

Validation of Expression Microarray Results via Quantitative RT-PCR
The microarray results were validated by performing qRT-PCR measurements of the 17 downregulated genes. Each of the 36 tumour samples included in the microarray analysis was tested using a qPCR based TaqMan Low Density Array.
For all of the 17 genes, the expression levels differed between melanomas with or without ulceration (more than a 2-fold change, p,0.05). These data showed a strong correlation with the microarray results. Six of the genes displayed a highly significant (p,0.01) expression difference between the melanoma groups ( Figure 2). Two of the genes (KLK11 and KLK5) belong to the protein-degrading kallikrein family. The KLF4 transcription factor showed decreased expression (by PCR as well) and is capable of both activating and repressing genes involved in cell-cycle regulation and differentiation. Furthermore, fibroblast growth factor receptor 2 (FGFR2); ARHGEF4, which acts as a guanine nucleotide exchange factor; and SDC1, a member of the syndecan proteoglycan family, were also expressed at significantly lower levels in ulcerated tumours. These genes are part of the ''Cancer, Cellular movement, Cellular Growth and Proliferation'' network and the ''Hair and Skin Development and Function'' network.

Comparison of Gene Expression Changes to Copy Number Alterations and DNA Methylation Patterns in Primary Melanomas
We performed a detailed aCGH analysis on 17 melanomas and collected integrated array CGH and gene expression data from the same tumour to characterise genetic alterations that might be associated with gene expression changes.
To determine the influence of copy number aberrations on the gene expression changes in primary melanomas, we sought to identify the genes whose expression was significantly correlated with copy number changes. By correlating the intensity of the aCGH ratios (median log2) with the expression levels of previously identified downregulated genes (987 genes), we identified 150 genes whose expression was significantly and positively correlated with copy number changes in ulcerated melanomas (p,0.05, Table S3). The gene list generated through this analysis was significantly enriched in genes mapping to chromosomes 6q and 10q (Table 5). Figure 3 shows the frequency distribution of copy number losses in ulcerated (red line) and non-ulcerated (green line) melanomas. We found more narrow regions that were deleted on the long arm of chromosome 6 (6q14.1-q14.2, 6q16.3-q21, 6q22.31-q22.32, 6q23.3 and 6q24.2), consisting of 9 deregulated genes (ELOV4, ME1, TPBG, AIM1, TPD52L1, IL20RA, HEPB2, PERP and UTRN). A significant association between genomic losses and the downregulation of genes on chromosome 10 was only observed for one gene (ABLIM1, localised to 10q25).
Lasso regression, performed by lol R-package, identified both copy number gains (Figure 4 B) and losses (Figure 4 A) as top   and the copy number gain of 17q22 (Figure 4 B).Cis-affecting copy number losses given in Table 5 were verified by Lassoregression as we revealed Score cis values more than 0.5 for each region.
To define the influence of the epigenetic alterations on melanoma ulceration, we estimated if the 98 overlapping CpG sites were methylated differently by univariate t-tests for each genes. Figure 5A plots a clustered heat map showing the lack of differences between the sample groups with an exception of JAK2 and IL1RN genes ( Figure 5B), however, significant differences did not remained after adjustment of multiple comparison ( Figure 5C). Furthermore, we correlated the gene expression data with the observed epigenetic changes (Illumina Golden Gate Cancer Panel 1) in the same tumour. Of the 987 downregulated genes identified using the Affymetrix microarray, we found 45 common genes represented by 98 different CpGs on the methylation array platform. A strong or medium inverse correlation (a medium inverse correlation was assessed if r#20.30 and a strong inverse correlation if r#20.50 by Pearson's correlation) between gene expression and methylation levels was detected in the case of 11 genes corresponding to 17 different CpG sites. Detailed list of these genes are shown in Table 6.

Discussion
Malignant melanoma develops from the malignant transformation of melanocytes, the pigment-producing cells in the basal epidermal layer of human skin, and is one of the most aggressive human cancers. Melanoma arises from the accumulation of different alterations in genes that are critical for cell proliferation, differentiation, and cell death [41,42]. To develop treatments for advanced melanoma and to increase survival in metastatic melanoma patients, it is important to characterise the genetic and gene expression changes leading to each progressive step of the disease.
During the last decade, microarrays have become the technology of choice for the selection of the genes responsible for the behaviour of malignant lesions. Previously published microarray results have shown highly different and frequently contradictory ''melanoma signatures'' containing many genes whose functional roles in melanoma progression have not been well characterised. Furthermore, the gene sets found in different microarray studies overlap poorly with each other [20,21,43,44]. In this study, we defined the gene expression profiles in 36 primary melanoma tissues and correlated these data with DNA copy number alterations and epigenetic changes in the selected samples. Using the Affymetrix microarray platform, we identified 1,080 genes with expression value that was significantly associated with the tumour surface ulceration of the primary tumours. The American Joint Committee on Cancer has recently reported that after Breslow thickness, the presence or absence of ulceration is the second most powerful independent predictor of survival for melanoma patients [35,45]. While the molecular background of melanoma ulceration is still unclear, according to the most recent study involved 4,661 patients, not only the presence but the extent of ulceration could be independent predictors of survival [46]. Furthermore, ulceration was reported to be significant predictive factor for outcome of adjuvant interferon treatment [47]. All of these results prompt further efforts in designing additional studies to obtain better insights into the molecular level of ulceration.
Our functional analysis revealed that the majority of the differentially expressed genes were downregulated (987 genes) and were involved in the p53, Nf-kappaB, and WNT/beta-catenin pathways. In agreement with previous studies, we also observed downregulation of pro-apoptotic genes (e.g., TPL73L and P53AIP1), members of the forkhead box gene family (e.g., FOXQ1) and members of the TNFRSF25 and TNFSF10 protein families [7].
According to the literature, osteopontin expression is closely associated with shortened relapse-free survival and other histological variables, including ulceration, that predict relapse and are associated with metastatic dissemination [48][49][50]. In this study, we have demonstrated that osteopontin is significantly overexpressed in ulcerated melanomas. Our finding agrees with previous studies and supports the assumption that osteopontin expression is a potential biomarker of melanoma progression [48,50].
For obtaining a better insights into the molecular mechanisms that might be responsible for the aggressive phenotype associated with the observed gene expression signatures of ulcerated melanomas, we extended our study to search for DNA copy number alterations and epigenetic changes that might influence transcription. Array CGH is currently the best tool for searching for non-random DNA copy number alterations in cancer genomes and finding new genes that harbour copy number gains or losses with prognostic relevance [51][52][53]. As a result of improved resolution, there is currently an abundance of CGH data available in both large genomic regions and for uniquely affected genes. However, these studies focused little on elucidating the genomic alterations that contribute to gene regulation. Valseria et al. were the first to highlight SCNA (somatic copy number alteration) genes. The downregulation and upregulation of these genes is accompanied by genomic losses or gains in melanoma [9]. The present study revealed new genes that might give better insight into copy number alteration (CNA)-induced gene regulation.
Because the experiments were based on a metastatic cell line, the aforementioned authors failed to find any progression-related clinical effects associated with these SNCA genes.
The main purpose of this part of the study was to give functional relevance to the copy number events by providing a statistically powerful integrated interpretation of gene expression changes and copy number alterations. Based on our aCGH analyses, we found more regions on chromosome 6q and only one region on chromosome 10q that showed significantly different loss of copy numbers between the two clinical subgroups. It is important to note that these regions were enriched for the downregulated transcripts and significantly correlated with our previously defined gene expression results. There are a total of 36 genes in these two regions, among which we identified 10 downregulated transcripts using the Affymetrix microarray. It is important to note that these genes include TPBG, PERP and UTRN, which are involved in cell-cell and cell-matrix adhesion. PERP also functions as a p53induced apoptosis effector molecule. Furthermore, TPD52L1 participates in apoptosis followed by nuclear fragmentation. According to the literature, IL20RA is a tissue-specific Interleukin Receptor that is highly expressed in normal skin [54]. The remaining 5 genes could not be directly associated with carcinogenesis; however, a detailed, functional analysis is needed to precisely define the role of each of the new genes we have identified. It is worth noting that previous high-resolution aCGH experiments have also provided convincing evidence of the importance of copy number alterations at 6q, as this region often experiences hemizygous deletion, with MYB1 being the only gene specified on 6q23 [53,55]. These experiments also suggest that copy number alterations in MYB1 are an important discriminator between melanomas and nevi, as validated by FISH in 123 melanomas and 110 nevi. This study supports the assumption that 6q23 can assist in the diagnosis of melanoma. In agreement with previous findings, we further support the idea that the deletion of 6q23 is an important alteration in melanoma progression. However, it should also be pointed out that additional new genes (IL20RA, HEPB2 and PERP) that are located in this region and are downregulated might make significant contributions to melanoma progression. In addition to 6q23, we found a significant association between the tendency towards gene deregulation and DNA sequence deletions at 6q14, 6q16, 6q22, 6q24 and 6q25. Somatic copy number alterations at 10q were characterised in Figure 5. Relationship between DNA methylation and melanoma ulceration. Clustered heatmap, performed to demonstrate promoter methylation patterns, shows lack of differences between ulcerated and non-ulcerated sample groups. The heatmap is based on the univariate t-tests performed for each sites (specific for 45 independent genes) that overlapped with the gene expression results. Two genes (IL1RN and JAK2) demonstrated increased methylation for the ulcerated sample group, however, significant differences did not remain after adjustment for multiple comparison. doi:10.1371/journal.pone.0054958.g005 many cases, but particularly attributed to the PTEN tumour suppressor gene [56].
While systematic correlation analysis among the copy number events and the corresponding genes captures cis-effects it is also important to measure trans-acting copy number alterations (DNA aberrations on one gene affect the mRNA expressions of other genes); for the latter purpose we used ''lol'' R package developed by Yuan et al. [39], which provided also an appropriate method for validating the cis-acting elements.
Regarding the trans-regulatory copy number elements, the top scoring (Score trans ) alterations included copy number gains and losses, however, both types of alterations were associated with transcriptomic deregulation. The week immunostaning of EZR protein in the majority of melanomas was mentioned by a single study, however, this phenomenon was not associated with ulceration [57]. In our study, we showed that EZR (6q25.3) deregulation was related to the copy number loss of 6q27 and the copy number gain of 17q22 suggesting the trans regulatory effects of these copy number alterations. Frequent loss of 9p21 region is a well characterized copy number aberration in melanomas [58,59], however, it has not been mentioned as a trans-regulatory element for Anoctamin 1 gene (11q13.3) so far. As our results demonstrated, the transcriptomic regulation of this gene can be possibly disturbed also by copy number loss of 9q21.11 and 7q11.23 regions.
Parallel to identifying genetic alterations via aCGH, we sought to examine our external dataset of methylation patterns of 1,505 CpG sites corresponding to 807 selected, cancer-related genes. Despite the substantial number of candidate genes that have been shown to be methylated in melanoma, only a few genes were identified as contributors to disease progression [60]. The major reason for this might be that most of these studies have been based on melanoma cell lines and have provided ambiguous results. Few direct experiments have focused mainly on elucidating the differences in methylation patterns between melanomas and benign skin lesions. However, these studies have had a greater impact regarding understanding what drives melanocytic lesions to develop into malignant tumours rather than on the progression of tumours [61].
By correlating 45 overlapping genes between the gene expression dataset and the methylation assay, we identified 11 genes exhibited inverse relationships. Two (EPHB3, EPHB6) of the 11 genes are members of the erythropoietin-producing hepatocyte kinase B (EPHB) receptor family. Both of these genes play a suggested role in tumour suppression via regulating cell adhesion and migration. In addition, these genes have downstream effects on several members of the kallikrein family, which we also found to be downregulated in ulcerated melanomas in line with other groups [15]. Based on our data, we assumed that the downregulation of 2 fibroblast growth factor receptors (FGFR2, FGFR3) probably resulted from promoter hypermethylation of the coding genes. In the literature, downregulation of the remaining genes (CDH13, DST, ETS2, JAG1, ITGA2, PTGS1, TIAM1) have not been mentioned due to promoter hypermethylation so far.
Even with no detection of direct relationship of ulceration and promoter methylation, the above mentioned inverse correlations support the claim that similarly to copy number variation, promoter hypermethylation also plays an important role in transcriptomic silencing. Notably, it was suggested in the literature that Knudson's two-hit hypothesis is often achieved through a combination of DNA methylation and copy number alteration to the same gene [62,63]. Due to limitation of our studies, parallel detection of both types of somatic alterations was not achieved in primary melanomas. Nevertheless, the inverse relation between gene expression and methylation of the corresponding genes, together with global presence of copy number alteration should be taken into consideration, which suggests that somatic aberrations do not act separately but represent different utility in an integrated apparatus that acts together and develops transcriptomic silencing in ulcerated melanomas.
In summary, this study has provided evidence that the gene expression signatures of primary melanomas are suitable for distinguishing patients with poor and favourable prognoses. We have also shown that different patterns of genetic and epigenetic aberrations associated with distinct molecular subtypes of the disease contribute to the specific transcriptomic profiles of these genes. We believe that our systematic correlation of gene regulation, methylation pattern and the DNA copy number alteration data in the same cohort of primary melanomas will be useful for finding new genes and will allow further functional utilisation of those genes. This study has important implications as we continue to develop a better understanding of the metastatic process associated with melanomas, which will allow us to identify specific genes for use as prognostic markers and, possibly, for targeted therapeutic approaches. Figure S1 Cell adhesion and Cadherin signalling network. The network was significantly enriched of those genes which decreased expression is correlated with ulceration (p,0.001). Red circles indicate the downregulated genes. (TIF)

(DOC)
Text S1 High quality total RNA and DNA isolation, quality controls.