Skip to main content

Advertisement

MEK inhibition remodels the active chromatin landscape and induces SOX10 genomic recruitment in BRAF(V600E) mutant melanoma cells

Epigenetics & Chromatinvolume 12, Article number: 50 (2019) | Download Citation

Abstract

Background

The MAPK/ERK signaling pathway is an essential regulator of numerous cell processes that are crucial for normal development as well as cancer progression. While much is known regarding MAPK/ERK signal conveyance from the cell membrane to the nucleus, the transcriptional and epigenetic mechanisms that govern gene expression downstream of MAPK signaling are not fully elucidated.

Results

This study employed an integrated epigenome analysis approach to interrogate the effects of MAPK/ERK pathway inhibition on the global transcriptome, the active chromatin landscape, and protein–DNA interactions in 501mel melanoma cells. Treatment of these cells with the small-molecule MEK inhibitor AZD6244 induces hyperpigmentation, widespread gene expression changes including alteration of genes linked to pigmentation, and extensive epigenomic reprogramming of transcriptionally distinct regulatory regions associated with the active chromatin mark H3K27ac. Regulatory regions with differentially acetylated H3K27ac regions following AZD6244 treatment are enriched in transcription factor binding motifs of ETV/ETS and ATF family members as well as the lineage-determining factors MITF and SOX10. H3K27ac-dense enhancer clusters known as super-enhancers show similar transcription factor motif enrichment, and furthermore, these super-enhancers are associated with genes encoding MITF, SOX10, and ETV/ETS proteins. Along with genome-wide resetting of the active enhancer landscape, MEK inhibition also results in widespread SOX10 recruitment throughout the genome, including increased SOX10 binding density at H3K27ac-marked enhancers. Importantly, these MEK inhibitor-responsive enhancers marked by H3K27ac and occupied by SOX10 are located near melanocyte lineage-specific and pigmentation genes and overlap numerous human SNPs associated with pigmentation and melanoma phenotypes, highlighting the variants located within these regions for prioritization in future studies.

Conclusions

These results reveal the epigenetic reprogramming underlying the re-activation of melanocyte pigmentation and developmental transcriptional programs in 501mel cells in response to MEK inhibition and suggest extensive involvement of a MEK-SOX10 axis in the regulation of these processes. The dynamic chromatin changes identified here provide a rich genomic resource for further analyses of the molecular mechanisms governing the MAPK pathway in pigmentation- and melanocyte-associated diseases.

Background

The mitogen-activated protein kinase (MAPK)/extracellular signal-regulated kinase (ERK) pathway is an evolutionarily conserved signaling cascade that relays multiple signals from the plasma membrane to the nucleus to regulate cell growth, differentiation, proliferation, and survival. In MAPK/ERK pathways, ligand binding to a receptor tyrosine kinase at the plasma membrane leads to activation of the small GTPase RAS. This is followed by a cascade of complex protein interactions and phosphorylation events, including RAF phosphorylation and activation of MAPK kinase (MEK), which directly phosphorylates ERK and regulates its nuclear translocation and transcriptional activity [1, 2]. MAPK pathway deregulation due to genomic mutations is common in many human cancers, including thyroid cancer, colon cancer, serous ovarian cancer, and melanoma [3, 4]. Approximately 50% of melanomas carry the activating BRAF(V600E) mutation, which causes constitutive hyperactivation of the MAPK pathway [5]. Small-molecule inhibitors of BRAF such as vemurafenib and dabrafenib have shown remarkable results as molecularly targeted therapies for melanoma, especially when used in combination with MEK inhibitors, including the MEK inhibitor AZD6244 (also known as selumetinib or ARRY-142886) [6,7,8,9,10,11,12,13,14]. However, clinical responses to BRAF and MEK inhibitors are often not durable due to acquisition of resistance after treatment [15,16,17], highlighting the need to identify novel melanoma vulnerabilities that could be exploited for new treatments.

One approach to address this need is through analysis of epigenetic changes that arise from altered MAPK/ERK pathway signaling. Chromatin-mediated processes play crucial roles in the signaling cascades that govern development, cancer progression, and drug resistance [18,19,20]; however, little is known about the chromatin signatures and epigenetic effects caused by MAPK inhibitors. In particular, altered chromatin modifications at enhancers—non-coding DNA elements that interact with transcription factors and chromatin remodelers—have been implicated in aberrant modulation of gene expression and drug resistance in many cancers, including melanoma [21,22,23,24,25,26]. Enhancers can be identified through profiling of specific biochemical signatures of chromatin, such as accessibility to transposase digestion, hypersensitivity to DNase I cleavage, and acetylation of histone H3 on lysine residue 27 (H3K27ac, [27,28,29]). In addition, large clusters of enhancers called super-enhancers have been implicated in the regulation of normal development [30, 31] as well as cancer progression and resistance to therapy [24, 32, 33]. Super-enhancers can undergo remodeling during development, in response to disease, or as a result of functional genomic and pharmacological perturbations, and correspondingly super-enhancers regulate genes that govern cell identity, development, and pathogenesis [24, 25, 30,31,32, 34,35,36]. Super-enhancers exhibit an exceptionally high density of binding for numerous factors, including the multiprotein Mediator complex and BRG chromatin remodeling complex, transcriptional cofactors, and lineage-determining transcription factors, including members of the SOX transcription factor family in embryonic stem cells (SOX2), neural progenitor cells (SOX2), chondrocytes (SOX9, SOX5/6), and hair follicle stem cells (SOX9) [30, 31, 37,38,39].

The SOX family member SOX10 plays a crucial role in neural crest and melanocyte development as well as human disease. SOX10 mutations are associated with the human neurocristopathies Waardenburg Syndrome (types 2 and 4) and PCWH (OMIM#609136, #611584, and #613266, https://omim.org/), in which abnormal neural crest/melanocyte development results in phenotypes that include hypopigmentation and enteric ganglia defects. SOX10 predominantly occupies distal (non-promoter) enhancer regions in melanocyte and melanoma genomes, colocalizes with MITF at enhancers, and interacts with BRG1 at melanocyte-specific enhancer regions [18, 40,41,42,43,44]. Furthermore, SOX10 is required for melanoma initiation, survival and progression, and cellular levels of SOX10 affect acquired resistance to melanoma inhibitor drugs [45,46,47,48,49].

In this study, we used a BRAF(V600E) mutant melanoma cell line as a model system to evaluate genome-wide changes in gene expression, enhancer chromatin marks, and transcription factor occupancy in response to the MEK inhibitor AZD6244. AZD6244 treatment caused hyperpigmentation as well as widespread alterations in transcription and the active chromatin landscape that were associated with numerous pigmentation genes. Distinct transcriptional motifs underlying AZD6244-induced epigenome remodeling were discovered, including SOX, MITF, ETV, and ATF motifs. A novel cohort of super-enhancers was identified that was enriched in melanocyte lineage-specific transcription factor motifs. Extensive recruitment of SOX10 occurred throughout the genome under MEK inhibition, including increased interactions at H3K27ac-associated enhancers and super-enhancers, with notable binding of SOX10 at super-enhancers associated with AZD6244-upregulated pigmentation genes. Furthermore, many of these chromatin regions directly overlapped GWAS loci for human pigmentation phenotypes. These data provide a roadmap of enhancer activity in a BRAF(V600E) mutant cell line under MEK inhibition, suggest the chromatin regions identified in this study regulate pigment- and developmentally related programs of gene expression, and implicate SOX10 in the regulation of these responses to MEK inhibition.

Results

MEK inhibition induces pigmentation in melanoma cells

To assess the impact of MEK inhibition on the epigenome in melanoma cells, MEK inhibitor-responsive human 501mel cells (heterozygous for the BRAF(V600E) mutation) were treated with AZD6244 (MEK inhibitor) or DMSO (vehicle control). As expected, treatment of 501mel cells with AZD6244 for 48 h or 72 h eliminated ERK phosphorylation (Fig. 1a); however, the cells remained proliferative (data not shown). Interestingly, 501mel cells became hyperpigmented following AZD6244 treatment, as evidenced by the dramatic change in color of the cell pellets (Fig. 1b). Similar hyperpigmentation effects of AZD6244 treatment were observed in two additional melanoma cell lines (A375 and mouse B16F10 melanoma cell lines, Additional file 1: Figure S1). The emergence of AZD6244-induced hyperpigmentation in melanoma cells suggested that melanocyte lineage-specific pathways may be re-activated in response to MEK inhibition.

Fig. 1
figure1

Inhibition of MEK induces hyperpigmentation and widespread gene expression changes in 501mel melanoma cells. a ERK and phosphorylated ERK (p-ERK) Western blots of whole-cell lysates derived from 501mel cells treated with DMSO or AZD6244 (200 nM). Tubulin was used as a protein loading control. b Cell pellets from 501mel cells after 72 h of treatment with DMSO or AZD6244. c Volcano plot of genes with significantly altered expression levels under AZD6244 treatment. Downregulated genes are in green, and upregulated genes are in red; the genes with an absolute fold-change > 2.4 and a p-adjusted value > 0.05 are labeled. d Gene ontology analysis on the genes significantly upregulated and downregulated in response to AZD6244 treatment. The top enriched processes for upregulated genes affect developmental pigmentation and autophagy, while those for downregulated genes affect protein targeting to membranes and protein translation

Full size image

Gene expression profiling by RNA-Seq was performed to identify genes whose expression changed in response to AZD6244 treatment in 501mel cells. Significant expression changes (false discovery rate (FDR)-adjusted p value < 0.05) were identified for 5526 genes following AZD6244 treatment compared to control cells (Fig. 1c), with similar numbers of upregulated (n = 2615) and downregulated (n = 2911) genes (Additional file 2: Table S1). Gene ontology (GO) analysis (geneontology.org, [50]) found distinct differences between the upregulated and downregulated genes (Fig. 1d). The upregulated genes showed enrichment of developmental pigmentation and autophagy processes, consistent with the reappearance of AZD6244-induced pigmentation (Fig. 1b). Correspondingly, 250 of the 5526 differentially expressed genes were known to affect pigmentation phenotypes in human, mouse, or zebrafish ([51], Additional file 2: Table S2), suggesting widespread alteration of pigmentation pathways in response to AZD6244 treatment. These genes included two of the genes with the largest increase in expression, OCA2 and IRF4 (Fig. 1c), which affect albinism and normal human pigmentation variation [52,53,54,55]. In contrast, downregulated genes were enriched for processes affecting protein targeting to membranes and protein translation, suggesting widespread attenuation of these pathways. These processes were driven by a cohort of 81 genes that encompassed 40S and 60S ribosomal proteins, signal recognition particle pathway proteins, and ribosomal cofactors, consistent with MAPK pathways regulating the ribosomal biogenesis necessary for proliferation and survival [56].

MEK inhibition remodels the active chromatin landscape

To evaluate the broad chromatin changes underlying the cellular response to MEK inhibition, chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-Seq) was performed for H3K27ac, an established histone modification mark of transcriptionally active regulatory elements such as promoters and enhancers [28, 29]. Treatment of 501mel cells with AZD6244 for 72 h resulted in a decrease in the total number of H3K27ac-enriched regions (peaks) relative to DMSO-treated (control) cells; 39,832 and 24,135 significantly enriched H3K27ac peaks were detected in control and AZD6244 groups, respectively (FDR < 0.05, Fig. 2a). In both control and AZD6244-treated samples, the majority of H3K27ac-enriched regions were located outside of proximal promoter regions (Fig. 2b), consistent with the known enrichment of H3K27ac in distal enhancer regions [28, 29]. The genome-wide reduction of H3K27ac in response to AZD6244 treatment was also reflected in the complete loss of 22,634 H3K27ac peaks, as compared to 1989 novel, newly gained H3K27ac regions.

Fig. 2
figure2

MEK inhibition induces genome-wide modulation of the active chromatin landscape in 501mel cells. a Venn diagrams showing overlap of H3K27ac peaks identified by ChIP-Seq in independent replicate samples of 501mel cells treated with DMSO (upper panel) and AZD6244 (lower panel). b Distribution of H3K27ac peaks in 501mel cells treated with DMSO (gray) or AZD6244 (red) at either promoter (< 2.5 kb from TSS) or non-promoter, distal enhancer regions. c, f Heatmaps showing relative H3K27ac signal in 501mel cells treated with DMSO or AZD6244. H3K27ac ChIP-Seq tag densities are shown for significantly differentially remodeled H3K27ac sites (FDR < 0.05) at regions c losing H3K27ac or f gaining H3K27ac. Each row represents a 10-kb region centered on differentially remodeled H3K27ac regions, sorted by H3K27ac signal enrichment. d, g Line plots show average H3K27ac ChIP-Seq signal intensities (y-axis) across regions d losing H3K27ac or g gaining H3K27ac over a ± 5 kb genomic window centered on differentially remodeled H3K27ac sites (x-axis). e, h The top significantly enriched transcription factor motifs identified by de novo motif discovery of sequences in regions e losing H3K27ac and h gaining H3K27ac. The top 5 enriched motifs, the fraction of H3K27ac regions containing the motifs, and the p values for the overrepresentation of each motif are shown in Additional file 2: Table S3

Full size image

MEK inhibitor-induced, differentially enriched H3K27ac regions contain distinct transcriptional motif signatures

DiffBind [57] was used to quantitatively compare H3K27ac peaks that were present in both DMSO-treated and AZD6244-treated cells, but exhibited differential enrichment of H3K27ac levels following AZD6244 treatment. This analysis found that 17,517 regions across the genome exhibited significant differential enrichment of H3K27ac (FDR < 0.05). More than 76% of these regions were located outside the proximal promoters of their putative target genes (> 5 kb from transcriptional start sites, TSS), again consistent with H3K27ac enrichment at distal regulatory elements/enhancers. Functional enrichment analysis of the top ~ 20% of differentially enriched H3K27ac regions (|FC| > 1.456, n = 3359) using the Genomic Region Enrichment of Annotation Tool (GREAT, [58]) showed enrichment of biological processes related to growth, differentiation, and development. The vast majority of the 17,517 differently enriched regions showed significantly decreased H3K27ac (16,476 regions, Fig. 2c, d), while the remaining 1041 regions exhibited significantly increased H3K27ac (Fig. 2f, g). Interestingly, identification of differentially expressed genes (Additional file 2: Table S1) that are also putative targets of the 17,517 H3K27ac-marked enhancers (identified using GREAT) showed similar numbers of both up-and downregulated genes in enhancer regions that were losing H3K27ac (1623 upregulated and 1993 downregulated genes) and in regions gaining H3K27ac (350 upregulated and 189 downregulated genes). This suggests that H3K27ac loss does not exclusively cause downregulation of gene expression, nor does H3K27ac gain only cause upregulation of gene expression.

H3K27ac-enriched chromatin regions harbor active enhancer elements and binding sites for transcription factors [28, 59,60,61]. To identify the transcriptional motif signatures underlying differentially enriched H3K27ac regions in AZD6244-induced MEK inhibition, the HOMER motif analysis algorithm [62] was applied to the DNA sequences at differentially enriched H3K27ac regions. This analysis revealed enrichment of motif signatures of signal/stress-dependent regulatory proteins in the regions losing H3K27ac, where the top two enriched motifs were ETV and ATF binding motifs (Fig. 2e, Additional file 2: Table S3). Interestingly, numerous ETV and ATF transcription factor family members with similar consensus motif recognition showed significantly decreased mRNA expression under AZD6244 treatment, including ETV1, ETV4, ETV5, ATF3, and ATF4 (Fig. 1c, Additional file 2: Table S1), suggesting MEK inhibition may attenuate transcriptional programs regulated by these factors. In contrast, the highest ranked motifs in regions that gained H3K27ac corresponded to E-box and SOX consensus motifs (Fig. 2h, Additional file 2: Table S3). E-box and SOX motifs are recognized by the lineage-determining factors MITF and SOX10, respectively, which are essential for melanocyte development and differentiation, normal melanocyte function, and melanoma initiation and proliferation [45, 63,64,65,66,67]. In addition, SOX10, SOX4, SOX11, and MITF all showed significantly increased mRNA expression under AZD6244 treatment, whereas SOX5, which acts to inhibit SOX10 activation of downstream targets in melanocyte development [68], was significantly downregulated (Additional file 2: Table S1). Therefore, the expression changes in these transcription factors are consistent with the re-activation of developmental and pigmentation pathways.

Lineage-specific and pigmentation genes are associated with differentially acetylated H3K27ac regions

Given the enriched motifs for the melanocyte regulatory transcription factors MITF and SOX10, functional enrichment analyses using GREAT were performed to assess the corresponding genes associated with differentially enriched H3K27ac regions. The GREAT-identified genes for all 17,517 of the differentially enriched regions included 195 (78%) of the 250 pigmentation genes that show significant expression changes under MEK inhibition (Additional file 2: Table S2). The subset of genes associated with the 1041 regions gaining H3K27ac showed enrichment for pathways regulating morphogenesis and development, including oligodendrocyte differentiation, the phosphatidylinositol 3-kinase cascade, regulation of gliogenesis, and cell-fate specification (Additional file 1: Figure S2). These pathways are intriguing because melanocytes and glia of the peripheral nervous system (Schwann cells) are both derived from the neural crest lineage, and also because SOX10 is known to regulate lineage-specific/developmental pathways in melanocytes, glia, and oligodendrocytes [69,70,71,72,73]. In contrast, genes associated with the 16,476 regions losing H3K27ac showed significant enrichment of pathways regulating protein localization and protein targeting to membranes and endoplasmic reticulum (Additional file 1: Figure S2). Taken together, these data suggest the H3K27ac-marked regions may regulate lineage-specific, developmental, and pigmentation programs in response to MEK inhibition.

Identification of super-enhancers in 501mel cells and their association with pigment genes

Super-enhancers are characterized by extremely high densities of H3K27ac and transcription factors and regulate genes involved in development and disease [24, 30,31,32,33]. To define the constitutive super-enhancers found in 501mel cells, the ROSE algorithm [30, 31] was applied to the H3K27ac ChIP sequences from control (DMSO-treated) 501mel cells; this approach identified 13,523 typical enhancers and 799 super-enhancers (Fig. 3a, Additional file 2: Table S4). DNA motif enrichment analysis was performed on the 799 super-enhancers to identify putative interacting transcription factors, and this revealed significant enrichment of binding motifs for MITF, SOX4/SOX10, and ETS (Fig. 3b), transcription factors that are essential for melanocyte lineage development, differentiation, and melanomagenesis [63, 65,66,67].

Fig. 3
figure3

Identification of super-enhancers in 501mel cells and analysis of super-enhancer-associated gene response to MEK inhibition. a Distribution of H3K27ac ChIP-Seq signal at enhancers (input-subtracted reads per million per enhancer region) in 501mel cells under DMSO treatment. Enhancers are ranked by decreasing H3K27ac ChIP-Seq signal; the gray dashed line demarcates the boundary between super-enhancers (red) and typical enhancers (gray). The transcription factors listed exhibit 3 notable characteristics: each of their genomic loci are associated with nearby super-enhancers (indicated by black lines), each show differential mRNA expression under AZD6244 treatment (shown in Additional file 2: Table S1), and each also have enriched motifs when all super-enhancer regions are examined (shown in b). Therefore, these transcription factors are not only associated with a neighboring super-enhancer, but themselves have the capacity to bind other super-enhancers throughout the genome, raising the possibility of autoregulatory loops involving these key factors. Of note, MITF is shown twice because the MITF genomic locus has two independent super-enhancers. b The top five transcription factor binding motifs significantly enriched in 501mel super-enhancers; the fraction of super-enhancer regions containing the indicated motifs and the p values for the overrepresentation of each motif are shown. c Venn diagram showing the overlap between all genes associated with super-enhancers and genes whose expression is significantly altered by AZD6244 in 501mel cells. GO analysis was performed using the PANTHER overrepresentation test, to identify top biological pathways enriched for upregulated (upper panel, red) and downregulated (lower panel, green) genes with super-enhancer signatures in 501mel cells. p values represent FDR-adjusted (Bonferroni method) significance of enrichment for pathway members. d Venn diagram showing the overlap of super-enhancers in control (DMSO-treated) and AZD6244-treated 501mel cells

Full size image

The ROSE program was used to link the 799 super-enhancers to their putative target genes, identifying 1080 super-enhancer-associated genes. GO analysis of these super-enhancer-associated genes revealed top enriched biological processes regulating melanocyte differentiation and pigmentation, as well as neural crest development, adhesion, and migration (Additional file 2: Table S5). Comparison of the 1080 super-enhancer-associated genes with the 5526 differentially expressed genes (Additional file 2: Table S1) identified a set of 463 genes that were both super-enhancer-associated and exhibited altered expression under MEK inhibition (Fig. 3c). Interestingly, these 463 genes included transcription factors whose binding motifs are enriched in the 799 super-enhancers, including MITF, SOX10, and ETS1 (Fig. 3a). This suggests that these transcription factors have the potential to act at their associated super-enhancers in autoregulatory feedback loops. Furthermore, the upregulated subset of the 463 genes was enriched in biological pathways linked to neural crest/melanocyte development and pigmentation, in contrast to the broad regulatory pathways affecting cellular component movement and localization associated with the downregulated gene subset (GREAT analysis, Fig. 3c). Correspondingly, 45 of these 463 differentially expressed, super-enhancer-associated genes are known pigmentation genes, and 71% (32 of 45) of these pigment genes were upregulated under AZD6244 treatment (Additional file 2: Table S6). Taken together, these data support a role for these super-enhancers in regulating pigmentation and melanocyte/neural crest-related developmental pathways in 501mel cells.

Super-enhancers in AZD6244-treated 501mel cells lose H3K27ac but retain similar genomic locations

Using the same approach employed for identifying super-enhancers in control 501mel cells, we identified 583 super-enhancers and 721 super-enhancer-associated genes in AZD6244-treated 501mel cells (Additional file 1: Figure S3, Additional file 2: Table S4). Along with a reduced total number of super-enhancers (583 vs. 799), H3K27ac signal density at super-enhancers was significantly decreased under AZD6244 treatment (Additional file 1: Figure S3), which was reflected in an overall decrease in enhancer scores and was consistent with MEK inhibitor-induced loss of H3K27ac throughout the genome (Fig. 1a). Interestingly, while many control super-enhancers were no longer present following AZD6244 treatment (42%; 342 out of 799), the majority of super-enhancers identified in AZD6244-treated cells overlapped with control super-enhancers (83%; 484 out of 583, Fig. 3d), indicating retention of most super-enhancers at similar genomic locations following MEK inhibition. Correspondingly, 85% of differentially expressed genes associated with AZD6244 super-enhancers were the same as the control 463 gene set (264 out of 311 genes), including persistence of the motif-enriched transcription factors SOX10, MITF, ETV5, and ETS1 (Additional file 1: Figure S3). Notably, this overlap also included super-enhancers linked to pigmentation genes, as 28 of the 45 super-enhancers associated with differentially expressed pigmentation genes persisted under AZD6244 treatment (Additional file 2: Table S6).

AZD6244 treatment stimulates widespread recruitment of SOX10 to chromatin

SOX10 plays critical roles in melanocyte development and function, including direct transcriptional regulation of a subset of pigmentation genes [74,75,76,77,78,79,80,81]. Having demonstrated that MEK inhibition induces pigmentation in 501mel cells, and that SOX binding motifs are overrepresented at enhancer and super-enhancer regions, we used ChIP-Seq analysis to examine the genomic binding of SOX10. DNA motif enrichment analysis of SOX10 ChIP-Seq peaks revealed significant enrichment of both dimer and monomer SOX consensus motifs in control and AZD6244-treated cells (Additional file 1: Figure S5). The regions directly flanking SOX10 ChIP-Seq peak summits (both in control and AZD6244-treated cells) also contained overrepresented sequences for binding motifs of HBP1 and LEF1, transcriptional regulators of the WNT/β-catenin pathway, as well as E-box binding motifs for the transcription factors TFE3, TFEB, and MITF (Additional file 1: Figure S5).

Remarkably, SOX10 showed a massive increase in binding sites in response to AZD6244 treatment: 5716 discrete peaks were occupied by SOX10 in AZD6244-treated cells, which is a nearly threefold increase in comparison with the 1908 SOX10 peaks identified in control cells (Fig. 4a). Most SOX10 peaks in control cells were retained in AZD6244-treated samples (1508, 79% of control peaks). Quantitative comparisons of SOX10 ChIP-Seq tag density in control and AZD6244-treated samples using DiffBind showed that 4602 peaks were differentially occupied by SOX10; over 70% of these are located more than 50 kb away from the TSSs of neighboring genes (Fig. 4b), as expected from the previously described distal enhancer binding profile of SOX10 [40].

Fig. 4
figure4

MEK inhibition induces extensive SOX10 chromatin binding in 501mel cells. a Venn diagram showing overlap of SOX10 ChIP-Seq peaks in DMSO- and AZD6244-treated 501mel cells. In both treatment groups, the number of SOX10 ChIP-Seq peaks represents the total that were replicated across two independent experiments. Significance cutoff for SOX10 ChIP-Seq peaks was p < 1e10−5. b Genomic distribution of differentially bound SOX10 peaks in 501mel cells (AZD6244 vs DMSO) with respect to the nearest TSS as revealed by GREAT (using the basal plus extension settings). Values on the Y-axis represent percent of peaks found in a given genomic bin relative to TSS, i.e., [0, 5 kb], [5 kb, 50 kb], [50 kb, 500 kb], [500 kb, infinity]. c Heatmap of SOX10 ChIP-Seq signals at genomic regions that lost H3K27ac (upper panel) or gained H3K27ac (lower panel) in response to AZD6244 treatment (in units of read per million of genomic region). Each row shows a ± 5 kb region centered on an H3K27ac region, and rows are ordered by SOX10 ChIP-Seq signal. d Line plots showing average SOX10 ChIP-Seq signal intensities (y-axis) across regions losing H3K27ac (upper panel) and gaining H3K27ac (lower panel) over a ± 5 kb genomic window centered on differentially remodeled H3K27ac sites (x-axis)

Full size image

Given the extensive genome-wide recruitment of SOX10, we examined the H3K27ac-marked enhancers that exhibited novel SOX10 binding in response to AZD6244 treatment, and found enriched SOX10 binding at H3K27ac-marked regions. A total of 1935 enhancers showed novel SOX10 binding sites in AZD6244-treated cells, a more than threefold increase over the 526 SOX10-bound enhancers that were only seen in control (DMSO-treated) cells. In addition, examination of SOX10 binding at the 17,517 differentially enriched H3K27ac regions showed increased SOX10 binding under AZD6244 treatment, with greater SOX10 binding at regions with reduced H3K27ac (Fig. 4c, d, Additional file 1: Figure S4). SOX10 binding at super-enhancers was also increased under MEK inhibition. In control cells, 36% of super-enhancers (285 out of 799) contained SOX10 ChIP binding sites, while in AZD6244-treated cells the super-enhancers containing SOX10 binding increased to 66% (385 of 583). A detailed description of the dynamic SOX10 binding seen at all H3K27ac-marked regions in control and AZD6244-treated cells is contained in Additional file 1: Figure S4. The enrichment of SOX10 binding at the H3K27ac-marked enhancers and super-enhancers under AZD6244 treatment suggests that SOX10 may function to regulate the pathways associated with these regions by GO analysis, which included melanocyte differentiation and pigmentation, neural crest development, adhesion, and migration.

SOX10 binds at enhancers near pigmentation genes

GREAT analysis on the 1935 enhancers with novel SOX10 binding showed top enriched biological processes of cell chemotaxis and pigmentation; correspondingly, many of the H3K27ac-marked regions with enriched SOX10 binding were associated with pigmentation genes that showed significant expression changes under AZD6244 treatment. Out of the 250 differentially expressed pigmentation genes associated with H3K27ac-enriched regions, 76 genes showed SOX10 binding at their associated H3K27ac-enriched regions under AZD6244 treatment (Additional file 2: Table S2), with 62 of these genes showing novel SOX10 binding sites. (The remaining 14 SOX10 binding sites were constant between DMSO- and AZD6244-treatment.) In contrast, only 21 genes showed novel SOX10 binding to their associated H3K27ac-enriched regions under control conditions. The H3K27ac-enriched regions with SOX10 binding included super-enhancers: 21 of the 29 super-enhancers associated with differentially expressed pigmentation genes under AZD6244 treatment were bound by SOX10 (72%), 18 of which exhibited novel SOX10 recruitment under AZD6244 (Additional file 2: Table S6). Furthermore, the majority of the pigmentation genes associated with these SOX10-bound super-enhancers showed upregulated gene expression under AZD6244 treatment (86%, Additional file 2: Table S6). This upregulated gene set included MITF, SNAI2, ZEB2, ETS1, and TFAP2A, all critical transcription factors of the melanocyte lineage (Additional file 1: Figure S6).

To assess SOX10’s functional role in AZD6244-induced pigmentation, we reduced SOX10 protein levels in 501mel cells by the stable transfection of two independent shSOX10 constructs, each in triplicate. Complete ablation of SOX10 in melanoma cells induces growth arrest and senescence [46], and consistent with this finding SOX10 protein was partially reduced by all 6 shSOX10-expressing cell lines (Additional file 1: Figure S7). When treated with AZD6244, all shSOX10-expressing cell lines became pigmented, similar to control cells stably transfected with a non-targeting shRNA (Additional file 1: Figure S7), indicating that the partial reduction in SOX10 levels that occurred in these cells did not prevent the ability of AZD6244 to induce pigmentation.

GWAS comparisons support functional relevance of active chromatin regions to pigmentation-associated genetic variants

The genetic variants for human diseases and traits identified from genome-wide association studies (GWAS) primarily occur in non-coding genomic regions [19, 82, 83]. We therefore wanted to determine which GWAS single nucleotide polymorphisms (SNPs) associated with melanocyte phenotypes resided within our collection of melanocyte lineage distal enhancer elements. The H3K27ac-marked regions in 501mel cells were compared to 814 SNPs that have been previously implicated in human pigmentation and/or melanoma through GWAS [84,85,86]. These comparisons identified 90 pigmentation/melanoma GWAS loci that overlap with regions of H3K27ac acetylation and/or SOX10 ChIP-Seq binding (Table 1, Additional file 2: Table S7), thus highlighting the variants located within these regions as high priority candidates for future studies of the function of these distal enhancers. The associated phenotypes included normal pigmentation of hair, skin, and eye as well as suntan or sunburn, nevus formation, and cutaneous melanoma (Table 1).

Table 1 Overlap of 90 pigmentation and melanoma GWAS SNPs with H3K27ac-marked regions and with SOX10 ChIP binding
Full size table

Discussion

MAPK/ERK signaling is crucial for regulation of proliferation, differentiation, and pigmentation of melanocytes, and abnormal MAPK/ERK signaling leads to uncontrolled cell proliferation and melanoma tumorigenesis [87, 88]. We used the BRAF(V600E) mutant 501mel human cell line as a model system to discover the chromatin alterations that underlie the re-activation of melanocyte-intrinsic pigmentation in the context of MEK inhibition. Our chromatin-based approach showed that MEK inhibition re-activates melanocyte lineage-specific transcriptional programs, consistent with a previous study suggesting re-activation of lineage-specific programs under MEK inhibition [89]. Our results also showed that MEK inhibition extensively alters the active chromatin landscape and induces widespread SOX10 binding to additional genomic targets.

Our analysis used H3K27ac, a marker of active regulatory DNA elements including enhancers, to reveal the epigenetic and transcriptional alterations underlying MEK inhibition. AZD6244 treatment resulted in a dramatic reduction in genomic H3K27ac levels, with over 20,000 regions completely lost, and reduced H3K27ac levels at 94% of the 17,517 enhancer regions with significant acetylation changes. We also examined super-enhancers, because of their predicted role in governing the expression of genes affecting cell development/identity as well as oncogenic pathways [30,31,32]. We identified 799 super-enhancers in 501mel cells, showed that a subset of these remain constant under MEK inhibition, and also found that over 50% of these are remodeled in a MAPK-dependent manner, again with widespread loss in H3K27ac acetylation levels. Intriguingly, the loss of H3K27ac did not result in widespread loss of gene expression, as might be expected since H3K27ac marks transcriptionally active chromatin [28, 29], but instead resulted in relatively equal numbers of genes with increased and decreased expression. This suggests 501mel cells show a complex transcriptional response to the extensive H3K27ac deacetylation induced by MEK inhibition that is not simply a widespread attenuation of gene expression. The enhancers and super-enhancers identified here and their associated changes under MEK inhibition provide a window into the tremendous genome-wide reorganization that occurs at these regulatory genome features.

Enriched motif analysis of these enhancer and super-enhancer regions suggested that transcription factors with the capacity to bind ETV/ETS, ATF, SOX, E-BOX, and TCF motifs may play a role in regulation of these regions. Interestingly, several factors that bind to these overrepresented motifs are expressed in 501mel cells and either upregulated (SOX10, SOX4, SOX11, and MITF) or downregulated (SOX5, ETV1, ETV4, ETV5, ATF3, and ATF4) in response to MAPK inhibition. These results correlate with previous studies showing that MITF, ETS1, and ETV1 are downstream effectors of MAPK signaling in normal melanocytes and in melanoma progression [90,91,92,93]. Interestingly, the genomic loci for SOX10, MITF, ETS1, ETV1, ETV4, and ETV5 are associated with super-enhancers, suggesting the possibility for autoregulation; future studies will be needed to determine the roles for each of these transcription factor candidates in regulating these enhancers and super-enhancers in melanocytes and melanoma. Importantly, our study also demonstrates that SOX10 directly binds to numerous H3K27ac-marked regions in 501mel cells under MEK inhibition, consistent with SOX10 acting as a downstream effector of the MAPK pathway.

SOX10 is a member of the SOX family of transcription factors, which all exhibit strong transcriptional regulation of developmental processes and cell-fate specification in a lineage-specific manner [94, 95]. Recent studies have provided evidence that SOX10 occupies a key regulatory node in the response to MAPK signaling in the nucleus: SOX10 emerged from a shRNA screen for genes that confer resistance to BRAF and MEK inhibitors via upregulation of EGFR and PDGFRB [47], and SOX10 is directly phosphorylated by ERK as part of a regulatory axis involving ERK1/2, SOX10, FOXD3, and ERBB3 that mediates adaptive resistance to RAF inhibitors in BRAF mutant melanoma cells [49]. Our analysis now provides additional data for SOX10 acting downstream of the MAPK signaling cascade: the epigenetic reconfiguration of the genome under MEK inhibition included both a dramatic increase in genome-wide SOX10 chromatin occupancy and extensive SOX10 binding at enhancer and super-enhancer regions with decreased H3K27ac levels. These SOX10-bound enhancer and super-enhancer regions are associated with genes that affect pigmentation as well as the development and differentiation of neural crest and melanocytes, consistent with SOX10 functioning at these loci in re-activation of pigmentation and developmental pathways under MAPK inhibition in 501mel cells. These results for SOX10 in melanocytes correlate with previous data showing other SOX family members can bind to super-enhancers and regulate lineage-specific genes in chondrocytes and hair follicle stem cells [37, 38]. Strikingly, the super-enhancer-associated genes bound by SOX10 included MITF, TFAP2A, ETS1, ZEB2, and SNAI2, each lineage-specific transcription factors known to play important roles in melanocyte development and pigmentation [41, 66, 67, 96,97,98,99,100]. Several of these factors exhibit co-binding and/or regulation of one another during neural crest and melanocyte development; for example, ETS1 and SOX10 interact to promote proper melanocyte development [40, 99], TFAP2A and MITF interact to regulate pigmentation [54, 100], and SOX10 activates MITF, while MITF has been suggested to have the capacity to both activate and repress SOX10 expression [74,75,76,77, 101, 102]. SOX10 binding at super-enhancers near these genes suggests SOX10 is integral to the complex regulatory cross talk exhibited by these lineage-determining transcription factors. Interestingly, the partial knockdown of SOX10 did not prevent the AZD6244-induced appearance of pigmentation in 501mel cells; this may reflect the ability of the remaining SOX10 protein to still undergo genomic recruitment or possible compensation by other cellular mechanisms.

SOX10 recognizes either monomeric or palindromic dimer motifs (as do the other SOXE family members SOX8 and SOX9) and causes substantial bending of the DNA helix via binding to the minor groove, a characteristic shared by all SOX proteins [40, 103,104,105,106]. Our consensus motif profile data from SOX10 ChIP showed enrichment of both monomer and dimer binding profiles, suggesting the potential for monomeric and dimeric SOX10 binding at H3K27ac-marked enhancers. Interestingly, SOX10 monomer and dimer binding causes different angles of DNA bending (70°–85° or 101°, respectively, [103]). Bending of enhancer DNA has been suggested to regulate enhancer structure and transcriptional regulation, and furthermore, specific bending angles may facilitate differential recruitment of transcriptional cofactors, as has been shown for SOX2 [95, 107]. Future functional studies will be needed to determine if monomeric and dimeric binding of SOX10 at these H3K27ac-marked enhancers confers distinct functional properties. SOX10 is known to act at enhancers to recruit the SWI/SNF chromatin remodeling complex protein BRG1 in melanocytes [18] and also to recruit the mediator complex (which is critical for recruitment of RNA Polymerase II) via binding to MED12 in oligodendrocytes and Schwann cells [108]. While further cofactors that interact with SOX10 at DNA have been suggested in glial cells [109, 110], the extent of SOX10 cofactor recruitment and interaction remains under-explored in melanocytes.

SOX proteins generally confer transcriptional regulation by co-binding with additional transcription factors [111], and candidate SOX10-interacting transcription factors in melanocytes are suggested by the significantly enriched consensus motifs present at the SOX10 ChIP peaks identified in this study. Under both DMSO and AZD6244 treatment, consensus motifs were found for HBP1 and LEF1, critical transcription factors for WNT/B-catenin signaling [112,113,114], as well as E-box motifs for the MITF/TFE family of basic helix-loop-helix (bHLH) transcription factors. Co-localization of consensus sequences and chromatin binding for MITF and SOX10 has been observed previously in melanocytes [40, 41]. In agreement with these findings, a comparison of our SOX10 ChIP binding with previously identified MITF binding sites in 501mel cells [115] found that 42.7% of our SOX10 peaks colocalize with MITF binding regions. This suggests the opportunity for cooperative binding and synergistic action of MITF and SOX10, similar to what has been observed at the DCT promoter [78, 79]. Additionally, novel SOX10-interacting candidates are suggested by the presence of consensus motifs for FOXG1, FOXJ2, and FOXJ3, factors whose roles in melanocyte development or function have not been analyzed to date.

GWAS studies have identified over 800 SNPs associated with human pigmentation variation and melanocyte-related phenotypes [84,85,86]. Most GWAS SNPs reside in non-coding regions [19, 82, 83] and are hypothesized to reside in cis-regulatory regions or in close linkage disequilibrium (LD) to causative SNPs located within cis-regulatory regions. Therefore, GWAS results implicate the region of LD surrounding each SNP in affecting the associated phenotype, and subsequent analysis is needed to directly ascertain causative alleles and/or haplotypes. The large number of SNPs associated with human pigmentation reflects variation in GWAS studies performed in diverse populations, the identification of two or more distinct SNPs within the same LD region by different studies, multiple enhancers conferring regulation at the same genomic locus, and the polygenic nature of pigmentation itself. Integration of our epigenetic roadmap with GWAS variation found that 90 of these pigmentation/melanoma GWAS SNPs are located within SOX10 binding sites and/or enhancers that undergo dynamic changes under MEK inhibition, highlighting SNPs contained within these regions as intriguing candidates for future analyses of enhancer regulatory regions affecting pigmentation. The variants (and their respective LD regions) identified in this study provide an epigenetic link between the MAPK pathway, variation in normal human pigmentation, and melanoma susceptibility, thus potentially expanding the applicability of these data from 501mel cells to a broad range of melanocyte phenotypes, with potential implications for disease progression in MAPK-related disorders and drug responses that alter MAPK signaling.

Conclusions

This study reveals the extensive genomic responses that result from MEK inhibition in 501mel melanoma cells, including widespread remodeling of the active chromatin landscape, increased recruitment of SOX10 binding to the genome, and re-activation of melanocyte pigmentation and developmental transcriptional programs. Beyond specific mechanistic insight into the MEK inhibitor response in BRAF(V600E) mutant melanoma cells, this study provides a comprehensive genomic resource for future studies to identify potential novel loci for pigmentation and drug resistance in melanocytes and melanoma.

Methods

Cell culture and Western blot

The human melanoma cell line 501mel (kindly provided by Dr. Yardena Samuels) was cultured in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS) and DMEM, respectively. The human melanoma cell line A375 (ATCC, Manassas, VA) and the mouse melanoma cell line B16F10 (ATCC) were cultured in DMEM medium supplemented with 10% FBS. All cell lines were maintained at 37 °C under 5% CO2/95% humidified air. For shSOX10 knockdown, stable transfections of either a non-targeting shRNA or two independent shSOX10 constructs were performed in triplicate; then, the cells were subjected to puromycin selection for 2 weeks. All shRNA constructs (non-targeting shRNA: catalog # RHS6848; shSOX10-2: catalog# RHS3979-201750192 and shSOX10-5: catalog # RHS3979-201750195) were purchased from Dharmacon (Lafayette, CO). All cell lines were treated with DMSO (vehicle control) or 200 nM AZD6244 (MEK inhibitor) for 48 h (Western blot analyses of phosphorylated ERK) or 72 h (all studies). Western blots were performed according to standard protocols.

ChIP-Seq and peak calling

SOX10 and H3K27ac ChIP-Seq assays were performed essentially as described previously [40, 116, 117]. All ChIP-Seq data were derived from two independent experimental replicates of DMSO-treated and AZD6244-treated 501mel cells. Total genomic DNA (input) derived from formaldehyde cross-linked samples was used as control for normalization during peak calling. Raw sequence reads that passed quality control were aligned to the human reference genome GRCh37/hg19 (available from the UCSC genome browser, http://genome.ucsc.edu/) using Burrows Wheeler Alignment (BWA; [118]. Peak calling on all ChIP-Seq data was performed using MACS v2.1 [119] by applying the default settings and significance thresholds (format = AUTO (fragment length detection from libraries), model fold = [10, 100], q value cutoff = 5.00e−02 or FDR < 0.05), and using input genomic DNA as background control.

Differential binding analysis of ChIP-Seq data

After identifying ChIP-Seq peaks with MACS2, significantly differentially bound regions across samples were detected by quantification of H3K27ac ChIP-Seq tag density using the Bioconductor software package DiffBind [57]; http://bioconductor.org/packages/release/bioc/html/DiffBind.html), which was implemented in R as follows. MACS detected peaks that overlapped in at least two libraries, i.e., only replicated peaks, were retained by setting minOverlaps = 2 in the dba.count function of DiffBind, and all binding sites identified as overlapping were merged. In order to identify peaks with significant differences in read densities (differential binding affinity), we first obtained “consensus peaks” as follows: (1) overlapping peaks between DMSO and AZD6244 were merged and the largest 5′ to 3′ footprint was retained; (2) regions with no overlaps were also retained and included in the consensus peaks to enable analysis of peaks that are only present in one sample. Next, the number of reads within each experimental condition overlapping the consensus peaks was counted, and the overlapping reads from the corresponding genomic background (input) were subtracted. Then, contrasts were made between the treatment conditions and significantly differentially enriched regions were determined by implementing DESeq2 within DiffBind, which performs a normalization procedure based on trimmed mean of M values (TMM) using the overall depth of sequencing in each library (ChIP read counts minus input read counts and full library size) and calculates enrichment p values. High-confidence differentially enriched regions were obtained by correcting p values for multiple testing using the Benjamini–Hochberg procedure, and only sites with an FDR < 0.1 were retained for downstream analysis.

Identification of super-enhancers

Super-enhancers were identified using the rank ordering of super-enhancers (ROSE) algorithm (http://younglab.wi.mit.edu/super_enhancer_code.html, [31]). Briefly, peaks were called on H3K27ace ChIP-Seq bam files using MACS2 as described above using the human reference genome GRCh37/hg19 (see peak-calling procedure). Peaks located within 2500 base pairs (bp) of TSSs were excluded, as a 5-kb TSS exclusion zone was applied. The remaining peaks were designated as putative enhancers. Enhancers located within 12,500 base pairs of each other were stitched together, scored, and ranked based on the H3K27ac ChIP-Seq signal (total normalized read count in the ChIP-Seq sample after input subtraction). Then, enhancer signal (H3K27ac density) was plotted against enhancer rank. Enhancer regions above the inflection point of the curve were designated as super-enhancers, and those below the inflection point were designated as typical enhancers. Super-enhancers were assigned to the nearest RefSeq gene or genes using the 3 default parameters of the ROSE algorithm: included were overlapping genes (transcript directly overlapping the super-enhancer), proximal genes (TSS within 50 kb of the center of the super-enhancer), and nearest gene.

Motif discovery

H3K27ac ChIP-Seq regions were scanned for both de novo and known motifs using the findMotifsGenome.pl functionality in HOMER (Hypergeometric Optimization of Motif EnRichment, v4.6, http://homer.ucsd.edu/homer/index.html), with the default settings and whole genome as background [62]. For de novo motif analysis of super-enhancer regions, a 1000-bp sequence surrounding the center of each H3K27ac site was searched for enriched motifs using HOMER. SOX10 ChIP-Seq regions, defined by the peak summit and including ± 100 bp, were evaluated for both de novo and known motifs using MEME-ChIP [120].

RNA extraction and RNA-Seq data analysis

501mel cells were harvested with TRIzol reagent (Invitrogen/Thermo Fisher Scientific, Inc.) followed by RNA extraction using the Direct-zol™ RNA Miniprep Kit (Zymo Research). RNA quality was assessed by Bioanalyzer (Agilent Technologies, Inc.) and samples that passed the quality control metric (RIN > 6.5) were used for RNA-Seq. Paired-end poly-A enriched RNA-Seq libraries were generated by the NIH Intramural Sequencing Center (NISC) using an Illumina TruSeq small RNA kit from ~ 200 ng total RNA. The RNA samples were sequenced to generate at least 54 million 101-base read pairs for each individual library. Following quality filtering, RNA-Seq reads were mapped to the hg19 genome with STAR [121] and differential expression analysis of genes was performed using DESeq2 ([122], https://bioconductor.org/packages/release/bioc/html/DESeq2.html). Further RNA-Seq read quality control and calculation of number of mapped reads were performed using RSeQC [123] and RSEM [124], respectively. The normalized expression and differentially expressed genes between DMSO treatment and AZD6244 treatment conditions were identified using DESeq2 with an FDR of 0.1. Only nonzero tag counts in either condition were considered for further analysis.

Pigment gene identification

The 250 pigmentation genes identified in this study have been described in a previously published list of 650 genes that show visible integument pigment phenotypes in human, mouse, or zebrafish [51]. The 250 pigmentation genes were identified by determining the overlap of the 638 human orthologs contained within this published list with our gene lists. Of note, the visible pigmentation phenotypes associated with the 650 genes include both increased and decreased pigmentation, developmental abnormalities affecting melanocyte specification and survival, and pigmentation in adults, such as graying.

Availability of data and materials

ChIP-Seq (H3K27ac and SOX10) and RNA-Seq datasets of 501mel samples supporting the conclusions of this article have been deposited to the NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra/) under the Bioproject Accession Number PRJNA515302.

References

  1. 1.

    Montagut C, Settleman J. Targeting the RAF-MEK-ERK pathway in cancer therapy. Cancer Lett. 2009;283(2):125–34.

  2. 2.

    Santarpia L, Lippman SM, El-Naggar AK. Targeting the MAPK-RAS-RAF signaling pathway in cancer therapy. Expert Opin Ther Targets. 2012;16(1):103–19.

  3. 3.

    Burotto M, Chiou VL, Lee J-M, Kohn EC. The MAPK pathway across different malignancies: a new perspective. Cancer. 2014;120(22):3446–56.

  4. 4.

    Imperial R, Toor OM, Hussain A, Subramanian J, Masood A. Comprehensive pancancer genomic analysis reveals (RTK)-RAS-RAF-MEK as a key dysregulated pathway in cancer: its clinical implications. Semin Cancer Biol. 2019;54:14–28.

  5. 5.

    Davies H, Bignell GR, Cox C, Stephens P, Edkins S, Clegg S, et al. Mutations of the BRAF gene in human cancer. Nature. 2002;417(6892):949–54.

  6. 6.

    Davies BR, Logie A, McKay JS, Martin P, Steele S, Jenkins R, et al. AZD6244 (ARRY-142886), a potent inhibitor of mitogen-activated protein kinase/extracellular signal-regulated kinase kinase 1/2 kinases: mechanism of action in vivo, pharmacokinetic/pharmacodynamic relationship, and potential for combination in preclinical models. Mol Cancer Ther. 2007;6(8):2209–19.

  7. 7.

    Dry JR, Pavey S, Pratilas CA, Harbron C, Runswick S, Hodgson D, et al. Transcriptional pathway signatures predict MEK addiction and response to selumetinib (AZD6244). Cancer Res. 2010;70(6):2264–73.

  8. 8.

    Ambrosini G, Pratilas CA, Qin L-X, Tadi M, Surriga O, Carvajal RD, et al. Identification of unique MEK-dependent genes in GNAQ mutant uveal melanoma involved in cell growth, tumor cell invasion, and MEK resistance. Clin Cancer Res. 2012;18(13):3552–61.

  9. 9.

    Flaherty KT, Infante JR, Daud A, Gonzalez R, Kefford RF, Sosman J, et al. Combined BRAF and MEK inhibition in melanoma with BRAF V600 mutations. N Engl J Med. 2012;367(18):1694–703.

  10. 10.

    Kim KB, Kefford R, Pavlick AC, Infante JR, Ribas A, Sosman JA, et al. Phase II study of the MEK1/MEK2 inhibitor Trametinib in patients with metastatic BRAF-mutant cutaneous melanoma previously treated with or without a BRAF inhibitor. J Clin Oncol. 2013;31(4):482–9.

  11. 11.

    Ascierto PA, Schadendorf D, Berking C, Agarwala SS, van Herpen CM, Queirolo P, et al. MEK162 for patients with advanced melanoma harbouring NRAS or Val600 BRAF mutations: a non-randomised, open-label phase 2 study. Lancet Oncol. 2013;14(3):249–56.

  12. 12.

    Sullivan RJ, Fisher DE. Understanding the biology of melanoma and therapeutic implications. Hematol Oncol Clin North Am. 2014;28(3):437–53.

  13. 13.

    Niezgoda A, Niezgoda P, Czajkowski R. Novel approaches to treatment of advanced melanoma: a review on targeted therapy and immunotherapy. BioMed Res Int. 2015;2015:851387.

  14. 14.

    Yaeger R, Corcoran RB. Targeting alterations in the RAF-MEK pathway. Cancer Discov. 2019;9(3):329–41.

  15. 15.

    Gopal YNV, Deng W, Woodman SE, Komurov K, Ram P, Smith PD, et al. Basal and treatment-induced activation of AKT mediates resistance to cell death by AZD6244 (ARRY-142886) in Braf-mutant human cutaneous melanoma cells. Cancer Res. 2010;70(21):8736–47.

  16. 16.

    Kakadia S, Yarlagadda N, Awad R, Kundranda M, Niu J, Naraev B, et al. Mechanisms of resistance to BRAF and MEK inhibitors and clinical update of US Food and Drug Administration-approved targeted therapy in advanced melanoma. Onco Targets Ther. 2018;11:7095–107.

  17. 17.

    Kozar I, Margue C, Rothengatter S, Haan C, Kreis S. Many ways to resistance: how melanoma cells evade targeted therapies. Biochim Biophys Acta Rev Cancer. 2019;1871(2):313–22.

  18. 18.

    Marathe HG, Watkins-Chow DE, Weider M, Hoffmann A, Mehta G, Trivedi A, et al. BRG1 interacts with SOX10 to establish the melanocyte lineage and to promote differentiation. Nucleic Acids Res. 2017;45(11):6442–58.

  19. 19.

    Sud A, Kinnersley B, Houlston RS. Genome-wide association studies of cancer: current insights and future perspectives. Nat Rev Cancer. 2017;17(11):692–704.

  20. 20.

    See YX, Wang BZ, Fullwood MJ. Chromatin interactions and regulatory elements in cancer: from bench to bedside. Trends Genet. 2019;35(2):145–58.

  21. 21.

    Strub T, Ghiraldini FG, Carcamo S, Li M, Wroblewska A, Singh R, et al. SIRT6 haploinsufficiency induces BRAFV600E melanoma cell resistance to MAPK inhibitors via IGF signalling. Nat Commun. 2018;9(1):3440.

  22. 22.

    Webster DE, Barajas B, Bussat RT, Yan KJ, Neela PH, Flockhart RJ, et al. Enhancer-targeted genome editing selectively blocks innate resistance to oncokinase inhibition. Genome Res. 2014;24(5):751–60.

  23. 23.

    Long HK, Prescott SL, Wysocka J. Ever-changing landscapes: transcriptional enhancers in development and evolution. Cell. 2016;167(5):1170–87.

  24. 24.

    Fontanals-Cirera B, Hasson D, Vardabasso C, Di Micco R, Agrawal P, Chowdhury A, et al. Harnessing BET inhibitor sensitivity reveals AMIGO2 as a melanoma survival gene. Mol Cell. 2017;68(4):731–9.

  25. 25.

    Gelato KA, Schöckel L, Klingbeil O, Rückert T, Lesche R, Toedling J, et al. Super-enhancers define a proliferative PGC-1α-expressing melanoma subgroup sensitive to BET inhibition. Oncogene. 2018;37(4):512–21.

  26. 26.

    Moran B, Silva R, Perry AS, Gallagher WM. Epigenetics of malignant melanoma. Semin Cancer Biol. 2018;51:80–8.

  27. 27.

    Wang X, Hayes JJ. Acetylation mimics within individual core histone tail domains indicate distinct roles in regulating the stability of higher-order chromatin structure. Mol Cell Biol. 2008;28(1):227–36.

  28. 28.

    Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci. 2010;107(50):21931–6.

  29. 29.

    Rada-Iglesias A, Bajpai R, Swigut T, Brugmann SA, Flynn RA, Wysocka J. A unique chromatin signature uncovers early developmental enhancers in humans. Nature. 2011;470(7333):279–83.

  30. 30.

    Hnisz D, Abraham BJ, Lee TI, Lau A, Saint-André V, Sigova AA, et al. Super-enhancers in the control of cell identity and disease. Cell. 2013;155(4):934–47.

  31. 31.

    Whyte WA, Orlando DA, Hnisz D, Abraham BJ, Lin CY, Kagey MH, et al. Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell. 2013;153(2):307–19.

  32. 32.

    Lovén J, Hoke HA, Lin CY, Lau A, Orlando DA, Vakoc CR, et al. Selective inhibition of tumor oncogenes by disruption of super-enhancers. Cell. 2013;153(2):320–34.

  33. 33.

    Cao F, Fang Y, Tan HK, Goh Y, Choy JYH, Koh BTH, et al. Super-enhancers and broad H3K4me3 domains form complex gene regulatory circuits involving chromatin interactions. Sci Rep. 2017;7(1):320.

  34. 34.

    Parker SCJ, Stitzel ML, Taylor DL, Orozco JM, Erdos MR, Akiyama JA, et al. Chromatin stretch enhancer states drive cell-specific gene regulation and harbor human disease risk variants. Proc Natl Acad Sci. 2013;110(44):17921–6.

  35. 35.

    Eliades P, Abraham BJ, Ji Z, Miller DM, Christensen CL, Kwiatkowski N, et al. High MITF expression is associated with super-enhancers and suppressed by CDK7 inhibition in melanoma. J Invest Dermatol. 2018;138(7):1582–90.

  36. 36.

    Yohe ME, Gryder BE, Shern JF, Song YK, Chou H-C, Sindiri S, et al. MEK inhibition induces MYOG and remodels super-enhancers in RAS-driven rhabdomyosarcoma. Sci Transl Med. 2018;10(448):4470.

  37. 37.

    Adam RC, Yang H, Rockowitz S, Larsen SB, Nikolova M, Oristian DS, et al. Pioneer factors govern super-enhancer dynamics in stem cell plasticity and lineage choice. Nature. 2015;521(7552):366–70.

  38. 38.

    Liu C-F, Lefebvre V. The transcription factors SOX9 and SOX5/SOX6 cooperate genome-wide through super-enhancers to drive chondrogenesis. Nucleic Acids Res. 2015;43(17):8183–203.

  39. 39.

    Thakurela S, Sahu SK, Garding A, Tiwari VK. Dynamics and function of distal regulatory elements during neurogenesis and neuroplasticity. Genome Res. 2015;25(9):1309–24.

  40. 40.

    Fufa TD, Harris ML, Watkins-Chow DE, Levy D, Gorkin DU, Gildea DE, et al. Genomic analysis reveals distinct mechanisms and functional classes of SOX10-regulated genes in melanocytes. Hum Mol Genet. 2015;24(19):5433–50.

  41. 41.

    Laurette P, Strub T, Koludrovic D, Keime C, Le Gras S, Seberg H, et al. Transcription factor MITF and remodeller BRG1 define chromatin organisation at regulatory elements in melanoma cells. Elife. 2015;4:e06857.

  42. 42.

    Verfaillie A, Imrichova H, Atak ZK, Dewaele M, Rambow F, Hulselmans G, et al. Decoding the regulatory landscape of melanoma reveals TEADS as regulators of the invasive cell state. Nat Commun. 2015;6:6683.

  43. 43.

    Kaufman CK, Mosimann C, Fan ZP, Yang S, Thomas AJ, Ablain J, et al. A zebrafish melanoma model reveals emergence of neural crest identity during melanoma initiation. Science. 2016;351(6272):aad2197.

  44. 44.

    Hejna M, Moon WM, Cheng J, Kawakami A, Fisher DE, Song JS. Local genomic features predict the distinct and overlapping binding patterns of the bHLH-Zip family oncoproteins MITF and MYC-MAX. Pigment Cell Melanoma Res. 2019;32(4):500–9.

  45. 45.

    Shakhova O, Zingg D, Schaefer SM, Hari L, Civenni G, Blunschi J, et al. Sox10 promotes the formation and maintenance of giant congenital naevi and melanoma. Nat Cell Biol. 2012;14(8):882–90.

  46. 46.

    Cronin JC, Watkins-Chow DE, Incao A, Hasskamp JH, Schönewolf N, Aoude LG, et al. SOX10 ablation arrests cell cycle, induces senescence, and suppresses melanomagenesis. Cancer Res. 2013;73(18):5709–18.

  47. 47.

    Sun C, Wang L, Huang S, Heynen GJJE, Prahallad A, Robert C, et al. Reversible and adaptive resistance to BRAF(V600E) inhibition in melanoma. Nature. 2014;508(7494):118–22.

  48. 48.

    Shaffer SM, Dunagin MC, Torborg SR, Torre EA, Emert B, Krepler C, et al. Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature. 2017;546(7658):431–5.

  49. 49.

    Han S, Ren Y, He W, Liu H, Zhi Z, Zhu X, et al. ERK-mediated phosphorylation regulates SOX10 sumoylation and targets expression in mutant BRAF melanoma. Nat Commun. 2018;9(1):28.

  50. 50.

    The Gene Ontology Consortium. The gene ontology resource: 20 years and still GOing strong. Nucleic Acids Res. 2019;47(D1):D330–8.

  51. 51.

    Baxter LL, Watkins-Chow DE, Pavan WJ, Loftus SK. A curated gene list for expanding the horizons of pigmentation biology. Pigment Cell Melanoma Res. 2018;23(2):171.

  52. 52.

    Mártinez-García M, Montoliu L. Albinism in Europe. J Dermatol. 2013;40(5):319–24.

  53. 53.

    Zhang M, Song F, Liang L, Nan H, Zhang J, Liu H, et al. Genome-wide association studies identify several new loci associated with pigmentation traits and skin cancer risk in European Americans. Hum Mol Genet. 2013;22(14):2948–59.

  54. 54.

    Praetorius C, Grill C, Stacey SN, Metcalf AM, Gorkin DU, Robinson KC, et al. A polymorphism in IRF4 affects human pigmentation through a tyrosinase-dependent MITF/TFAP2A pathway. Cell. 2013;155(5):1022–33.

  55. 55.

    Crawford NG, Kelly DE, Hansen MEB, Beltrame MH, Fan S, Bowman SL, et al. Loci associated with skin pigmentation identified in African populations. Science. 2017;358(6365):eaan8433.

  56. 56.

    Mahoney SJ, Dempsey JM, Blenis J. Cell signaling in protein synthesis ribosome biogenesis and translation initiation and elongation. Prog Mol Biol Transl Sci. 2009;90:53–107.

  57. 57.

    Ross-Innes CS, Stark R, Teschendorff AE, Holmes KA, Ali HR, Dunning MJ, et al. Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature. 2012;481(7381):389–93.

  58. 58.

    McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28(5):495–501.

  59. 59.

    Arvey A, Agius P, Noble WS, Leslie C. Sequence and chromatin determinants of cell-type-specific transcription factor binding. Genome Res. 2012;22(9):1723–34.

  60. 60.

    Pham T-H, Benner C, Lichtinger M, Schwarzfischer L, Hu Y, Andreesen R, et al. Dynamic epigenetic enhancer signatures reveal key transcription factors associated with monocytic differentiation states. Blood. 2012;119(24):e161–71.

  61. 61.

    Liu Y, Pelham-Webb B, Di Giammartino DC, Li J, Kim D, Kita K, et al. Widespread mitotic bookmarking by histone marks and transcription factors in pluripotent stem cells. Cell Rep. 2017;19(7):1283–93.

  62. 62.

    Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38(4):576–89.

  63. 63.

    Harris ML, Baxter LL, Loftus SK, Pavan WJ. Sox proteins in melanocyte development and melanoma. Pig Cell Mel Res. 2010;23(4):496–513.

  64. 64.

    Vachtenheim J, Ondrušová L, Borovanský J. SWI/SNF chromatin remodeling complex is critical for the expression of microphthalmia-associated transcription factor in melanoma cells. Biochem Biophys Res Commun. 2010;392(3):454–9.

  65. 65.

    Shakhova O, Cheng P, Mishra PJ, Zingg D, Schaefer SM, Debbache J, et al. Antagonistic cross-regulation between Sox9 and Sox10 controls an anti-tumorigenic program in melanoma. PLoS Genet. 2015;11(1):e1004877.

  66. 66.

    Kawakami A, Fisher DE. The master role of microphthalmia-associated transcription factor in melanocyte and melanoma biology. Lab Invest. 2017;97(6):649–56.

  67. 67.

    Seberg HE, Van Otterloo E, Cornell RA. Beyond MITF: multiple transcription factors directly regulate the cellular phenotype in melanocytes and melanoma. Pigment Cell Melonoma Res. 2017;30(5):454–66.

  68. 68.

    Stolt CC, Lommes P, Hillgärtner S, Wegner M. The transcription factor Sox5 modulates Sox10 function during melanocyte development. Nucleic Acids Res. 2008;36(17):5427–40.

  69. 69.

    Britsch S. The transcription factor Sox10 is a key regulator of peripheral glial development. Genes Dev. 2001;15(1):66–78.

  70. 70.

    Stolt CC, Rehberg S, Ader M, Lommes P, Riethmacher D, Schachner M, et al. Terminal differentiation of myelin-forming oligodendrocytes depends on the transcription factor Sox10. Genes Dev. 2002;16(2):165–70.

  71. 71.

    Mollaaghababa R, Pavan WJ. The importance of having your SOX on: role of SOX10 in the development of neural crest-derived melanocytes and glia. Oncogene. 2003;22(20):3024–34.

  72. 72.

    Bremer M, Fröb F, Kichko T, Reeh P, Tamm ER, Suter U, et al. Sox10 is required for Schwann-cell homeostasis and myelin maintenance in the adult peripheral nerve. Glia. 2011;59(7):1022–32.

  73. 73.

    Hornig J, Fröb F, Vogl MR, Hermans-Borgmeyer I, Tamm ER, Wegner M. The transcription factors sox10 and myrf define an essential regulatory network module in differentiating oligodendrocytes. PLoS Genet. 2013;9(10):e1003907.

  74. 74.

    Bondurand N, Pingault V, Goerich DE, Lemort N, Sock E, Le Caignec C, et al. Interaction among SOX10, PAX3 and MITF, three genes altered in Waardenburg syndrome. Hum Mol Genet. 2000;9(13):1907–17.

  75. 75.

    Lee M, Goodall J, Verastegui C, Ballotti R, Goding CR. Direct regulation of the Microphthalmia promoter by Sox10 links Waardenburg-Shah syndrome (WS4)-associated hypopigmentation and deafness to WS2. J Biol Chem. 2000;275(48):37978–83.

  76. 76.

    Potterf SB, Furumura M, Dunn KJ, Arnheiter H, Pavan WJ. Transcription factor hierarchy in Waardenburg syndrome: regulation of MITF expression by SOX10 and PAX3. Hum Genet. 2000;107(1):1–6.

  77. 77.

    Verastegui C, Bille K, Ortonne JP, Ballotti R. Regulation of the microphthalmia-associated transcription factor gene by the Waardenburg syndrome type 4 gene, SOX10. J Biol Chem. 2000;275(40):30757–60.

  78. 78.

    Ludwig A, Rehberg S, Wegner M. Melanocyte-specific expression of dopachrome tautomerase is dependent on synergistic gene activation by the Sox10 and Mitf transcription factors. FEBS Lett. 2004;556(1–3):236–44.

  79. 79.

    Jiao Z, Mollaaghababa R, Pavan WJ, Antonellis A, Green ED, Hornyak TJ. Direct interaction of Sox10 with the promoter of murine Dopachrome Tautomerase (Dct) and synergistic activation of Dct expression with Mitf. Pigment Cell Res. 2004;17(4):352–62.

  80. 80.

    Murisier F, Guichard S, Beermann F. A conserved transcriptional enhancer that specifies Tyrp1 expression to melanocytes. Dev Biol. 2006;298(2):644–55.

  81. 81.

    Murisier F, Guichard S, Beermann F. The tyrosinase enhancer is activated by Sox10 and Mitf in mouse melanocytes. Pigment Cell Res. 2007;20(3):173–84.

  82. 82.

    Maurano MT, Humbert R, Rynes E, Thurman RE, Haugen E, Wang H, et al. Systematic localization of common disease-associated variation in regulatory DNA. Science. 2012;337(6099):1190–5.

  83. 83.

    Weinhold N, Jacobsen A, Schultz N, Sander C, Lee W. Genome-wide analysis of noncoding regulatory mutations in cancer. Nat Genet. 2014;46(11):1160–5.

  84. 84.

    MacArthur J, Bowler E, Cerezo M, Gil L, Hall P, Hastings E, et al. The new NHGRI-EBI Catalog of published genome-wide association studies (GWAS Catalog). Nucleic Acids Res. 2017;45(D1):D896–901.

  85. 85.

    Hysi PG, Valdes AM, Liu F, Furlotte NA, Evans DM, Bataille V, et al. Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability. Nat Genet. 2018;50(5):652–6.

  86. 86.

    Morgan MD, Pairo-Castineira E, Rawlik K, Canela-Xandri O, Rees J, Sims D, et al. Genome-wide study of hair colour in UK Biobank explains most of the SNP heritability. Nat Commun. 2018;9(1):5271.

  87. 87.

    Herraiz C, García-Borrón JC, Jiménez-Cervantes C, Olivares C. MC1R signaling. Intracellular partners and pathophysiological implications. Biochim Biophys Acta. 2017;1863(10 Pt A):2448–61.

  88. 88.

    Wellbrock C, Arozarena I. The complexity of the ERK/MAP-kinase pathway and the treatment of melanoma skin cancer. Front Cell Dev Biol. 2016;4:33.

  89. 89.

    Johannessen CM, Johnson LA, Piccioni F, Townes A, Frederick DT, Donahue MK, et al. A melanocyte lineage program confers resistance to MAP kinase pathway inhibition. Nature. 2013;504(7478):138–42.

  90. 90.

    Jané-Valbuena J, Widlund HR, Perner S, Johnson LA, Dibner AC, Lin WM, et al. An oncogenic role for ETV1 in melanoma. Cancer Res. 2010;70(5):2075–84.

  91. 91.

    Plotnik JP, Budka JA, Ferris MW, Hollenhorst PC. ETS1 is a genome-wide effector of RAS/ERK signaling in epithelial cells. Nucleic Acids Res. 2014;42(19):11928–40.

  92. 92.

    Wellbrock C, Arozarena I. Microphthalmia-associated transcription factor in melanoma development and MAP-kinase pathway targeted therapy. Pigment Cell Melanoma Res. 2015;28(4):390–406.

  93. 93.

    Xie Y, Cao Z, Wong EW, Guan Y, Ma W, Zhang JQ, et al. COP1/DET1/ETS axis regulates ERK transcriptome and sensitivity to MAPK inhibitors. J Clin Invest. 2018;128(4):1442–57.

  94. 94.

    Kamachi Y, Kondoh H. Sox proteins: regulators of cell fate specification and differentiation. Development. 2013;140(20):4129–44.

  95. 95.

    Julian LM, McDonald AC, Stanford WL. Direct reprogramming with SOX factors: masters of cell fate. Curr Opin Genet Dev. 2017;46:24–36.

  96. 96.

    Sánchez-Martín M, Rodríguez-García A, Pérez-Losada J, Sagrera A, Read AP, Sánchez-García I. SLUG (SNAI2) deletions in patients with Waardenburg disease. Hum Mol Genet. 2002;11(25):3231–6.

  97. 97.

    Sánchez-Martín M, Pérez-Losada J, Rodríguez-García A, González-Sánchez B, Korf BR, Kuster W, et al. Deletion of the SLUG (SNAI2) gene results in human piebaldism. Am J Med Genet A. 2003;122A(2):125–32.

  98. 98.

    Denecker G, Vandamme N, Akay O, Koludrovic D, Taminau J, Lemeire K, et al. Identification of a ZEB2-MITF-ZEB1 transcriptional network that controls melanogenesis and melanoma progression. Cell Death Differ. 2014;21(8):1250–61.

  99. 99.

    Saldana-Caboverde A, Perera EM, Watkins-Chow DE, Hansen NF, Vemulapalli M, Mullikin JC, et al. The transcription factors Ets1 and Sox10 interact during murine melanocyte development. Dev Biol. 2015;407(2):300–12.

  100. 100.

    Seberg HE, Van Otterloo E, Loftus SK, Liu H, Bonde G, Sompallae R, et al. TFAP2 paralogs regulate melanocyte differentiation in parallel with MITF. PLoS Genet. 2017;13(3):e1006636.

  101. 101.

    Elworthy S, Lister JA, Carney TJ, Raible DW, Kelsh RN. Transcriptional regulation of MITFA accounts for the sox10 requirement in zebrafish melanophore development. Development. 2003;130(12):2809–18.

  102. 102.

    Greenhill ER, Rocco A, Vibert L, Nikaido M, Kelsh RN. An iterative genetic and dynamical modelling approach identifies novel features of the gene regulatory network underlying melanocyte development. PLoS Genet. 2011;7(9):e1002265.

  103. 103.

    Peirano RI, Wegner M. The glial transcription factor Sox10 binds to DNA both as monomer and dimer with different functional consequences. Nucleic Acids Res. 2000;28(16):3047–55.

  104. 104.

    Antonellis A, Huynh JL, Lee-Lin S-Q, Vinton RM, Renaud G, Loftus SK, et al. Identification of neural crest and glial enhancers at the mouse Sox10 locus through transgenesis in zebrafish. PLoS Genet. 2008;4(9):e1000174.

  105. 105.

    Srinivasan R, Sun G, Keles S, Jones EA, Jang S-W, Krueger C, et al. Genome-wide analysis of EGR2/SOX10 binding in myelinating peripheral nerve. Nucleic Acids Res. 2012;40(14):6449–60.

  106. 106.

    Huang Y-H, Jankowski A, Cheah KSE, Prabhakar S, Jauch R. SOXE transcription factors form selective dimers on non-compact DNA motifs through multifaceted interactions between dimerization and high-mobility group domains. Sci Rep. 2015;5:10398.

  107. 107.

    Scaffidi P, Bianchi ME. Spatially precise DNA bending is an essential activity of the sox2 transcription factor. J Biol Chem. 2001;276(50):47296–302.

  108. 108.

    Vogl MR, Reiprich S, Küspert M, Kosian T, Schrewe H, Nave K-A, et al. Sox10 cooperates with the mediator subunit 12 during terminal differentiation of Myelinating Glia. J Neurosci. 2013;33(15):6679–90.

  109. 109.

    Weider M, Reiprich S, Wegner M. Sox appeal-Sox10 attracts epigenetic and transcriptional regulators in myelinating glia. Biol Chem. 2013;394(12):1583–93.

  110. 110.

    Lopez-Anido C, Sun G, Koenning M, Srinivasan R, Hung HA, Emery B, et al. Differential Sox10 genomic occupancy in myelinating glia. Glia. 2015;63(11):1897–914.

  111. 111.

    Kondoh H, Kamachi Y. SOX-partner code for cell specification: regulatory target selection and underlying molecular mechanisms. Int J Biochem Cell Biol. 2010;42(3):391–9.

  112. 112.

    Larue L, Delmas V. The WNT/Beta-catenin pathway in melanoma. Front Biosci. 2006;1(11):733–42.

  113. 113.

    Santiago L, Daniels G, Wang D, Deng F-M, Lee P. Wnt signaling pathway protein LEF1 in cancer, as a biomarker for prognosis and a target for treatment. Am J Cancer Res. 2017;7(6):1389–406.

  114. 114.

    Bollaert E, de Rocca Serra A, Demoulin J-B. The HMG box transcription factor HBP1: a cell cycle inhibitor at the crossroads of cancer signaling pathways. Cell Mol Life Sci. 2019;76(8):1529–39.

  115. 115.

    Strub T, Giuliano S, Ye T, Bonet C, Keime C, Kobi D, et al. Essential role of microphthalmia transcription factor for DNA replication, mitosis and genomic stability in melanoma. Oncogene. 2011;30(20):2319–32.

  116. 116.

    Lee TI, Johnstone SE, Young RA. Chromatin immunoprecipitation and microarray-based analysis of protein location. Nat Protoc. 2006;1(2):729–48.

  117. 117.

    Gorkin DU, Lee D, Reed X, Fletez-Brant C, Bessling SL, Loftus SK, et al. Integration of ChIP-seq and machine learning reveals enhancers and a predictive regulatory sequence vocabulary in melanocytes. Genome Res. 2012;22(11):2290–301.

  118. 118.

    Li H, Durbin R. Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics. 2010;26(5):589–95.

  119. 119.

    Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137.

  120. 120.

    Ma W, Noble WS, Bailey TL. Motif-based analysis of large nucleotide data sets using MEME-ChIP. Nat Protoc. 2014;9(6):1428–50.

  121. 121.

    Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.

  122. 122.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

  123. 123.

    Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012;28(16):2184–5.

  124. 124.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12(1):323.

Download references

Acknowledgements

We thank all members of the Pavan laboratory for helpful discussions and careful reading of this manuscript.

Funding

This research was supported by the Intramural Research Program of the NIH, NHGRI. This study was funded by National Human Genome Research Institute (1ZIAHG000136-20).

Author information

Affiliations

  1. Genetic Disease Research Branch, National Human Genome Research Institute, National Institutes of Health, Bethesda, MD, 20892, USA

    • Temesgen D. Fufa
    • , Laura L. Baxter
    • , Julia C. Wedel
    • , Stacie K. Loftus
    •  & William J. Pavan

  2. Ophthalmic Genetics and Visual Function Branch, National Eye Institute, National Institutes of Health, Bethesda, MD, 20892, USA

    • Temesgen D. Fufa

  3. Computational and Statistical Genomics Branch, National Human Genome Research Institute, National Institutes of Health, Bethesda, MD, 20892, USA

    • Derek E. Gildea

  4. NIH Intramural Sequencing Center, National Human Genome Research Institute, National Institutes of Health, Bethesda, MD, 20892, USA

    Authors

    1. Search for Temesgen D. Fufa in:

    2. Search for Laura L. Baxter in:

    3. Search for Julia C. Wedel in:

    4. Search for Derek E. Gildea in:

    5. Search for Stacie K. Loftus in:

    6. Search for William J. Pavan in:

    Consortia

    NISC Comparative Sequencing Program

    Contributions

    TF and JW designed and performed experiments. TF, LB, DG, NISC Comparative Sequencing Program, and SL analyzed data. TF, LB, SL, and WP wrote and revised the manuscript. WP supervised and obtained funding for the project. All authors read and approved the final manuscript.

    Corresponding author

    Correspondence to William J. Pavan.

    Ethics declarations

    Ethics approval and consent to participate

    Not applicable.

    Consent for publication

    Not applicable.

    Competing interests

    The authors declare that they have no competing interests.

    Additional information

    Publisher's Note

    Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

    Additional files

    13072_2019_297_MOESM1_ESM.pptx

    Additional file 1: Figure S1. MEK inhibition induces hyperpigmentation in mouse B16F10 and human A375 melanoma cell lines. Cells were treated with 200 nM AZD6244 for 72 h. Cell suspensions are shown above cell pellets, with equal numbers of cells present in samples with and without AZD6244 treatment. Note that the B16F10 cells are pigmented initially, and then exhibited a further darkening of pigmentation under AZD6244 treatment. Figure S2. Enriched biological processes identified by GREAT analyses on the 16,476 regions that were losing H3K27ac and the 1041 regions that were gaining H3K27ac under AZD6244 treatment. Figure S3. Effect of MEK inhibition on super-enhancers in 501mel cells. A) Distribution of H3K27ac ChIP-Seq signal at enhancers (input-subtracted reads per million per enhancer regions) in AZD6244-treated 501mel cells. Enhancers are ranked by decreasing H3K27ac ChIP-Seq signal; the gray dashed line demarcates the boundary between super-enhancers (red) and typical enhancers (gray). Of the six super-enhancer-associated transcription factors with motifs enriched in 501mel super-enhancer regions (Fig. 3b), four (MITF, SOX10, ETS1 and ETV5) persisted as super-enhancer-associated under AZD6244 treatment. B) Box-and-scatter plots of super-enhancer signal density (rpm/bp) in control (SE_DMSO) and AZD6244 treated (SE_AZD) 501mel cells. Significance of the difference between distributions determined using a two-tailed t test. ***p value ≤ 0.001. Figure S4. SOX10 shows increased genome-wide binding under AZD6244 treatment, including recruitment to H3K27ac-marked regions. The numbers of SOX10 binding sites from SOX10 ChIP analysis throughout the genome and in relation to H3K27ac-marked regions in 501mel cells are shown, both under control (DMSO-treated) and AZD6244-treated conditions. Figure S5. Significantly enriched transcription factor motifs identified under SOX10 ChIP-Seq peaks. Enriched motifs are shown for (A) control (DMSO-treated) conditions and (B) AZD6244-treated conditions. Figure S6. Transcription factors that regulate melanocyte lineage and pigmentation processes undergo significant expression changes under AZD6244 and are associated with SOX10-bound super-enhancers. UCSC genome browser tracks are shown for the transcription factors MITF, ETS1, ZEB2, SNAI2, and TFAP2A. Colored tracks indicate control and AZD6244-treated regions with H3K27ac binding (pink and gray), with differentially acetylated H3K27ac levels (orange and brown), super-enhancers (light and dark green), and SOX10-ChIP genomic binding (teal and blue). Note the novel recruitment of SOX10 to the super-enhancer regions associated with these genes resulting from AZD6244 treatment (teal sites that do not overlap blue sites). Figure S7. Partial reduction of SOX10 protein using shRNA knockdown in 501mel cells does not prevent the ability of AZD6244 to induce pigmentation. We performed stable knockdowns in triplicate experiments using two different shSOX10 constructs (shSOX10 #2 and shSOX10 #5), and these cells were then subjected to either DMSO or AZD6244 treatment for 72 h. Western blotting for SOX10, phospho-ERK and GAPDH (loading control) showed that each shSOX10 induced a partial knockdown of SOX10 protein, and AZD6244 treatment reduced phosphorylated ERK levels. Visual inspection of the cells indicated that while all cells under DMSO treatment remained unpigmented (indicated by “−” sign), pigmentation was induced under AZD6244 treatment for all cells (indicated by “+” sign), including in the context of shSOX10 treatment. NT = non-targeting.

    13072_2019_297_MOESM2_ESM.xlsx

    Additional file 2: Table S1. Genes exhibiting significantly different expression changes in 501mel cells in response to AZD6244 treatment. Table S2. The 250 pigmentation genes with differential expression in 501mel cells under AZD6244 treatment. Genes associated with H3K27ac loss or gain regions by GREAT analysis, or with SOX10 binding sites under AZD6244 are also indicated. Table S3. Top 5 significantly enriched transcription factor motifs identified by de novo motif enrichment analysis of genomic regions that were losing or gaining H3K27ac in 501mel cells in response to AZD6244 treatment. Table S4. The 799 super-enhancers identified in 501mel cells treated with DMSO only, and the 583 super-enhancers identified in 501mel cells treated with AZD6244. Table S5. Top enriched biological processes from GO analyses of 1080 genes associated with super-enhancers in control (DMSO-treated) 501mel cells. Table S6. Super-enhancer-associated pigmentation genes with significant expression changes following AZD6244 treatment. Table S7. GWAS loci associated with pigmentation and/or melanoma phenotypes that overlap with H3K27ac regions and/or with SOX10 binding sites.

    Rights and permissions

    Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

    Reprints and Permissions

    About this article

    Verify currency and authenticity via CrossMark

    Keywords

    • Pigmentation
    • Melanocyte
    • MAPK pathway
    • MEK inhibitor
    • Enhancers
    • SOX10

    Epigenetics & Chromatin

    ISSN: 1756-8935