Zebrafish MITF-Low Melanoma Subtype Models Reveal Transcriptional Subclusters and MITF-Independent Residual Disease

Melanoma is the most lethal form of skin cancer, and its incidence continues to rise rapidly. Despite the significant progress in targeted therapy against BRAF mutations and MEK activity, and in immune therapy, most patients with metastatic melanoma still succumb to the disease because they do not respond to therapy or they develop drug resistance (1). Genetic mutations define subgroups of melanoma [BRAF, RAS or NF1 mutated, or triple-wild type (WT)] that can help predict a positive response to targeted therapy; however, overall this genomic classification cannot predict patient survival (2).

In contrast, gene expression profiling of melanomas has revealed transcription-based melanoma subtypes that have significance for disease outcome and may be important for patient care (3). The Cancer Genome Atlas (TCGA) melanoma dataset revealed three transcriptionally distinct subtypes that are enriched for immune genes, for keratinocyte genes, or have low expression of the melanocyte-inducing transcription factor (MITF) and its targets as well as high expression of neuronal genes (called the MITF-low signature; ref. 2).

The TCGA MITF-low signature in primary and metastatic melanomas shares over 75% overlap with an MITF-low proliferation transcriptional subtype identified in the Lund dataset of stage IV metastatic melanomas (4). The TCGA MITF-low and Lund MITF-low proliferation groups are further found in the primary (Leeds), stage III (Lund), and metastatic (Bergen) transcriptome melanoma datasets (5–7). Together, these data support the concept that MITF-low transcriptional activity defines a melanoma biological subtype (8). Importantly, outcomes for patients with melanomas that express the TCGA MITF-low signature or the Lund MITF-low proliferative signature genes are among the poorest for survival and metastasis (8). However, whether the MITF-low signature is simply a marker for poor outcomes or whether it is causative remains an open question.

Here, we use zebrafish genetics coupled with analysis of human melanoma transcriptomes to determine the biological significance of the MITF-low classification and reveal a preexisting melanoma cell subpopulation that is resistant to loss of MITF activity and has very low-to-no MITF activity toward its target genes (termed MITF-independent cells). These cells remain at the site of residual disease and could be an important drug target.

Zebrafish were maintained in accordance with UK Home Office regulations, UK Animals (Scientific Procedures) Act 1986, under project license 70/8000 and P8F7F7E52. All experiments were approved by the Home office and AWERB (University of Edinburgh Ethics Committee).

Zebrafish were genotyped using DNA extracted from fin clipped tissue by PCR to establish the mutant allele status tp53^M214K or mitfa-BRAF^V600E and grown in temperature-controlled systems as described (9). Briefly, for melanoma induction, or recurrence, fish were transferred to glass tanks at room temperature (24–26°C), whereas for melanoma regression, fish were transferred to tanks with heaters that kept the temperature at 32°C.

Fish samples were collected, fixed, and processed as described (9). The detailed method is described in Supplementary Materials and Methods. Antibodies and their dilutions are listed in Supplementary Table S1.

For low cell number and single-cell RNA sequencing (RNA-seq), the Smart-seq2 protocol was used with minor modifications as described previously (10) and detailed in Supplementary Materials and Methods.

Reads were aligned to the Ensembl Danio rerio reference genome version GRCz11 (Ensembl V92, 20178) that included coding sequences for EGFP and BRAF^V600E using STAR (V 2.5.1b; ref. 11). For quality control and preprocessing, quantification of mapped reads per gene was calculated using HTseq (12). Genes that were not expressed in any cells were excluded. Data were not normalized for any variable based on histologic evaluation. The raw read counts were loaded in the R DEseq2 package (V 1.20.0, ref. 13). For bulk samples, genes with fewer than 10 reads across the samples were excluded. Quality control exclusion criteria for low cell number RNA-seq samples were the following: samples with fewer than 100,000 unique-mapped reads and fewer than 1,000 genes detected (minimum one read per gene) were excluded. For differential gene expression analysis, genes displaying fewer than ten counts across all samples for bulk RNA-seq or that did not express with at least one count in a minimum of three samples for low number RNA-seq were excluded.

Single-cell RNA-seq analysis was performed using adapted pipeline from Lun and colleagues (14). Detailed description of the analysis is found in Supplementary Materials and Methods.

Genes upregulated from DESeq2 (bulk samples) or scDE (single cell) were selected based on FDR < 0.05 and ranked by Log₂FC (Log₂FC>0) for input to the R ClusterProfiler package (15). Pathway overrepresentation analyses were based on available literature datasets (abbreviated by L: Supplementary Table S2; all genes manually converted to zebrafish orthologs using zfin.org), or using Kyoto Encyclopedia of Genes and Genomes (KEGG) (K:) and GO-BP (G:) databases. Pathways with P.adj < 0.05 (hypergeometric test, BH adjusted) were considered to be significant, listed in Supplementary tables with subsets visualized using dot-plots.

Gene set functional enrichment analysis (GSEA; ref. 16) identified enriched pathways at FDR < 0.05 (Kolmogorov–Smirnov test). The matrix of all read counts was used to compare against literature-based datasets (Supplementary Table S2) using gene set permutation settings. FDR values and leading-edge genes of enrichment plots shown in figures are listed in Supplementary Table S3.

Detailed method is described in Supplementary Materials and Methods with primers listed in Supplementary Table S4.

Data availability: RNA-seq data are deposited in the GEO database (GSE130037 and GSE136900).

For additional Materials and Methods, please see Supplementary Information.

The MITF gene is highly conserved, and in zebrafish, there are two orthologous mitf genes (a and b) whereby mitfa is required for the development of body melanocytes (17). A genetic temperature-sensitive mitfa mutation (mitfa^vc7) in zebrafish has enabled us to test the impact of altering MITF activity in development and to show that low MITF activity cooperates with BRAF^V600E in melanomagenesis (9, 18, 19). The conditional mitfa^vc7 mutant allele results from a single-nucleotide T-to-A mutation in intron 6 that disrupts splicing in a temperature-sensitive manner (18, 20). Analysis of the splicing variants of mitfa^vc7 indicates that at permissive temperatures (25°C) lower levels of WT mitfa RNA are generated (defined as a MITF-low state), and that these transcripts are further reduced with increasing water temperatures (Fig. 1A; refs. 9, 20). Thus, this conditional control of mitfa RNA enables us to regulate WT mitfa expression in zebrafish development and cancer.

In the TCGA data, an MITF-low transcriptional signature is found in 18% of cutaneous melanomas and across all four TCGA cutaneous melanoma genomic groups (Fig. 1B). TP53 (p53) is frequently mutated in human melanoma (18% of cutaneous melanoma; TCGA), and these mutations cooperate with BRAF^V600E mutations to promote melanoma in zebrafish and mouse models (21, 22). We aimed to model the MITF-low transcriptional state using genetic mutation combinations present in MITF-low melanomas in patients. We explored the human melanoma genetic and transcriptomic TCGA dataset and found TP53 mutations or low levels of TP53 expression in BRAF-mutant and non–BRAF-mutant melanomas, including a minor subset that has the MITF-low transcriptional signature (Fig. 1B). To test the potential for p53 mutations to cooperate with MITF-low activity in nevi and melanoma, we generated zebrafish homozygous for mitfa^vc7 and tp53 (p53^M214K) mutations. The impact of the conditional mitfa^vc7 activity on zebrafish stripe formation can be seen in Fig. 1C (left panels), whereby zebrafish grown at <25°C (MITF-low activity state) develop fewer melanocytes in a stripe pattern, whereas animals grown at 32°C lack sufficient functional MITF activity to develop body melanocytes.

Strikingly, we found that the mitfa^vc7;p53^M214K animals developed large darkly pigmented melanocyte lesions by 9 months that progressed into melanoma with 73% incidence (n = 94/129; Fig. 1C and D). These data reveal an unexpected cooperation between MITF-low activity and p53 mutations in melanoma initiation that is independent of MAPK activation. To address the impact of the BRAF^V600E mutation on mitfa^vc7;p53^M214K zebrafish, we crossed the mitfa^vc7;p53^M214K zebrafish to zebrafish expressing human BRAF^V600E from the mitfa promoter (21). This promoter is not regulated by Mitfa protein levels, and thus BRAF^V600E expression is not altered by levels of Mitfa in the mitfa^vc7 mutant (9). BRAF^V600E dramatically accelerated the onset of melanoma and gave rise to melanomas with both superficial growth and exophytic nodular growth patterns, and could be either pigmented or unpigmented (Fig. 1C–E). Taken together, these results reveal that BRAF^V600E, MITF-low activity, and p53 mutations have an independent and cooperative impact on melanoma initiation and growth.

Next, we examined the pathology of the melanoma subtypes. Histologic examination showed all tumors to be highly cellular and express Melan A (Supplementary Fig. S1A and S1B; n = 15/15). Melanomas from the mitfa^vc7;p53^M214K fish originated from the skin, grew in a superficial growth pattern, and were highly invasive into the muscle, penetrating the tissue to deeper levels, around the internal organs and toward the spinal cord (Fig. 2A–C). We term these melanomas double superficial (DS) because they have two engineered genetic mutations and grow in a superficial growth pattern. The melanomas in the Tg(mitfa:BRAF^V600E); mitfa^vc7;p53^M214K fish grew in either a superficial or a highly proliferative nodular growth pattern [referred to as triple superficial (TS) or triple nodular (TN); Supplementary Fig. S1C and S1D]. Both the superficial and nodular growth melanomas were highly invasive, but they differed in the nature of invasion. We measured from the level of the background contour of the fish to the deepest invasive cell and found no difference between the superficial and nodular growth melanomas in invasion depth. In contrast, there was a marked difference in the mode of invasion: superficial growth melanomas showed a single-to-few cell-wide stream of highly infiltrative cancer cells between muscle fibers, whereas nodular growth melanomas showed a more confluent invasive growth pattern often leading to deterioration of surrounding tissues (Fig. 2D–F).

Notably, despite heavy pigmentation in the DS and TS tumors at the surface site, melanoma cells of the invasive region were unpigmented (Fig. 2A–C). DS melanomas showed little-to-no activation of the MAPK-signaling pathway in the superficial part of the tumor, but the unpigmented invasive compartments expressed high levels of phospho-ERK in a heterogeneous pattern (Fig. 2G and H). To relate this finding to human melanoma, we asked whether MAPK signaling pathways are associated with an increase in invasive properties in TCGA triple-WT melanomas. We found a significant positive correlation between the invasive signature genes and MAPK-signaling signature genes (Supplementary Fig. S1E), indicating that MAPK signaling is associated with melanoma invasion in both zebrafish and human triple-WT melanomas. This is consistent with the findings of Shain and colleagues (23) in which genetic and transcriptomic activation of the MAPK signaling pathway increases with disease progression.

To analyze the three melanoma models, we dissected the DS, TS, and TN tumors from melanoma-bearing fish and performed bulk RNA-seq (Fig. 3A; Supplementary Fig. S2A–S2C). Principal component analysis divided the melanoma transcriptomes into three groups (Fig. 3B). Along the PC1 axis, melanomas separated according to the direction of growth (superficial or nodular), and along the PC2 axis, melanomas separated according to the presence or absence of BRAF^V600E. Although pigmented nodular growth melanomas are found between the BRAF^V600E superficial growth and the BRAF^V600E unpigmented nodular growth melanomas in the PCA (Fig. 3B), they clustered together with BRAF^V600E unpigmented nodular growth melanomas by hierarchical clustering (Supplementary Fig. S3A).

We quantified mRNA expression of mitfa between DS and TS superficial growth melanomas and found no significant difference in its expression (Fig. 3C). Indeed, we found that most MITF target gene expression levels (24) were similar in both DS and TS melanomas, suggesting that the BRAF^V600E mutation alone does not significantly alter the overall MITF-dependent transcriptome of MITF-low tumors (Fig. 3C; Supplementary Fig. S3B). In contrast, there was a marked increase in mitfa expression in TN growth melanomas (Fig. 3C). This change was also accompanied by an increase in BRAF^V600E mRNA expression from the mitfa promoter in addition to an increase in expression of MITF target genes involved in proliferation and differentiation (cdk2, dct, tyr, mc1r, tyrp1a/b).

Given that the mitfa^vc7 mutant is a splice-site mutation, the increase in mitfa expression by RNA-seq might simply reflect an increase in aberrantly spliced nonfunctional transcripts. To test this, we designed primers that specifically amplified correctly spliced mitfa exons and found a significant increase in functional WT mitfa transcripts in the TN melanomas compared with DS and TS melanomas (Fig. 3D and E). This was associated with a modest increase in Mitfa protein in TS melanomas, and a robust increase in Mitfa protein in the TN melanomas (Fig. 3F and G). Based on these data, we classified DS melanomas as MITF-low melanomas and TN melanomas as MITF-high melanomas. Some TS melanomas had higher RNA and protein levels of Mitfa; however we also classified these as MITF-low melanomas because they shared a strong similarity of MITF target gene expression with DS melanomas and clustered together by unsupervised clustering (Fig. 3C; Supplementary Fig. S3A and S3B). Further, MITF has been proposed to function as a repressor at some gene promoters (25), and we find derepression of some of these genes in MITF-low melanomas, including genes involved in stem cell states (klf4, klf5a, klf9, erbb3) and metabolism (acsf2, apodb, mboat1; Supplementary Fig. S3B). These results suggest that DS and TS MITF-low melanomas are transcriptionally distinct from TN MITF-high melanomas in MITF activity.

In human melanoma, MITF gene expression is controlled by neural crest transcription factors including PAX3 and BRN2 (also known as POU3F2) that act in a reciprocal manner to control high and low MITF gene expression, respectively (26). Increased mitfa gene expression in BRAF^V600E nodular growth melanomas was accompanied by an increase in pax3a expression, whereas mitfa-low expression level in BRAF^V600E melanomas with superficial growth was associated with upregulation in brn2b (pou3f2b) expression. Both pax3a and brn2b (pou3f2b) were downregulated in (non-BRAF^V600E) MITF-low p53^mut melanomas (DS, Fig. 3C). This suggests that a Brn2b-Pax3a reciprocal mechanism might be active only in BRAF^V600E melanomas and be responsible for a switch from low to high mitfa RNA expression.

We sought to relate our findings in zebrafish melanoma to MITF-low melanomas from patients. We established the genes that were differentially expressed between superficial growth (low mitfa expression) and nodular growth (higher mitfa expression) zebrafish melanomas, and then compared these with the genes that were differentially expressed between 40 human TCGA melanomas with the lowest and 40 melanomas with the highest MITF mRNA expression (Fig. 3H). MITF-low superficially growing zebrafish melanomas (both with or without BRAF^V600E) shared marked and significant similarity to human melanomas with low MITF mRNA expression (Fig. 3I) such as mesenchymal and invasive programs, retinoic acid metabolic processes, as well as neural crest cell development, migration, and differentiation signatures that include genes involved in stem cell and mesenchymal fates (Fig. 3I; Supplementary Tables S2 and S5). These results show that zebrafish MITF-low melanomas model human MITF-low transcriptional signature melanomas at the molecular level.

Nodular growth melanomas showed less resemblance to the 40 human melanomas with the highest MITF mRNA expression levels; however, some pathways were significantly shared between nodular zebrafish melanomas and human melanomas with the highest MITF RNA expression level, including the Lund MITF-high pigmentation cluster (Supplementary Fig. S4A and S4B; Supplementary Table S6). These data indicate that zebrafish MITF-low superficial growth melanomas model the TCGA MITF-low melanomas, whereas zebrafish MITF-high nodular melanomas more closely align with the Lund MITF-high pigmentation melanomas at the transcriptional level.

Next, we performed a pairwise comparison of differential mRNA expression from each zebrafish melanoma subtype to identify the molecular signatures that define them (Fig. 4A; Supplementary Table S7). DS melanomas expressed high levels of WNT and mesenchymal genes (Fig. 4A–D), as well as the vascular genes pdgfra and kdr that are orthologous to human PDGFRA and KDR and are characteristic of human Triple-WT melanomas (Supplementary Fig. S4C and S4D). Surprisingly, the TS melanomas showed upregulation of neuronal genes, including the neuronal gene set specifically upregulated in MITF-low TCGA transcriptional clusters (Fig. 4A and D). TN melanomas that were enriched for MITF program activation expressed proliferation genes and a neural crest formation program (Fig. 4A–D). This MITF-high state includes expression of sox10, foxd3, pax3, and tfap2 genes, as well as crestin (a zebrafish neural crest marker gene) that is expressed in melanoma-initiating cells in zebrafish BRAF^V600E p53^mut melanomas (Fig. 4B–D; ref. 27).

Neural crest genes were expressed in all melanoma subtypes. Nevertheless, the levels and gene signatures differed between them. The superficial growth melanomas expressed neural crest signatures that are broadly stem-like (e.g., erbb2, erbb3a/b) and mesenchymal (e.g., sox9, twist1a/1b/3), whereas the nodular growth melanomas expressed genes involved in developmental programs (e.g., sox10, foxd3, mitfa; Supplementary Fig. S4E). We validated this finding using a zebrafish-specific Sox10 antibody and found highest Sox10 protein expression in the TN melanomas (Fig. 4E and F). These results indicate that more than one neural crest cell lineage specification program can be active in melanoma, perhaps reflecting activation of some of the recently described gene networks that characterize neural crest stage and fate decisions (28).

MITF is a lineage-dependent oncogene (29), and we have previously demonstrated that BRAF^V600E MITF-low melanomas depend on MITF activity for survival (9). We asked whether our three new melanoma models presented here also depended on MITF activity for survival. Indeed, tumors in all three melanoma models rapidly regressed by apoptosis when MITF was eliminated by switching to the restrictive temperature (32°C), and recurred at the permissive temperature (<26°C; Fig. 5A and B; Supplementary Fig. S5A and S5B). In all fish, the melanomas recurred at the same sites as the initial tumor site. These experiments indicate that MITF is a lineage survival oncogene for the bulk of melanoma with and without BRAF^V600E mutations, and in melanomas with MITF-low or MITF-high activity.

The dramatic loss of tumor followed by regrowth in the same site led us to hypothesize that a subpopulation of melanoma cells remains at the regression site and can resist loss of MITF activity. To address this, we crossed our melanoma-prone zebrafish (TS and TN) to zebrafish that carry a transgene that expresses GFP from the mitfa promoter in order to follow GFP-positive cells during melanoma regression. As explained previously, we chose this mitfa promoter because it is not regulated by Mitfa protein, and therefore the expression of BRAF^V600E or GFP does not change as a result of altering Mitfa protein activity in our mitfa^vc7 conditional allele (9). We found GFP expression was clearly detected in the bulk of the melanoma prior to regression (Fig. 5C).

Strikingly, despite no overt evidence of melanoma on the animal, confocal imaging revealed small clusters of GFP⁺ melanoma cells at the original tumor site following regression (Fig. 5D). Clusters of GFP⁺ cells were specific to the regression site and were not found in matched healthy tissue from the same animal (Fig. 5D), and non-melanoma melanocyte progenitor GFP⁺ cells were only sparsely located elsewhere on the fish skin (Supplementary Fig. S5C and S5D). Sectioning the fish at the melanoma regression site confirmed the presence of small clusters of GFP⁺ cells (Fig. 5E).

These results show that MITF is a lineage oncogene in BRAF-dependent and -independent melanomas, and reveal subpopulations of melanoma cells that remain at the regression site and are resistant to MITF loss. This residual disease cell population (termed MITF-independent cells) likely contributes to melanoma recurrence, and may represent an important new target for therapeutic interventions in melanoma.

We next sought to understand the nature of residual disease cells remaining at the tumor site after melanoma regression. To this end, we sorted small cell clusters (<50 cells) of GFP⁺ cells from primary tumor and residual disease from TS and TN growth melanomas (Fig. 6A; Supplementary Fig. S6A and S6B) and generated RNA-seq libraries using the Smart-Seq2 protocol (10). Analysis of primary and regressed cell clusters showed a number of differentially expressed genes with significant upregulation in the primary tumor (Fig. 6B). Expression of MITF target genes was strongly suppressed in melanoma cells residing at the regression site (Fig. 6C). Mitfa protein expression was almost undetectable at the regression site, whereas residual disease cells maintained the expression of Sox10 protein (Fig. 6D and E). GSEA analysis revealed that MITF-independent residual disease cells are enriched in mesenchymal and invasive gene expression signatures (e.g., as defined in refs. 30–32; Supplementary Tables S2 and S3), and in the neural crest stem cell (NCSC) signature identified in residual disease following BRAF plus MEK inhibitor treatment (Fig. 6F and G; Supplementary Table S3; ref. 33).

To establish whether these cell states are relevant to human melanoma, we analyzed an RNA-seq dataset of patients before BRAF plus MEK inhibitor treatment and following resistance (34, 35), and found that the MITF-independent residual disease cells share similar molecular expression signatures with human therapy–resistant melanomas (Fig. 6G). This enrichment is specific for melanoma samples because matched GFP⁺ melanocyte progenitor cells from unaffected tissue of the same fish are not enriched for these gene sets (Fig. 6G). Thus, the MITF-independent residual disease cell state is specific to melanoma, does not exist in non-melanoma melanocytes, and is not a result of the temperature-shift to turn off MITF activity.

Minimal residual disease cells are rare, and so we performed deep plate–based (Smart-Seq2) sequencing to investigate the residual disease state at the single-cell resolution. To investigate the residual disease state at the single-cell resolution, we profiled the transcriptomes of 332 GFP⁺ single cells from TS and TN growth melanomas and residual disease from corresponding subtypes, as well as from TN melanoma recurring disease (n = 43, 80, 63, 61, and 85 single cells, respectively; Supplementary Fig. S6C). We clustered cells based on similarities of their transcriptional profiles using the unsupervised Louvain algorithm on principal component reduction and visualized using Unifold Manifold Approximation and Projection (UMAP). We identified 6 clusters (numbered 0–5; Fig. 7A). Clusters 0 and 2 to 5 were primarily composed of cells from a predominant but not exclusive single experimental condition (Fig. 7B). In contrast, cluster 1 contained cells from all experimental conditions revealing a cell state common to all melanoma subtypes and stages (Fig. 7B).

All residual disease states had very low-to-no MITF target gene expression confirming the MITF-independent state of these cells. We investigated whether the MITF-independent cell state was specific to the residual disease cells or whether an MITF-independent state preexisted in the primary tumor. Strikingly, although MITF target genes were generally upregulated in the primary tumor cells (even at low levels; Fig. 7C), a subpopulation of cells in the primary tumors showed markedly reduced MITF target gene expression (Fig. 7C), indicating that the MITF-independent state is already present in the primary tumors. A preexisting MITF-independent subpopulation was also present in melanomas and residual disease from a smaller additional single cell dataset we generated (Supplementary Figs. S6D and S7A–S7B). All cells continued to express sox10 perhaps as a means to bypass the dependence on MITF activity for survival (Fig. 7D; Supplementary Fig. S7C). Thus, although the bulk of the tumor is dependent on MITF activity for survival, an MITF-independent state is present as a subpopulation and is resistant to loss of the MITF lineage factor.

Interestingly, we identified distinct cell populations in residual disease that were specific to melanoma subtype. Superficial melanoma residual disease cells (cluster 4) were enriched for invasive and mesenchymal genes, and expressed some of the genes of the Rambow-NCSC signature, compared with the superficial primary tumor (cluster 3; Fig. 7E; Supplementary Fig. S7D; Supplementary Table S8). In contrast, nodular melanoma residual disease (cluster 5) was enriched for genes involved in purine metabolism and the neural crest developmental program (Kaufman_zfin_NC; ref. 27) compared with the primary nodular tumor (cluster 0; Fig. 7F; Supplementary Table S8). Enhanced expression of and dependency on nucleotide metabolism is a feature of the neural crest and melanoma (36, 37). Expression of purine metabolism genes may represent a Sox10+ MITF-independent neural crest metabolic state. Alternatively, purine metabolism genes are highly expressed in xanthophores, a closely related zebrafish neural crest–derived pigment cell type (38), and expression of purine metabolism pathways may represent differentiation toward this neural crest cell fate to maintain survival.

Finally, we explored the cell state in nodular melanomas in the early stages of recurrence (10 days following reactivation of Mitfa). We found recurrent melanoma cells enriched for MITF target genes and protein translation genes (Fig. 7G; Supplementary Table S8). Despite the apparent shift away from the neural crest and purine metabolism pathways expressed in residual disease and toward MITF pathway activity, recurrent disease did not show gene expression associated with the cell cycle that is present in the primary nodular melanomas (Fig. 7G). Thus, early recurrent disease may need to first activate MITF transcriptional and protein translation pathways prior to activation of a proliferative phase to contribute to tumor regrowth.

Strikingly, all experimental conditions (primary, residual, and recurrent disease) shared a common cell state (cluster 1) that expresses neural crest genes including enrichment for erbb2 and zeb2a. Further, we found particular enrichment of genes recently described in a G₀-like neural stem cell state (39) such as the epigenetic factors arid1b, kmt2c, the cell-cycle inhibitor p27 (cdkn1bb), and ribosomal components (Fig. 7H and I). Some cells in cluster 1 also express high levels of the melanoma slow cycling cell marker JARID1B (encoded by the kdm5ba gene in zebrafish; ref. 40). These experiments show that residual disease generated through MITF loss leads to multiple MITF-independent cell states both exclusive and shared with primary tumors. Further, MITF depletion leads to a common neural crest G₀-like state (NC_G0-like_ribosomal) that exists in primary, regressed, and recurrent disease stages that may be vital for melanoma survival.

Gene expression signatures defined by the TCGA and Lund subtypes are predictive for survival across multiple independent melanoma datasets, indicating that gene expression profiling is a clinically relevant biomarker for patient stratification (2, 4–8). Moreover, immune-based transcriptional signatures can be predictive for response to immune therapies (41, 42). Of particular importance, the Lund MITF-low proliferation melanoma subtype and the TCGA MITF-low signature predict among the worst outcomes for patients (8). However, animal models of melanomas defined by these transcriptional signatures are lacking in part because little is known about how these subtypes are established.

Here, we have used a temperature sensitive conditional mitfa^vc7 mutant allele to model the MITF-low melanoma transcriptional state. We discover that the MITF-low state is highly invasive (mesenchymal) and exhibits plasticity toward a proliferative MITF-high state. Further, we discover a shared MITF-independent cell state in residual disease that preexists in primary melanomas. Together, these cell states may partially explain the poor outcomes for patients with MITF-low melanomas.

We define this MITF-low state in zebrafish by reduced levels of WT mitfa RNA transcript and Mitfa protein (Fig. 3; refs. 9, 20). The zebrafish MITF-low activity results in reduced expression of pigmentation genes and shares stem and mesenchymal gene expression patterns with the human TCGA MITF-low and Lund MITF-low proliferation gene expression clusters. Although the TCGA MITF-low transcriptional signature includes expression of neuronal genes, we demonstrate that the MITF-low state can also be characterized by a vasculogenic signature uncoupled from neuronal genes (i.e., DS melanomas). Taken together, zebrafish MITF-low melanomas are relevant models for the human MITF-low melanoma transcriptional gene expression cluster.

BRAF^V600E is not sufficient on its own to increase mitfa expression, and many BRAF^V600E-positive melanomas maintain a MITF-low gene expression signature. How mitfa expression is significantly increased in only some BRAF^V600E melanomas in our zebrafish models is not known, but our results demonstrate that nodular growth melanomas have the potential to increase sox10 and mitfa RNA and protein expression, and express genes similar to the Lund MITF-high pigmentation subtype. We term these MITF-high melanomas and find an increase in WT mitfa mRNA, and Mitfa protein expression is accompanied by an increase in BRAF^V600E expression because the BRAF^V600E transgene is under the control of a mitfa promoter fragment. Although we cannot differentiate between the impact of Mitfa and BRAF^V600E activity on nodular growth pathology, our models reveal the important finding that the MITF-low state is unstable with the potential to increase mitfa and Mitfa target gene expression.

Although melanomas can express high or low levels of MITF overall, cells can be heterogeneous for MITF within the tumor itself (43), and cells with low MITF are resistant to MAPK pathway inhibition (44). Single-cell RNA analysis indicates that cells that express high or low MITF RNA levels can coexist within an individual melanoma tumor supporting the concept that cells can switch their MITF transcriptional states in vivo (33, 34, 45–47). Our results suggest a Pax3a-Brn2b rheostat may regulate this switch in MITF gene expression from a low to high state specifically in BRAF^V600E melanomas (26, 43, 48). Understanding how MITF-low melanomas shift to MITF-high melanomas, and vice versa, is important because the Lund MITF-low proliferation transcriptional signature is enriched in melanomas post therapy, and this state confers intrinsic and acquired resistance to MAPK pathway inhibition in cell lines and in patient tumors (5, 44). Notably, patients can relapse on targeted therapy with individual melanoma clones that are both MITF-high and MITF-low, underscoring the importance of the transcriptional state switch in patients (44).

Our work shows that neural crest programs are active in both the superficial and nodular growth subtypes, but that the neural crest signatures that characterize the superficial growth melanomas (MITF-low melanomas) are enriched in mesenchymal and invasive genes (e.g., sox9b, twist 1a/b/3), whereas nodular growth melanomas (MITF-high melanomas) are enriched in developmental program genes (e.g., mitfa, sox10, foxd3, pax3a). Recent single-cell RNA-seq studies of the mouse neural crest reveal early migratory transcriptional states followed by branched fate transcriptional trajectories that coexpress bipotent cell fates prior to fate commitment (28). Although these states are not as markedly defined in our melanomas as found in development, our results extend the understanding of neural crest pathways in melanoma by showing that more than one neural crest gene expression program can be active depending on melanoma subtype.

MITF activity is considered a drug target because MITF is a lineage survival oncogene and turning off MITF activity has little impact upon other cell types (29). Regression followed by recurrence in our zebrafish MITF on-off-on models has similarities to the regression and recurrence of patient tumors following BRAF inhibitor therapy in which both diseases progress through undetectable minimal residual disease. Through live imaging and IHC of zebrafish at the melanoma regression site, we visualized small-cell clusters of residual disease. Our model predicts these cell populations would be resistant to MITF inhibitors. Further, it seems likely that these cells contribute to disease recurrence, although we cannot exclude that other melanoma cells exist at the regression site that do not express GFP and may be transformed in the context of the tumor regression microenvironment for expansion and growth (e.g., as reported in ref. 49). Lineage-tracing studies in this model system, and more broadly in the study of residual disease, will be critical to determine which cells at the residual disease site maintain malignancy and contribute to disease relapse.

Our results here show that residual disease has very low-to-no MITF transcriptional activity, very low-to-undetectable Mitfa protein, and is resistant to MITF as a lineage survival oncogene. Recently, Rambow and colleagues (33) identified four subpopulations of cells in residual disease following BRAF and MEK inhibitor treatment including invasive cells and NCSCs with "very low to undetectable" expression of MITF or its target genes. MITF-independent cells express genes that are enriched for the invasive and NCSC states, indicating that these signatures can define both the BRAF plus MEK inhibitor–resistant cells and MITF-independent cells. One explanation for this result is that MAPK pathway inhibitors lead to the direct loss of MITF activity in residual disease and therefore resemble the same cell state as experimental loss of MITF as seen in our models. An alternative explanation is that tumor regression itself (caused by BRAF/MEK inhibitors or MITF loss) changes the microenvironment such that the extant MITF-independent subpopulation undergoes a transcriptional shift toward invasive/NCSC programs (e.g., TS residual disease), or even toward differentiation programs (e.g., TN residual disease) depending on the primary tumor subtype.

The MITF-independent state that is shared between all disease conditions (cluster 1) maintains expression of neural crest markers sox10 and foxd3 with specific enrichment of erbb2 and zeb2a. Further, cluster 1 shows upregulation of ribosomal genes together with the cycle-cell inhibitor p27 (cdkn1bb) and arid1b, which together are enriched in a G₀-like state (39). Notably, a few cells in cluster 1 strongly express the slow-cycling cell gene jarid1b (kdm5ba; ref. 40). This G₀-like state may have important clinical implications for patients because G₀-like cells can be enriched following therapy (50) and have been reported to require JARID1B to exert a stem cell–like function (51). Lineage tracing experiments will be important to directly test if the extant cluster 1 cells in the primary tumor are the same cells in residual disease, and how they contribute to recurrent disease. Future cell transplantation studies will enable the characterization of phenotypic differences among residual disease cell subpopulations, such as tumor initiation or invasion potential, or differentiation states.

In summary, our results prove the MITF-low transcriptional signature reflects an underlying biological state and reveals melanoma transcriptional heterogeneity within this subtype, which may contribute to poor outcomes for patients. Further, given MITF is an important drug target, we reveal an extant melanoma cell subpopulation that resists loss of MITF activity and may contribute to disease perdurance. Some of the MITF-independent residual disease may be sensitive to therapy. For example, RXR-gamma inhibitors can target melanoma cells in the NCSC state (33), and mitochondrial respiration inhibitors can target the slow-cycling JARID1B^High cells (52). Our work shows further complexity for the targeting of residual disease because we find residual disease states that are shared by all disease stages and those that are exclusive to tumor subtype. These findings underscore the importance of targeting both the bulk of the tumor and residual disease cells for long-term survival outcomes.

No potential conflicts of interest were disclosed.

Conception and design: J. Travnickova, S. Wojciechowska, D.V. Brown, J.A. Lister, E.E. Patton

Development of methodology: J. Travnickova, S. Wojciechowska, A. Brombin

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Travnickova, D.V. Brown, T. Lefevre, A. Capper, R. Dilshat, T. Voet, E.E. Patton

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Travnickova, S. Wojciechowska, A. Khamseh, P. Gautier, A. Brombin, A. Ewing, M. Spitzer, R. Dilshat, C.A. Semple, M.E. Mathers, E. Steingrimsson, C.P. Ponting, E.E. Patton

Writing, review, and/or revision of the manuscript: J. Travnickova, S. Wojciechowska, A. Khamseh, A. Ewing, C.A. Semple, M.E. Mathers, J.A. Lister, E. Steingrimsson, T. Voet, C.P. Ponting, E.E. Patton

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J. Travnickova, A. Khamseh

Study supervision: C.P. Ponting, E.E. Patton

We are grateful to Yuting Lu for assistance with TCGA analysis, Craig Nicol and Connor Warnock for assistance with figure preparation, Elisabeth Freyer and Stacey Thomson for FACS analysis, Tamir Chandra and Angela Salzano for Smart-seq2 experimental advice, Lauren Saunders for single-cell analysis advice, and Cameron Wyatt and Wei Qing for zebrafish husbandry. We thank the following funders for supporting this work: S. Wojciechowska: CRUK PhD studentship; A. Khamseh is a cross-disciplinary postdoctoral fellow supported by funding from the University of Edinburgh and Medical Research Council (core grant to the MRC Institute of Genetics and Molecular Medicine); D.V. Brown: [PEGASUS]2 Marie Skłodowska-Curie Fellowship (12O5617N); C.A. Semple: MRC (MC_UU_00007/16); A. Ewing: UKRI Innovation Funds (MR/R026017/1); E. Steingrimsson: Icelandic Research Fund (184861-051); R. Dilshat: University of Iceland PhD studentship; T. Voet: Foundation Against Cancer grant (2015-143); C.P. Ponting: MRC (MC_PC_15075, MC_UU_00007/15); E.E. Patton: MRC (MC_UU_00007/9), European Research Council (ZF-MEL-CHEMBIO-648489), and L’Oreal-Melanoma Research Alliance (401181).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.