This article is mentioned in:

Introduction

Cutaneous melanoma is the most type of aggressive skin cancer due to its ability to rapidly metastasize and its radio- and chemotherapy resistance. Therefore, cutaneous melanoma is associated with a high mortality rate if it is not treated at early stages (1). Furthermore, the frequency of cutaneous melanoma has increased fourfold worldwide over the last 30 years (2). Melanoma may develop directly from transformed melanocytes or through a stepwise progression after first forming a benign nevus (3,4). Findings show that both paths result from mutations in genes of the mitogen-activated protein (MAP) kinase pathway and subsequent changes in other genes that affect regulatory pathways, including p53 (5) and phosphoinositide 3-kinase/AKT signaling pathways (6).

Melanocytes, which generate melanin in mammals (7), are uniformly distributed in the skin. These cells may form nevi, which are melanocyte clusters in specific areas. In most cases, nevi are caused by BRAF mutations or mutations in genes of the RAS family (8,9). These mutations constitutively activate the MAP kinase pathway, which is involved in cell division. The resultant nevus consists of melanocytes with limited growth that are in a ‘senescence-like’ state (10). Additional mutations in genes, such as cyclin dependent kinase inhibitor 2A (11) or phosphatase and tensin homolog (5), are known to be required for a cell in a nevus to develop into melanoma. This mechanism is observed in ~30% of melanoma cases (12).

Mutations may alter gene expression. Therefore, a variety of studies have addressed the classification and staging of melanoma through transcriptome analysis such as microarray technology. Some of those gene expression microarray studies on cancer transformation compare melanoma biopsies with skin as ‘normal’ sample or with nevus (13,14). However, strategies based on a meta-analysis that include a collection of such studies have proven to be useful identifying genes and pathways not only involved in the transformation but also in describing disease progression through gene profiling (15).

Melanoma is considered one of the most aggressive malignancies, and as such, patients with advanced melanoma tend to have a poor prognosis (1), and there is an urgent need to identify biomarkers. The work presented here proposes a list of biomarkers associated with malignant transformation in melanoma based on a transcriptome meta-analysis. Our approach to reveal these biomarkers was to compare the gene expression profiles of primary melanocyte cultures, nevus samples, and stages I, II, and III of melanoma samples according to the tumor-node metastasis (TNM) staging system, seventh edition (1) [primary tumor (T), regional lymph nodes (N) and distant metastasis (M)]. Microarrays downloaded from Gene Expression Omnibus repository (GEO) (16) were carefully selected. The present study focused on two paths: Malignant transformation directly from melanocytes and transformation through an intermediary nevus. By analyzing data carefully collected from different studies, the aim of the present study was to identify novel mechanisms underlying the molecular changes observed during malignant transformation. Outcomes of the study involved the identification of new pathways and biomarkers with potential application in clinical practice.

Materials and methods

Data extraction

Gene expression data obtained from a collection of different microarrays studies conducted prior to March 2015 were used for the present study. All data are publicly available through the GEO (16), which incorporated samples from melanocytes [accession no., GSE38312 (17), n=5], nevi [accession no., GSE53223 (18), n=5], and melanoma stage I, II and III tumors [accession no., GSE15605 (19), n=40].

Inclusion criteria for microarray samples

The inclusion criteria for microarray samples are as follows: i) Human biopsies; ii) cutaneous melanoma; iii) samples corresponding to primary melanoma sites; iv) samples classified as stages I, II or III according to TNM; v) nevus samples taken from biopsy; vi) melanocyte samples from primary cultures.

Exclusion criteria

The exclusion criteria are as follows: i) Samples from patients already under treatment; ii) metastatic samples; iii) sub-stage IIIC samples equivalent to in-transit metastasis due to an already intermediate path between the primary site and the invasive stage; iv) cell lines; v) melanoma in situ; vi) dysplastic nevi; vii) nevi in patients with a family history of melanoma (as these patients have a higher risk for melanoma). The selection of data files used in the present study is depicted in Fig. 1.

Figure 1. Selection of gene expression datasets for meta-analysis. The following terms were used for the initial search: ‘Homo sapiens’, ‘expression transcription profiling’, ‘array assay’ and ‘melanoma’. GEO, Gene Expression Omnibus.

Quality assessment and data extraction

All statistical analyses were conducted using the statistical software R (20) and a variety of libraries from Bioconductor (21). Raw data were background corrected using Robust Multiarray Average (22) and were normalized with quantile normalization (23) using the Affy package. As all microarrays involved were processed at raw level, no further correction for batch effect was required. Any batch effect was addressed during the normalization steps, which ultimately aims at making all samples comparable and at removing systematic variation.

Statistical analysis

Differential expression was analyzed using linear models and the Limma library (24). Two main summary statistics were used to classify genes as being differentially expressed: Fold-change on a log2 scale and B-statistic, which is the Bayesian equivalent to an adjusted P-value. According to the design of the present study, the analysis focused on 3 comparisons: i) Nevi vs. melanocytes; ii) melanoma vs. nevi and iii) melanoma vs. melanocytes. Target genes were selected based on a log fold-change >1 for upregulation and log fold-change <-1 for downregulation. Statistical significance was defined as a B-statistic >3, which corresponds to a 95.3% confidence level. Particular interest was focused on the gene set at the intersection of two comparisons: Melanoma vs. nevi and melanoma vs. melanocytes; as these genes represent possible biomarkers for melanoma transformation.

Unsupervised hierarchical clustering

Results from two comparisons of interest were used: Melanoma vs. nevi and melanoma vs. melanocytes. The two lists of differentially expressed genes were further reduced using Median Absolute Deviation (25) to capture the most variable set of genes. Using an unsupervised hierarchical clustering approach with the R-library ggplots (26), clustering patterns of genes (a) at the intersection of the two lists and (b) disjoint genes were searched for.

Gene Set Enrichment Analysis (GSEA)

GSEA was used to determine whether an a priori defined set of genes derived from annotation or pathway databases exhibited statistically significant differential expression between two models. The analysis was conducted using the Ingenuity® Pathway Analysis tool (IPA®, Qiagen, Inc., Valencia, CA, USA). The following thresholds were used for gene selection: −log(P-value) >1.3 (P-value <0.05) and a |z-score|>1. P<0.05 and a |z-score|>2 were used for analysis of biological processes.

Results

Eligible studies for the meta-analysis

The initial search yielded 315 datasets. Next, by screening the full characteristics and the quality control, 312 datasets were excluded. Therefore, only 3 datasets from GEO were selected, with samples from melanocytes [GSE38312 (17), n=5], nevi [GSE53223 (18), n=5], and melanoma stages I, II and III [GSE15605 (19), n=40]. Primary melanomas at stage IV and metastatic melanoma samples were excluded to simplify the analysis of gene expression changes that occur during the early transformation of non-cancerous cells (melanocytes or nevi) to melanoma cells. Samples that are clinically identified to have metastatic capacity according to staging were not included. By the same reasoning, sub-stage IIIC samples equivalent to in-transit metastasis were excluded. A flow chart of the selection process is illustrated in Fig. 1. The samples of each microarray in the meta-analysis are described in Table I.

Table I. Description of studies included in the meta-analysis.

Table I.

Description of studies included in the meta-analysis.

Accession no.	Year of publication	Summary	Total samples	Type of sample selected	Number of sample selected	(Refs.)
GSE38312	2014	Early-passage (<10 passages) cultures of melanoma cells from metastatic lymph node lesions and normal adult melanocytes explanted in parallel from the adjacent, non-involved skin of 5 patients	10	Melanocytes	5	(17)
GSE53223	2016	Excisional biopsies consisting of 5 common melanocytic nevi, 7 dysplastic nevi and 6 normal skin	18	Nevi	5	(18)
GSE15605	2013	Excisional biopsies consisting of 46 primary melanomas, 12 melanoma metastases and 16 normal skin	74	Melanoma stages I, II and III	40	(19)

Differential expression

Differentially expressed genes were identified in three comparisons: i) Nevi vs. melanocytes (2,492 genes); ii) melanoma vs. nevi (1,881 genes); and iii) melanoma vs. melanocytes (4,112 genes). The proportion of genes that are exclusive to each comparison and those that overlap is indicated in the Venn diagram in Fig. 2A. These comparisons explain the transformation from melanocytes to melanoma through two paths. One path represents a direct transformation process, with 1,385 downregulated and 946 upregulated genes (which is exclusive to this comparison), and another path develops from an intermediary nevus with 360 downregulated and 528 upregulated genes (which is also exclusive to this comparison) (Fig. 2B). It is important to note that even though differentially expressed genes obtained from the comparison of nevi vs. melanocytes were of no direct use, they were required to filter out those genes at the different intersections. For instance, the melanoma vs. nevi comparison with 1,881 differentially expressed genes comprises only 888 unique genes, because 311 genes were also present in the progression from melanocytes to nevi. Therefore, these genes must be excluded from the list, since they represent false positives in terms of explaining malignant transformation. Genes at the intersection of melanoma vs. nevi with melanoma vs. melanocytes (682 genes) were analyzed as potential biomarkers, as described below, since these are associated with melanoma development either directly from melanocytes or from nevi.

Figure 2. Differentially expressed genes. (A) Venn diagram of the differential expression results for the comparisons among nevi, melanocytes and melanoma Green circle, the number of genes with different expression levels between nevi and melanocytes. Yellow circle, the number of genes with different expression levels between melanoma vs. nevi. Orange circle, the number of genes with different expression levels between melanoma vs. melanocytes. Colors in the Venn diagram match the color of each comparison shown in B. (B) Diagram illustrating how the expression of genes change as they follow a path that ultimately results in melanoma. Blue and red lines indicate downregulation and upregulation, respectively. Black lines indicate no change. Genes in samples that developed melanoma directly from melanocytes are shown in the top panel, and the number of genes that was used to assess changes is indicated. Genes in samples in which melanoma developed from an intermediary nevus are shown at the bottom, and the number of genes that was used to assess changes from nevus to melanoma are indicated. As indicated in this diagram, 3 possible paths were identified in melanoma samples that had directly transformed from melanocytes, and 9 possible paths were identified in melanoma samples that transformed via an intermediary nevus. A given gene may have 10 possible transformation paths at the transcriptome level when removing the 2 paths that did not show changes in expression. (C) Hierarchical clustering based on the 50 most discriminating genes that are unique to each comparison: Melanoma vs. melanocytes (2,331 genes) and melanoma vs. nevi (888 genes). Melanocyte (light green label), nevus (dark green label) and melanoma (yellow label) samples fell into two major groups. (D) Hierarchical clustering based on the 50 most discriminating genes that are common to each comparison: Melanoma vs. melanocytes and melanoma vs. nevi (682 genes). Melanocyte (light green label), nevus (dark green label) and melanoma (yellow label) samples fell into two major groups.

In a search for a minimum gene set that could discriminate between the 3 groups (melanocytes, nevi and melanoma), the top 50 genes following median absolute deviation were selected. This procedure was followed by an unsupervised hierarchical clustering analysis. The results from this analysis are illustrated in Fig. 2C, which represents a heat map of the 50 most discriminating genes associated with cancer that is based on the comparisons: Melanoma vs. melanocytes and melanoma vs. nevi (2,331 and 888 genes, respectively; Fig. 2A). These top 50 genes were able to identify a transition profile from melanocytes to melanoma using nevi as intermediates. The clustering patterns for the genes at the intersection of melanoma vs. melanocytes and melanoma vs. nevi (682 genes, Fig. 2A) are represented by the heat map in Fig. 2D. The 50 genes appearing in both comparisons were non-randomly partitioned into two major groups and clearly describe malignant transformation, as melanocytes and nevi were clearly separated from melanoma.

Use of gene interaction enrichment and network analysis to identify dysregulated pathways

To interpret gene expression data based on the functional annotation of the differentially expressed genes and their association with biological processes or molecular functions, an enrichment analysis of genes that resulted from the malignant transformation was performed through nevus analysis as well as to the direct transformation obtained from the melanocyte analysis. The results are indicated in Fig. 3.

Figure 3. Gene Set Enrichment Analysis. (A) Signaling pathways and biological functions that are involved in cancer transformation via a nevus; three of the pathways exhibiting decreased activity are involved in cell cycle checkpoints. Conversely, a number of the activated pathways are associated with cell cycle progression. (B) Signaling pathways and biological functions that are involved in cancer transformation from melanocytes. The samples exhibited decreased cAMP synthesis, which is consistent with the decreased activity of this pathway. Numbers in parenthesis represent statistical significance (P-values) as assigned by Ingenuity® Pathway Analysis. BF, biological functions; PPARα, peroxisome proliferator activated receptor α; SP, signaling pathways.

Malignant transformation via nevi

Malignant transformation through an intermediary nevus was characterized by a total of 1,881 differentially expressed genes (comprising 737 downregulated and 1,144 upregulated genes). Only the 888 unique genes based on the Venn diagram in Fig. 2A were used to explain malignant transformation via nevi. An enrichment analysis of these genes identified 13 dysregulated pathways; 6 of these (mostly associated with cell cycle control) exhibited decreased activity, and 7 (associated with cell cycle progression) exhibited increased activity. In addition, 6 biological functions were activated, and 4 of these were associated with the hallmarks of cancer: Proliferation, migration, cell movement and autophagy (Fig. 3A).

Direct transformation from melanocytes

The model based on direct transformation from melanocytes included 4,112 differentially expressed genes (2,310 downregulated and 1,802 upregulated). However, only the 2,331 genes unique to malignant transformation from melanocytes were used to perform the gene enrichment analysis for signaling pathways and biological processes. The results show that specifically cAMP-mediated signaling was the only pathway to exhibit decreased activity. By contrast, pattern recognition receptors (C1QA, C1QB, C1QC, CCL5 and IFNA4) were upregulated. Biological functions associated with the phosphorylation of L-amino acids and cAMP synthesis exhibited decreased activity (Fig. 3B).

Potential biomarkers

To identify transformation biomarkers in melanoma, the present study focused on the genes at the intersection of melanoma vs. melanocytes and melanoma vs. nevi comparisons (682 genes, Fig. 2A). The analysis began with a classification of these genes into intracellular and extracellular biomarkers with the aim of capturing effects at both the transcriptional and protein levels. Genes that appeared upregulated in both comparisons and whose product was a secreted protein were selected as potential secreted biomarkers (22 genes) and are listed in Table II. Previously characterized cancer biomarkers, cytokines C-X-C motif chemokine ligand (CXCL)8 and CXCL16, were included in this group. Another biomarker, that was included in the list, matrix metallopeptidase 19 (MMP19), is inhibited by the drug Marimastat.

Table II. Potential secreted biomarkers associated with melanoma transformation.

Table II.

Potential secreted biomarkers associated with melanoma transformation.

			Differential expression levels
Symbol	Gene name	Cancer-associated functions	Melanoma vs. nevi	Melanoma vs. melanocytes	(Refs.)
C10orf95	Chromosome 10 open reading frame 95	NR	54.57	64.89
ADAMDEC1	ADAM-like, decysin 1	Tumor progression in craniopharyngioma	14.72	24.76	(27)
IK	IK cytokine, down-regulator of HLA II	Immune response inhibition in cutaneous T cell lymphomas	14.52	12.64	(28)
NENF	Neudesin neurotrophic factor	Breast tumorigenesis, growth and metastasis in liver cancer	11.24	9.25	(29,30)
MMP19	Matrix metallopeptidase 19	Melanoma development and migration of tumor cells	8.88	6.63	(31)
TNFSF13	Tumor necrosis factor superfamily member 13	CLL leukaemogenesis	8.88	8.51	(32)
APOL6	Apolipoprotein L6	Inducer of apoptosis in cancer cell lines	6.59	7.67	(33)
TFIP11	Tuftelin interacting protein 11	NSCLC progression	6.45	6.06	(34)
ASXL2	Additional sex combs like 2, transcriptional regulator	Prognostic marker in breast cancer	5.98	3.58	(35)
WDR73	WD repeat domain 73	NR	5.94	7.73
CXCL8	C-X-C motif chemokine ligand 8	Marker for urinary bladder cancer, autophagy biomarker and metastatic potential in melanoma	5.66	5.66	(36,37)
BIVM	Basic, immunoglobulin-like variable motif containing	NR	5.62	4.41
TCTN2	Tectonic family member 2	Poor clinical outcome in glioblastoma	5.62	3.34	(38)
ZFC3H1	Zinc finger C3H1-type containing	Regulation of CXCL8 transcription in cancer cell lines	5.46	5.70	(39)
CXCL16	C-X-C motif chemokine ligand 16	Invasion and migration of prostate cancer cells	4.56	13.64	(40)
FAM3C	Family with sequence similarity 3 member C	Phenotype switching in melanoma cells	4.17	5.31	(41)
FAM49B	Family with sequence similarity 49 member B	Potential biomarker for oral cancer	3.61	3.68	(42)
RSPH10B2	Radial spoke head 10 homolog B2	NR	3.56	5.62
LGALS8	Lectin, galactoside binding soluble 8	Proliferation and survival of glioblastoma cells	3.16	4.08	(43)
TXLNA	Taxilin α	Potential biomarker for aggressive behavior in RCC	3.12	3.05	(44)
PLTP	Phospholipid transfer protein	Potential biomarker for higher grade of glioma	2.87	4.59	(45)
FAM60A	Family with sequence similarity 60 member A	Cell cycle regulator upregulated in NSCLC	2.81	6.59	(46)

[i] NR, not reported; CLL, chronic lymphocytic leukemia. NSCLC, non-small cell lung cancer; RCC, renal cell carcinoma.

There were more intracellular biomarkers in the list, which might be used as cancer biomarkers in biopsies. In total, there were 48 genes comprising 24 downregulated and 24 upregulated genes. Literature reviews on each gene consistently identify these biomarkers in human pathology (9,47–78). The list is shown in Table III.

Table III. Potential intracellular biomarkers associated with melanoma transformation.

Table III.

Potential intracellular biomarkers associated with melanoma transformation.

			Differential expression levels
Symbol	Gene name	Cancer-associated functions	Melanoma vs. nevi	Melanoma vs. melanocytes	(Refs.)
KRAS	Kirsten rat sarcoma viral oncogene homolog	Oncogene	−6.96	−10.56	(9)
STAT5B	Signal transducer and activator of transcription 5B	Favorable prognostic marker in ovarian cancer	−4.59	−4.92	(47)
GRPR	Gastrin releasing peptide receptor	Angiogenesis during neuroblastoma progression	−3.73	−5.66	(48)
ANTXR1	Anthrax toxin receptor 1	Tumor angiogenesis, including melanoma	−3.73	−5.66	(49)
CACNA1G	Calcium voltage-gated channel subunit α1 G	Tumor suppressor, low expression in numerous cancers	−3.73	−4.59	(50)
PDE5A	Phosphodiesterase 5A	Regulation of melanoma cell invasion	−3.73	−4.59	(51)
RTN4	Reticulon 4	Proapoptotic in cancer cell lines	−3.48	−4.29	(52)
NR4A1	Nuclear receptor subfamily 4 group A member 1	Marker for survival in lung tumors	−2.83	−3.48	(53)
C1QTNF5	C1q and tumor necrosis factor related protein 5	NR	−2.83	−3.48
TCF7L2	Transcription factor 7 like 2	Tumor growth in skin squamous cell carcinoma	−2.83	−6.06	(54)
CLCNKB	Chloride voltage-gated channel Kb	Underexpressed in kidney tumors	−2.83	−3.25	(55)
DRD2	Dopamine receptor D2	NR	−2.64	−3.48
SEPT9	Septin 9	Biomarker of colorectal cancer	−2.64	−2.46	(56)
DISC1	Disrupted in schizophrenia 1	NR	−2.64	−2.00
CDH6	Cadherin 6	Tumor suppressor in cholangiocarcinoma	−2.64	−2.46	(57)
AQP2	Aquaporin 2	NR	−2.46	−2.64
TUBA1A	Tubulin α1a	NR	−2.46	−2.83
MYOD1	Myogenic differentiation 1	Regulation of gastric cancer cell apoptosis	−2.46	−2.46	(58)
TPM1	Tropomyosin 1 (α)	Radiation resistance in glioblastoma cells	−2.46	−4.59	(59)
BCAR1	BCAR1, Cas family scaffolding protein	Growth and invasion of several types of cancer cells	−2.30	−4.59	(60)
MLH1	MutL homolog 1	Risk factor of colorectal cancer	−2.30	−2.46	(61)
UBE2I	Ubiquitin conjugating enzyme E2 I	NR	−2.14	−3.03
SARNP	SAP domain containing ribonucleoprotein	NR	−2.14	−2.64
TNNT2	Troponin T2, cardiac type	NR	−2.00	−2.30
TCF3	Transcription factor 3	NR	26.72	82.14
RANBP1	RAN binding protein 1	NR	15.24	14.12
KIF14	Kinesin family member 14	Tumor progression and poor prognosis in prostate cancer	14.22	5.78	(62)
SOCS3	Suppressor of cytokine signaling 3	Anti-tumor effect in lung cancer and gastric tumors	14.12	5.86	(63)
FGFR2	Fibroblast growth factor receptor 2	Poor prognosis biomarker in most solid tumors	12.91	11.00	(64)
CTLA4	Cytotoxic T-lymphocyte associated protein 4	Poor prognosis biomarker in melanoma	11.79	14.72	(65)
CXCR4	C-X-C motif chemokine receptor 4	Prognostic biomarker for stage II primary melanomas	10.85	257.78	(66)
TXN	Thioredoxin	NR	9.85	10.93
TP53	Tumor protein p53	Strong protein expression in melanoma tissue indicates aberrations of p53 activity	8.46	4.69	(67)
LMCD1	LIM and cysteine rich domains 1	Cell migration and tumor metastasis in hepatocellular carcinoma	6.96	8.34	(68)
B4GALT1	β-1,4-galactosyltransferase 1	Development of multidrug resistance of leukemia cells	5.86	7.94	(69)
EIF4E	Eukaryotic translation initiation factor 4E	Melanoma proliferation and invasion	5.31	2.97	(70)
TNFRSF4	Tumor necrosis factor receptor superfamily member 4	NR	4.69	4.69
HINT1	Histidine triad nucleotide binding protein 1	Tumor suppressive activity	4.66	4.66	(71)
CPT1B	Carnitine palmitoyltransferase 1B	NR	4.59	5.62
GNA12	G protein subunit α12	Cell invasion in prostate cancer cells	4.56	4.00	(72)
STK11	Serine/threonine kinase 11	Tumor invasion and migration	4.26	5.21	(73)
GAA	Glucosidase α, acid	NR	4.06	4.44
CTNNB1	Catenin β1	Promoter of melanoma development	3.84	2.53	(74)
PRKCI	Protein kinase C ι	Tumor progression in some lung cancer	3.41	5.10	(75)
TSNAX	Translin-associated factor X	NR	2.87	4.38
GLI3	GLI family zinc finger 3	Tumorigenesis of glioblastoma and medulloblastoma	2.77	3.27	(76)
NCOA3	Nuclear receptor coactivator 3	Tumorigenesis in multiple types of cancer	2.35	3.18	(77)
USP9X	Ubiquitin specific peptidase 9, X-linked	Promotion of glioma cell proliferation and survival	2.33	2.79	(78)

[i] NR, not reported.

Discussion

In this meta-analysis, which included 3 independent datasets with 50 gene expression microarrays from GEO, the gene expression of two molecular pathways for melanoma were compared. One pathway transforms melanocytes directly to melanoma and another pathway that passes through an intermediate state (nevi). Data were collected based on samples taken from melanocytes, nevi, which were taken as an intermediary stage, and melanoma stages I, II and III.

A number of previous studies compared samples in the radial growth phase vs. samples in the vertical growth phase (79), whereas other studies analyzed melanoma but without considering clinical stage (13). A number of studies also used cell lines, nevi or skin as ‘normal’ samples (14). In the present study, the selection of melanocytes as the ‘normal’ sample and the selection of biopsies obtained only from primary melanomas (stages I, II or III) are major strengths of our work.

The results show clear differences between the two pathways of melanoma development according to the cell of origin. The transformation through nevi model shows that the list of differentially expressed genes are associated with biological functions that are consistent with the development and progression of cancer, including proliferation, migration, cell movement (80) and autophagy (81). In addition, three pathways associated with cell cycle control that exhibited declined activity were identified: Cell cycle checkpoint control involving CHK proteins, DNA damage checkpoint regulation at G2/M and regulation of G1/S checkpoint. Another pathway exhibiting a decreased activity was peroxisome proliferator activated receptor α/RXRα activation, which has been reported to exert anti-proliferative effects on melanocytes and melanoma (82,83). By contrast, pathways with increased activity are cyclins and cell cycle regulation, and estrogen-mediated S-phase entry, which promote cell cycle progression. These pathways are consistent with malignant transformation.

Shepelin et al (15) recently identified 25 metabolic pathways and 19 signaling pathways as being implicated in the transformation to melanoma using skin samples as the starting stage and nevi as an intermediary stage. An increased function of pathways was reported, including hypoxia inducible factor 1 α subunit, ATM (DNA repair), putrescine biosynthesis II, and breast and ovarian cancer susceptibility protein 1. A decreased activity was reported for biosynthesis of histamine, eumelanin and allopregnanolone and degradation of valine. However, these results were not consistent with our findings.

Melanocytes represent only a minority in the cell population within the basilar epidermis since the melanocyte to keratinocyte ratio is 1:5 and they are neural crest-derived cells (84). The present analysis included transcriptomes from primary cultures of melanocytes whilst previous studies have used skin as ‘normal’ samples (13,14). On the other hand, continuous cell lines have already evaded senescence control as they have increased cell proliferation rate. The authors consider that finite cultures of melanocytes represent a better biological model since these cells are manipulated in vitro for shorter periods (<10 passages) (17). We also consider that this is the main reason why the results showed a large number of differentially expressed genes. Enrichment analysis only showed that the decreased activity of the cAMP-mediated signaling pathway was significantly associated with transformation. Previous studies report that the cAMP signaling pathway inhibits the development of melanoma (85); this finding contrasts with another report that claims that this pathway promotes the growth of melanoma (86).

So far, lactate dehydrogenase (LDH) has been useful as a secretory biomarker in the tumor-node metastasis (TNM) staging of melanoma. High levels of LDH are associated with poor prognosis at advanced stages of the disease (M1c, stage IV) (1). In addition to the use of LDH as a biomarker, the identification of more molecules as good biomarkers of melanoma is necessary to identify not only stage IV disease but also early stages of the disease and high-risk populations. The present meta-analysis identified 22 potential secreted biomarkers, of which CXCL16 and CXCL8 have been previously reported as biomarkers in patients with Crohn's disease (87) and in several diseases, including urinary bladder cancer, non-Hodgkin's lymphoma and pulmonary infections (88). A literature review indicated that CXCL8 is linked to autophagy in melanoma cell lines (35) and has been reported in metastasis (37). Here, CXCL8 appeared in the top 50 of modulated genes (Fig. 3C). Another potential secreted biomarker, MMP19, has been identified in progression and invasiveness events in melanoma (31). To examine the basal abundance of these proteins, the Protein Atlas (89) was searched, and the search identified the levels of CXCL8, CXCL16, MMP19, basic, immunoglobulin-like variable motif containing (BIVM) and family with sequence similarity 60 member A (FAM60A) to be low in other tissues. BIVM is a poorly studied protein, whereas FAM60A is a cell cycle regulator that is reportedly upregulated in non-small cell lung cancer (46). These results are interesting because low abundance of CXCL8, CXCL16, MMP19, BIVM and FAM60A in other tissues might correspond to low serum concentrations in healthy individuals. Therefore, if these proteins were abundant in melanoma, this would facilitate their detection in the serum as biomarkers. Other genes that encode secreted proteins include ADAM like decysin 1 (ADAMDEC1), neudesin neurotrophic factor (NENF) and apolipoprotein L6 (APOL6), which are included in the list of top 50 modulated genes (Fig. 3D). In brief, ADAMDEC1 encodes a protein with metalloendopeptidase activity and has been identified as a biomarker of recurrence and poor prognosis in craniopharyngioma (28). NENF has been reported as an oncogene in breast cancer (29) and as a potential therapeutic target in liver cancer (30). APOL6 encodes an apolipoprotein BCL-2 homolog that acts as an apoptosis inducer (33).

Several tentative intracellular biomarkers were also identified in the analysis, 24 of which were downregulated and 24 of which were upregulated. Among previously reported genes, it was observed that tumor protein P53 (TP53) is upregulated in melanoma compared with melanocytes and nevi. The upregulation of TP53 in melanoma has been reported, albeit with abnormal function and regulation (67). Other findings include MDM2, which was upregulated in melanoma compared with melanocytes, and MDM4, which was upregulated in melanoma compared with nevi. Both MDM2 and MDM4 are inhibitors of p53. Similarly, the downregulation of various inhibitors of cyclin-dependent kinases [cyclin dependent kinase inhibitor (CDKN)1A, CDKN1C, CDKN2C and CDKN2D], which interact with p53 or inhibit the activity of various cyclin-dependent kinases, was observed. It was also observed that C-X-C motif chemokine receptor 4 (CXCR4), a chemokine receptor involved in metastasis, was upregulated. CXCR4 has a prognostic value in stage II melanoma (66). Another potential intracellular biomarker is cytotoxic T-lymphocyte associated protein 4 (CTLA4), which is upregulated during the transformation to melanoma, and in the list of the top 50 modulated genes (Fig. 2D). This protein has gained increasing attention in recent years as the target of ipilimumab, the first immune checkpoint inhibitor approved to treat unresectable melanoma (90). CTLA-4, a CD28 homologue, was originally described as a glycoprotein that is membrane bound on activated T-cells and regulatory T (Treg)-cells and which has a high binding affinity with the molecule B7, which is primarily expressed on antigen-presenting cells (APCs) (91). The expression and function of CTLA4 have been previously observed for various types of non-T cells, including several malignant solid tumors such as melanoma (92). This suggests that CTLA-4 might be involved in controlling functions other than the widely described T-cell response inactivation (93). In the tumor microenvironment, it is well documented that CTLA4 is upregulated in activated T cells, therefore providing inhibitory co-stimulation; it may lead to apoptosis and/or a state of immunologic anergy. Treg cells, which constitutively express high levels of CTLA4, also regulate effector T cells by downregulating B7 in APCs, leading to the induction and maintenance of T-cell tolerance by reducing CD28 co-stimulation (94). In addition, B7 ligation by CTLA-4 on Treg cells induces immunosuppressive APC by enhancing the indoleamine 2,3-dioxygenase tryptophan catabolic enzyme cascade, thereby preventing T-cell activation by depriving T cells of tryptophan. Proapoptotic tryptophan catabolites can also exert strong immunosuppression and indirectly suppress effector T cell response by favoring differentiation of Treg cells, which in turn would repeat CTLA-4 interactions (95). For all those reasons, CTLA-4 is considered the ‘leader’ of the immune checkpoint. At present, the physiological role of CTLA4 expression in melanoma cells remains unknown, but a recent study reported that high levels of CTLA4 in tumor cells in melanoma patients were associated with poor clinical outcomes (65). Other studies have demonstrated that the upregulation of CTLA4 is a factor predictive of poor prognosis in nasopharyngeal (96), esophageal (97) and breast carcinomas (98).

In this context, it is worth noting that even soluble forms of CTLA-4 (99) and cell transfectants that upregulate membrane-anchored CTLA-4 (100) can signal through B7 in a similar manner to CTLA-4 on T cells and resemble the same tolerogenic effects. Indeed, breast cancer cells that express CTLA-4 are involved in this novel tumor-evading mechanism (93). A number of types of cancer cells express CTLA4 (92), so it is likely that in order to evade immune surveillance the same mechanism is used by melanoma cells.

We would also like to address some limitations of the present study. First of all, there are a small number of studies with TNM staging and therefore, the number of included microarrays was not as large as one would expect in a meta-analysis. However, statistical power might be undermined by a large collection of heterogeneous samples than it could be for a small but well characterized set of samples. Secondly, the data from melanocytes comes from a monoculture setting. It is well known that keratinocytes significantly impact the mRNA expression profile of melanocytes; therefore, the profile may not be the same for melanocytes in situ. Primary culture is the best resembling natural tissue but is largely heterogeneous. We believe that this is the main reason why transition from melanocyte to melanoma showed only 2 dysregulated signaling pathways based on the enrichment analysis of genes identified as differentially expressed. Lastly, there are technical biases that are due to the fact that the samples were processed in different labs, by different people and at different times. This was well addressed during the normalization process and batch effect was removed by pre-processing all samples at raw level.

In conclusion, the present study provided a list of differentially expressed genes, molecular pathways and potential biomarkers associated with melanoma for direct transformation from melanocytes and from nevus as intermediary stage. The data suggest that molecular mechanisms underlying these two pathways of transformation into melanoma are different; however, both pathways share molecular changes that might be important for the identification of early biomarkers that are specific to melanoma. Overall, testing for the high-level expression of CXCL8, CXCL16, MMP19, BIVM, FAM60A, ADAMDEC1, NENF and APOL6, whose products are secreted proteins, might detect malignant transformation at its earliest presence, a requisite for improving preventive interventions and reducing cancer incidence. Other intracellular biomarkers found might have prognostic or therapeutic target potential. The present meta-analysis found in silico evidence that upregulation of CTLA4 may have a role in the transformation towards melanoma, possibly related to immune tolerance.

Acknowledgements

Not applicable.

Funding

The present study was supported by grants from Consejo Nacional de Ciencia y Tecnología (CONACYT), which was award to Claudia H. González-De la Rosa (grant no. 181534) and Consejo Nacional de Ciencia y Tecnología (CONACYT), which was awarded to Daniel Ortega Bernal (fellowship no. 372614).

Availability of data and materials

The datasets analyzed during the current study are available in the Gene Expression Omnibus (GEO) database repository, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38312, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53223 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE15605.

Authors' contributions

DO, CG and CR were responsible for the conception and design of the study. DO, EA and CR developed themethodology, and DO was responsible for the acquisition of data. DO, CG, EA, MA, NS and CR were involved in the analysis and interpretation of data. DO, CG, EA, MA, NS and CR wrote, reviewed and revised the manuscript. CG and CR were responsible for the supervision of the study.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References