Microenvironment-derived factors driving metastatic plasticity in melanoma

Cellular plasticity is a state in which cancer cells exist along a reversible phenotypic spectrum, and underlies key traits such as drug resistance and metastasis. Melanoma plasticity is linked to phenotype switching, where the microenvironment induces switches between invasive/MITFLO versus proliferative/MITFHI states. Since MITF also induces pigmentation, we hypothesize that macrometastatic success should be favoured by microenvironments that induce a MITFHI/differentiated/proliferative state. Zebrafish imaging demonstrates that after extravasation, melanoma cells become pigmented and enact a gene expression program of melanocyte differentiation. We screened for microenvironmental factors leading to phenotype switching, and find that EDN3 induces a state that is both proliferative and differentiated. CRISPR-mediated inactivation of EDN3, or its synthetic enzyme ECE2, from the microenvironment abrogates phenotype switching and increases animal survival. These results demonstrate that after metastatic dissemination, the microenvironment provides signals to promote phenotype switching and provide proof that targeting tumour cell plasticity is a viable therapeutic opportunity.

P revailing models of metastasis have implicated changes in cellular differentiation as a key driver of the process. These models have suggested that undifferentiated/invasive cells initiate metastasis, and these cells then seed secondary organs to establish successful macrometastases. Two models for cell state transitions in melanoma have been proposed: the cancer stem cell model (hierarchical, irreversible) and the phenotype switching/plasticity model (reversible). An unanswered question prevails in the existing models: how do the slowly proliferative, invasive cells which arrive at new sites exit that state to then resume proliferation in secondary sites? In melanoma, these changes in cell state are regulated by the transcription factor MITF, which has been proposed to act as a molecular rheostat: cells with low levels of MITF are invasive, whereas those with higher levels of MITF are proliferative [1][2][3] . Because MITF also drives differentiation, one implication of this model is that cells which are ultimately successful in metastatic sites should be simultaneously proliferative and differentiated. As switches between the invasive/MITF LO versus proliferative/ MITF HI states are posited to be due to microenvironmental factors, this would paradoxically imply that microenvironments which promote the differentiated cell state would be most strongly associated with metastatic success 4 . The identity of such differentiation inducing factors in melanoma metastasis remains largely unknown, yet has important consequences in understanding the forces that drive macrometastatic colonization. To address this idea, we have used a zebrafish model of melanoma to monitor changes in cell differentiation during metastatic engraftment, and find that although early metastases are seeded by undifferentiated cells, over time these cells enact a gene program of melanocytic differentiation that is strongly associated with proliferation. The microenvironment induces this differentiated metastatic state in part through the developmental morphogen EDN3. Prevention of this differentiated cell state via EDN3/ECE2 inactivation doubles survival of the animals, suggesting that the acquisition of a differentiated cell state is necessary for metastatic success in melanoma, and can be prevented by interrupting microenvironment-melanoma cross-talk. Our data would imply that solely targeting undifferentiated cell populations is likely to miss an important component of metastatic lesions.

Results
Zebrafish imaging demonstrates that metastases differentiate. We utilized a zebrafish model of melanoma to monitor changes in differentiation during metastatic spread. We created transgenic zebrafish in which the melanocyte-specific mitfa promoter drives the human BRAF V600E gene 5,6 along with a mitfa-green fluorescent protein (GFP) reporter cassette. From this animal, we derived a zebrafish-specific melanoma cell line, ZMEL1-GFP, which can be transplanted into the transparent casper strain of fish and metastatic patterns visualized using in vivo imaging 7 . Similar to most human melanoma cell lines, the ZMEL1-GFP line was completely unpigmented in vitro and has a highly mesenchymal appearance (Fig. 1a, left). This allowed us to monitor differentiation using both the mitfa-GFP transgene as well as melanin pigmentation of the tumours. Whereas melanin in mammals is not a consistent marker of differentiation, zebrafish produce a very dark variant of melanin that is easily visualized, because fish carry a polymorphism of the pigmentation gene SLC24A5 that is typically associated with darker skin in individuals of African descent 8 . Upon transplantation into subcutaneous sites (Fig. 1a, top), akin to an in-transit metastasis in humans, we first observed that local tumours were mitf-GFP þ but completely unpigmented. With a latency of 7-14 days, 100% of the locally engrafted skin tumours regained pigmentation, a reflection of melanocytic differentiation. The secondary subcutaneous metastases that then developed were also initially mitf-GFP þ and unpigmented, but 53% of the subcutaneous metastases then became visibly pigmented by day 14 (Fig. 1a and Supplementary  Movies 1 and 2). This suggested that metastases are initially undifferentiated but became differentiated after engraftment. To confirm this, we directly transplanted B50 unpigmented ZMEL1-GFP cells into the vasculature of larval casper zebrafish, bypassing the initial skin site (Fig. 1a, bottom). We used larvae for this assay, rather than adults, due to the greater accessibility of the vasculature. Cells initially circulated widely, but within 24 hours extravasated in several areas including the skin, caudal haematopoietic territory and the eye. Over the next 7-28 days, these animals formed widespread macrometastases that were also visibly pigmented (Fig. 1a, bottom), the appearance of which is dominated by the superficial subcutaneous metastases. In both the adult and larval assays, this pattern of metastatic spread is analogous to Stage IV (Tx/Nx/M1a) disease. Histological analysis of a widely disseminated tumour confirmed clear areas of pigmentation in muscle invasive disease (Fig. 1c) adjacent to the overlying skin. To ensure that our findings were not unique to the ZMEL1 line, we performed similar experiments with an additional zebrafish melanoma cell line, ZCREST1. This line is derived from a BRAF;p53 tumour but the reporter GFP gene is driven by the crestin neural crest reporter, as previously described 9 . Similar to the to the ZMEL1 results, the initially engrafted ZCREST tumours were GFP þ but unpigmented, but had become pigmented by day 14-44 post-transplant ( Supplementary Fig. 2). These data are consistent with mouse studies showing acquisition of a pigmented phenotype after dissemination and engraftment 10 . Pigmentation of melanocytes typically coincides with the acquisition of a dendritic phenotype 11 . To further confirm that the initially undifferentiated ZMEL1-GFP cells could actively differentiate in vivo, we used time-lapse microscopy of cells injected into the vasculature of the fish to analyse their morphology. Within 57 hours post extravasation, the melanoma cells acquire a highly dendritic appearance characteristic of differentiated melanocytes, and by 5 days post transplant, 83% had acquired a dendritic appearance (defined as greater than two projections away from the cell body) compared with none on the day of transplant (Fig. 1c,d, Supplementary Fig. 1 for quantification of dendricity and Supplementary Movie 3 for time-lapse imaging movie), confirming that post-extravasation differentiation is an early step in metastatic engraftment.
Metastases exhibit a differentiation signature. To quantify the extent of this differentiation program, we performed RNA-seq on the metastatic ZMEL1-GFP cells as compared with the parental cells in vitro, generating a list of B2,000 genes that were altered after dissemination ( Fig. 2a and Supplementary Data 1). We found that canonical differentiation genes such as PMEL, TYR, TYRP1 and DCT were all upregulated in the disseminated cells (Fig. 2b), all known targets of MITF 12 . We found only a small increase in MITFB and no increase in MITFA (the two zebrafish orthologues of MITF), suggesting that it is increased activity, rather than expression, of MITF that was leading to the differentiation program. We also found a significant increase in markers of proliferation such as MYC and MYCL (Supplementary Fig. 3), consistent with prior observations that differentiation and proliferation are tightly linked in melanoma 2 , in part by regulating levels of MITF activity 13 . We used Gene set enrichment analysis (GSEA) to compare the zebrafish metastatic signature to a panel of known human differentiation genes (from ES-derived melanocytes) 12 . We found a significant enrichment between the disseminated zebrafish cells and human differentiated melanocyte genes (NES ¼ 1.276, FDR ¼ 0.009), including PMEL, TYR and TYRP1 ( Fig. 2c and Supplementary Data 2). We also used GSEA to compare our disseminated ZMEL1-GFP signature to a panel of human melanomas that have been previously classified as 'invasive' versus 'proliferative' 14 , and found that the metastatic zebrafish cells were strongly enriched for the 'proliferative' signature that is dominated by differentiation-related target genes of MITF such as PMEL, TYR, SLC45A2 and GPR143, (Fig. 2c and Supplementary Data 3, NES ¼ 1.317, FDR ¼ 0.08). Taken together, this data demonstrates that successful macrometastatic engraftment is associated with the acquisition of a differentiated and proliferative cell state.
Differentiation genes predict worse outcome in humans. This data suggests that the gene programs associated with metastatic cells phenotypically switching to a more differentiated state should predict prognosis in human metastatic melanoma. To examine this question, we generated a list of 83 differentiationassociated genes (that is, PMEL, TYR, TYRP1, DCT, SLC45A2, MLPHA, MLANA) derived from the GSEA analysis above as well as the ZFIN database (http://zfin.org/), and then assessed how their expression correlated with survival in the TCGA melanoma cohort. We stratified patients by those who presented with stage I (localized) versus stage III/IV (lymph node or distant metastases). In those presenting with stage I disease (even if they eventually went on to develop metastases, as did most patients in the TCGA cohort), none of the differentiation associated genes were associated with prognosis ( Fig. 2e and Supplementary Figs 4 and 5), which is consistent with prior reports indicating no relationship between pigmentation of primary tumours with prognosis 15 . In contrast, in stage III/IV disease, a cluster of genes indicative of differentiation state all portended a significantly worse prognosis: PMEL, EDNRB, TYR, GPR143, SLC45A2, MLPH and BACE2 ( Fig. 2f and Supplementary Figs 4 and 5). Examination of the Kaplan-Meier curves for the most strongly associated gene, PMEL (Fig. 2g,h), shows that patients with the highest level expression, and not intermediate or low levels, had the worst prognosis. This extends to other differentiation associated genes such as TYR and TYRP1 ( Supplementary Fig. 6). This data indicates that these pigmentation genes stratify those who present with advanced stage III/IV disease and is associated with macrometastatic growth once the cell has arrived at a secondary location.
Microenvironmental factors promoting plasticity. We wanted to identify the microenvironmental triggers of this differentiation program. We analysed our RNA-seq data using Ingenuity Pathway Analysis (IPA), including the canonical and upstream regulatory analyses. IPA makes statistical predictions (z-scores) of what upstream molecules or pathways are likely to have produced that gene signature. From this analysis, we identified 253 potential molecules that could yield the disseminated melanoma signature, which we then categorized into six pathways that we thought likely activated by extracellular/secreted factors: endothelin, dopamine, IGF, cholesterol, fatty acid and c-kit signalling (Fig. 3a, left panel). We predicted that if we applied agonists of these pathways to our cells in vitro, we would recapitulate the differentiated/proliferative state we had seen in vivo. Using the InCell6000 high-throughput/high-content imaging system, we performed a small molecule in vitro screen (Fig. 3b, middle and right panels and Supplementary Data 4) to identify molecules that could induce both differentiation (read out by elongation morphology) as well as proliferation (read out by cell number). We tested these in the zebrafish ZMEL1 melanoma cells, along with p53 À / À ) were orthotopically transplanted into the ventral skin of the transparent casper strain and then imaged using brightfield and GFP over 28 days. At days 1 and 7 post-transplant, the tumours are GFP þ but devoid of pigmentation, a marker of differentiation. By day 14 post transplant, the primary tumour mass in the orthotopic site has become deeply pigmented and the animal has developed small anterior metastases that are GFP þ but unpigmented (red arrowhead). By day 28, these metastases have now enlarged and are clearly pigmented (black arrowhead), consistent with metastatic differentiation. (a, bottom) ZMEL1-GFP cells were transplanted into the vasculature of larval casper recipients to assess direct differentiation capacity at sites of metastatic colonization, bypassing the primary skin site. Similar to what is seen in the orthotopic transplantation, the cells are initially unpigmented at days 1 and 7, but become increasingly pigmented at days 14 and day 28, indicating that cells can directly differentiate after extravasation. human A375 and SK-Mel28 melanoma cells to ensure that the effects we see are not fish-specific. In the ZMEL1 line (Fig. 3b), we found that 5/7 tested molecules (EDN1, EDN3, IGF1, SCF and LDL-cholesterol) significantly increased proliferation, whereas dopamine/L-DOPA caused a marked decrease in cell proliferation. In contrast, we found that only EDN1, EDN3 and Dopamine/L-DOPA caused a significant elongation phenotype (Fig. 3d) that the elongation morphology represented bona fide differentiation, we treated ZMEL1-GFP cells with either EDN3 or L-DOPA and then directly measured melanin content using spectrophotometric absorbance (Fig. 3f) and examination of the cell pellets (Fig. 3g). This confirmed that both of these molecules significantly increase melanin content/differentiation. In the human cells, we found that EDN3 caused a significant increase in both proliferation and elongation in A375 (Fig. 3c,e) and SK-Mel28 ( Supplementary Fig. 7) cells, whereas dopamine/ L-DOPA only caused a block in proliferation and did not cause the elongation phenotype that was seen in the fish cells. When taken as a whole, this data suggests that melanoma cells use a complex set of interconnected signals to enact the differentiated/proliferative state. This state relies on high levels of endothelins/SCF/IGF1 to promote proliferation, but low levels of dopamine/L-DOPA to allow for differentiation but not at a concentration that would adversely affect proliferation (Fig. 3h). Among these various factors, our data indicate that microenvironmental EDN3 is uniquely capable of inducing a phenotype switch to yield melanoma cells that are both proliferative and differentiated across both zebrafish and human conditions.
Loss of microenvironmental EDN3/ECE2 abrogates phenotype switching. Endothelins such as EDN3 are vasoactive peptides that are expressed in the zebrafish epidermis 16 , likely by the keratinocytes which normally surround the melanocytes 17,18 . The endothelins are known to act as key inducers of normal embryonic and adult melanocyte formation by binding to the EDNRB receptor 19 and have been previously implicated in melanoma growth by some but not all groups [20][21][22] . We hypothesized that inactivation of EDN3 in the microenvironment would abrogate phenotype switching in the nearby melanomas. To test this, we created CRISPR mutantrecipient zebrafish in which we introduced loss-of-function mutations in EDN3b, one of the two zebrafish orthologues of EDN3. In addition, we also created CRISPR mutants against endothelin-converting enzyme 2, ECE2b. This is the proteolytic enzyme recently shown to synthesize EDN3b in zebrafish 16 and has been shown to be one of the enzymes involved in endothelin synthesis in keratinocytes 23 . ECE2 was also one of the differentiation-associated genes we found associated with worse prognosis in human patients. We performed these experiments using the AB strain of zebrafish, rather than casper, because we wished to ascertain the strength of the EDN3b or ECE2b mutation itself on normal melanocyte development, and to ensure that the CRISPR was leading to a true loss of function ( Supplementary Fig. 8). We found that melanocyte development in both the EDN3b mutant and ECE2b mutant was severely abnormal (Fig. 4a-c), with loss of the dorsal and ventral melanocyte stripe, an overall reduction in the number of mature pigment cells, and a misshapen appearance of the residual melanocytes. This is consistent with a role for EDN3b and ECE2b in promoting melanocyte differentiation and growth 16,24,25 . We then transplanted ZMEL1-GFP cells subcutaneously into either WT, EDN3b or ECE2b-deficient recipients to ascertain the microenvironmental influence on melanoma proliferation and differentiation (Fig. 4d). Because in the wild-type (WT) background, we cannot reliably quantify distant metastases, we transplanted the ZMEL1-GFP sites into multiple subcutaneous sites, essentially to mimic a patient with multiple subcutaneous metastases. Melanomas that developed in both the EDN3b or ECE2b-deficient backgrounds were significantly smaller than those that developed in control animals ( Fig. 4e-g, quantified in 4h), suggesting that EDN3b supports melanoma proliferation in this site. Moreover, the tumours that developed in the EDN3b or ECE2b-deficient backgrounds were markedly less pigmented with tumours typically having a greyish appearance, and exhibited less mitf-GFP fluorescence compared with WT recipients (Fig. 4i), consistent with lack of differentiation. These phenotypic effects translated to a significant improvement in survival of the animals (Fig. 4j). At 50 days post-transplant, we found that only 6.25% of the WT animals had survived, in contrast to 42.1% of the ECE2b and 54.1% of the EDN3b-recipient animals (P ¼ 0.0035, log-rank test). To further ensure that these results were not unique to the AB strain, we performed additional transplants in the casper strain as well, and found that the EDN3b-deficient recipients similarly developed smaller tumours that were less pigmented ( Supplementary Fig. 9).
Collectively, these data demonstrate that microenvironmental EDN3b, acting via ECE2b, promotes phenotype switching of melanoma cells by simultaneously inducing a cell capable of both differentiation and proliferation. We find that abrogation of phenotype switching in the melanoma makes the tumours ultimately less lethal to the animal.

Discussion
The capacity for cells to successfully colonize secondary organs during metastatic spread has a strong temporal component. Cells that initially arrive at new sites are invasive, yet poorly proliferative, since these two programs tend to exist at opposite ends of a phenotype spectrum 26 . After initial extravasation, these nascent metastatic cells must inevitably exit this invasive state and begin to proliferate again, but how cells make that transition in melanoma has remained elusive. Our data would suggest that the reacquisition of a proliferative state in metastatic sites is tightly linked to the acquisition of a differentiation gene program. This likely reflects, in part, the dual capacity of the melanocyte master regulator MITF to induce both proliferation and differentiation. In support of this idea, we find that zebrafish melanomas harbour clear characteristics of a differentiated, pigmented cell state after extravasation. Most importantly, we find that human metastatic tumours that harbour a similar differentiation gene signature portend a much worse prognosis, likely because those tumours are more capable of metastatic proliferation. We have identified multiple microenvironmental factors inducing this differentiation gene signature. Our data  It is one factor that can promote a cell state capable of both proliferation and differentiation, which strongly promotes metastatic outgrowth. Recent studies have implicated enhanced endothelin signalling in the propensity of melanoma to metastasize to the brain, and EDN3 overexpression in mice leads to more metastatic and hyperpigmented melanomas 27 . This may be in part due to increased phenotype switching towards a differentiated state. The mechanism by which lack of microenvironmental endothelins lead to smaller, less pigmented tumours is likely due to a direct effect on tumour cells (via the EDNRB receptor on the melanoma cell) in addition to possible vascular effects from its role in endothelial cells (as evidenced by central necrosis in our tumours). However, clinical trials of endothelin antagonists have thus far been disappointing 28 . In part, this may be due to the relatively low potency of EDNRBspecific antagonists. We tested one such antagonist, IRL2500 (ref. 29), and found that it only inhibited melanoma growth at high concentrations that are unlikely to be clinically attainable ( Supplementary Fig. 10). The development of more potent EDNRB antagonists are being developed 30   illustrated in Fig. 6). For example, we find that dopamine (an intermediate in the synthesis of melanin) may promote differentiation, and a recent study using labelled betaCIT, a cocaine derivative that reflects the dopamine uptake transporter 31 , inadvertently identified an undiagnosed melanoma brain metastasis. The plasma ratio of L-DOPA to tyrosine is associated with increased metastatic propensity 32 . The mechanism by which L-DOPA leads to enhanced differentiation (at the cost of proliferation) remains to be determined in future studies, but may be related to its known role in synthesis of melanin itself since it is the substrate for the TYR gene product, tyrosinase. Second, while prevention of phenotype switching towards differentiation may counteract macrometastatic growth, it may paradoxically leave behind a pool of highly invasive melanoma cells untouched, with subsequent negative consequences for the patient. Recent work has identified Wnt/ROR kinase and TGFb signalling as potential microenvironmental supporters of this invasive state 10,33 . We would predict that simultaneous targeting of the microenvironmental forces that promote both the invasive (via Wnt/TFGb) and differentiated (via EDN3/IGF/SCF/ dopamine) cell states will be required for long term efficacy, a concept which remains to be explored in future studies.

Methods
Cell culture. All cells were maintained in DMEM with 10% FBS/pencillin/streptomycin/glutamine. The ZMEL1-GFP cells were derived as previously described 7 and came from a transgenic zebrafish with the mitf-BRAF V600E and p53 À / À background. The A375 and SKMel28 cells were obtained directly from the American Type Culture Collection.
Adult transplantation. The zebrafish melanoma cell line ZMEL1-GFP was used as the parental population. This line contains a transgene in which the mitfa promoter drives human BRAF V600E and also harbours a p53 À / À loss of function mutation. The GFP reporter cassette in the parental line is also driven by the mita promoter. For the transplantation, adult casper-recipient zebrafish were irradiated with a total of 30 Gy (15 Gy split dose on days À 4 and À 3 before transplant). In all, 300,000-500,000 cells were resuspended in 1-5 ml of 0.9X PBS, and then transplanted into either the ventral or dorsal skin (using a pulled glass micropipette). The fish were imaged within 1d to ensure transplantation success, and then serially imaged over a period of 30 days using brightfield and fluorescence imaging.
Larval transplantation. Adult casper fish were incrossed, and embryos were cleaned and dechorionated. At 2 days post fertilization, approximately 50 ZMEL1 cells resuspended in 1nl of 0.9 Â PBS were directly transplanted into the circulation via the Duct of Cuvier technique. The fish were imaged immediately after imaging and any fish that did not receive the proper number of cells was discarded. The fish were imaged as above until approximately day 30 post transplant.
No irradiation was used in this assay.
Imaging and image processing. Fish were imaged using an upright Zeiss Discovery V16 equipped with a motorized stage, brightfield, GFP and tdTomato filter sets. To acquire images, fish were lightly anaesthetized with Tricaine 4 mg ml À 1 . Images were acquired with the Zeiss Zen software v1, and the post image processing was done using MATLAB or ImageJ. Quantification of pigmentation in the larval stage images was done using brightfield imaging after thresholding for black versus white areas of the fish, and then calculated as the percentage area of the tail occupied by black pigmented cells, see paragraph further down for quantification of pigmentation in adult fish images.
RNA-seq analysis. Larval casper zebrafish were transplanted with B50 ZMEL1-GFP cells and allowed to form widespread, pigmented metastases for 18-21 days. At that point, the fish were completely disaggregated using Liberase, and a single cell suspension was made by resuspending the cells in DMEM with 2% FBS. Simultaneously, ZMEL1-GFP cells maintained in culture were trypsinized and resuspended in DMEM/2% FBS. Both populations underwent fluorescence-activated cell sorting sorting using an Aria system, and sorted directly into DMEM with 20% FBS to ensure higher cell viability. The cells were then pelleted and RNA was isolated using RNA isolation kits (Zymo Research). Because of small amount of RNA recovered from this process, all RNA samples underwent an RNA amplification step using the NuGen system, followed by Illumina library preparation and attachment of barcodes. Samples were sequenced using the HiSeq2500 with approximately 20 million reads per sample, using 50 bp single end reads. Reads from each RNA-Seq run were mapped to the zebrafish reference genome version danRer7 from the UCSC Genome Browser using GSNAP and quantified on the gene level using HTSeq and Ensembl version 75. Differential expression analysis was performed using DESeq2. The zebrafish gene symbols were mapped to their human orthologues using the DIOPT tool (http://www.flyrnai.org/cgi-bin/DRSC_orthologs.pl). Genes with a corrected P value of o0.05 were considered significant.
Gene set enrichment analysis. Comparison of zebrafish melanoma to human melanocytes: We used previously published data in which human ES cells were induced to become either neural crest or mature melanocytes 34 . Illumina microarrays were performed on pre-differentiation human emrbyonic stem cells, day 11 neural crest cells, and day 26 mature melanocytes. From this analysis, normalized expression values were generated, and ranked gene lists for neural crest versus melanocyte (compared with day 0 cells) were input into GSEA PreRank as RNK files. From the zebrafish RNA-seq data set, we created a genelist of either upregulated (2040) or downregulated (342) genes in the disseminated ZMEL1-GFP versus in vitro ZMEL1-GFP cells and these were input into GSEA as GMX files. For all situation in which there were duplicate gene identifiers, we collapsed these by selecting the most highly expressed of the probeset or transcript ID. We then ran GSEA PreRank using the default weighted statistic. Comparison of zebrafish melanoma to human melanomas: We used previously published data in which human melanomas were categorized into invasive versus proliferative subtypes 3,14 . From this data, we created a ranked list of genes in either category, which were then input into GSEA as RNK files. The zebrafish data was the same as that described above. GSEA PreRank was similarly run using the default parameters.
Human TCGA survival analysis. RNA-sequencing and clinical data for cutaneous melanoma patients was downloaded using the R package TCGA-assembler 35  Ingenuity pathway analysis. From the RNA-seq analysis, the top 1,000 up or downregulated genes were used as input into IPA (http://www.ingenuity.com/ products/ipa). The Core Analysis module was selected, using Genes and Endogenous Chemicals. Only genes with a significant FDRo0.05 were used for the analysis. For identification of likely extracellular factors causing the gene expression signature, we used both the Canonical Pathways as well as Upstream Regulator analysis modules, with the P-values shown in Fig. 4 derived from this analysis.
Drug screening assay. For convenience, we combined a 1:1 mixture of 25,000 ZMEL1-GFP along with ZMEL1 cells overexpressing a SLUG-tdTomato transgene (for an unrelated question) and seeded these in black wall, clear bottom 96-well plates (Corning, catalog #3340) in 100 ml of DMEM/10% FBS. For this manuscript, the two cells were analysed separately and only the ZMEL1-GFP cells used. The cells were allowed to adhere overnight, and then media changed to drug containing media at the indicated concentrations for each agent. These were done at either low (0.5%) or high (10%) serum conditions. Cells were grown for a total of 4 days after the addition of drug, with media/drug refreshed at day 2. At the end of the treatment, the cells were fixed with paraformaldehyde (4%), and stained using Hoechst 33342 (4 g ml À 1 ).
The agonists in the screen described in Fig. 3a  High-content image analysis of ZMEL cells. High-content image analysis was used to measure the effects of pharmacological agents on ZMEL cell proliferation and differentiation. Using a GE INCell 6000, 2 fields were captured per well using a Nikon 10 Â Plan Apo objective, 0.45NA. Images of the Hoechst signal was captured via the 4,6-diamidino-2-phenylindole channel with 0.1 and 0.3 seconds exposure in order to identify the nucleus and cell body respectively. The FITC channel was used for GFP with a 1 second exposure, and the dsRed channel used for tdTomato, with a 2 second exposure. Using GE Workstation software, nuclei were identified using a top-hat segmentation of the first (shortest exposure time) Hoechst image and cell bodies using a global threshold of the second (longer exposure time) Hoechst image, with identified nuclei as seed objects. GFP-positive cells were classified as having an average intensity within the nuclear mask of 4500 in the GFP channel and tdTomato-positive cells with an average intensity 4360 in the dsRed channel. Measurements were captured for all cells as well as GFP and tdTomato cells separately. Cell elongation is defined as the ratio of the minor to major axis of the cell. Major axis is the longest straight-line interval that can be drawn in the cell mask, while the minor axis is the longest straight line perpendicular to the major axis.
High-content image analysis of human cells. A 1:1 mixture of 2500 A375 cells (RFP positive) and SK-Mel28 cells (no fluorescent protein) were seeded together in 96-well plates and treated under the same conditions as the Zebrafish cells except now low serum conditions were at (0.5%) and high (5%).
Using a GE INCell 6000, 2 fields were captured per well using a Nikon 10 Â Plan Apo objective, 0.45NA. Images of the Hoechst signal was captured via the 4,6-diamidino-2-phenylindole channel with 0.1 and 0.3 seconds exposure in order to identify the nucleus and cell body, respectively. The dsRed channel used for RFP, with a 2-second exposure. Nuclei and cell bodies were identified as per the Zebrafish cells. Cells were classified as either RFP positive (A375; having an average intensity within the nuclear mask of greater than 283 in the dsRed channel) or negative (SK-Mel28). Measurements were captured for all cells as well as RFP positive and non-fluorescent cells separately.
Composite Z-score used in heatmaps. Our composite Z-score is similar to the one defined in (ref. 36). The Z-score for a measured quantity (for example, the number of cells) x i in one replicate drug treated well is defined as: where m c and s c is the mean and standard deviation of the measurements in the control wells with the exact same conditions as in the i'th replicate well except no drug is present. (In our case these are control wells containing same vehicle and on the same plate, with cells from the same reservoir as the drug treated well). Thus, each replicate well receives a signed Z-score corresponding to the number of standard deviations it falls above or below the mean of a well-defined mock-treatment distribution. The composite Z-score is then defined as: where u¼ u 1 ; u 2 ; :::; u i ; :::; u n ½ , (and u i ¼1= ffiffiffi n p for all i) is a vector of unit length in R N and Z¼ Z 1 ; Z 2 ; :::; Z i ; :::; Z n ½ is a vector where the components are the Z-scores for n individual replicate measurements from n wells.
In the limit of perfect reproducibility all replicate wells have the same Z-score Z 1 ¼Z 2 ¼:::¼Z i ð Þ and Z and u are parallel. Z comp is the length of the projection of Z onto the line spanned by u -the direction in replicate space, which represents perfect reproducibility.
Melanin content assay. Treated cells were grown in standard tissue culture flasks and then trypsizined. Equal numbers of cells from each condition (that is, dimethylsulphoxide versus drug treated) were pelleted at 500 g for 5 min, then resuspended in 200 ml of 1N NaOH and heated at 80°C. The pellets were vigorously pipetted to homogenize them, and then the cell extract transferred to a 96-well plate to be read at 405 nM. On each day of experiments, a standard curve of melanin using commercially available purified melanin (Invitrogen) was used as a control. Values are expressed relative to the dimethylsulphoxide control treated cells.
Creation of EDN3b/ECE2b CRISPR mutants. We designed multiple guide RNAs (gRNAs) against the EDN3b or ECE2b coding sequence, and synthesized these via in vitro transcription. The gRNAs were incubated with purified recombinant Cas9 protein (PNA Bio), and these were coinjected (75 ng ml À 1 of RNA and 500 ng ml À 1 of protein) into single cell zebrafish embryos. A portion of the injected clutch was separated at 24 hpf and genomic DNA extracted. PCR amplification of the region flanking the expected cut site was performed, followed by the T7E endonuclease assay to identify those gRNAs which gave the expected cut size. The remaining fish from that clutch were raised to adulthood, and the T7E assay repeated to identify founders capable of passing on the mutation in the germline. After identification of two such founders, these were then incrossed to yield trans-heterozygous F1 fish which were examined for phenotype of the endogenous melanocytes. Because these F1s will be a combination of WT, true heterozygotes and transheterozygotes, we only performed transplants into fish that had clearly abnormal melanocyte patterns. The control WT transplants were done into a separate group of fish that had been injected only with Cas9 protein and no gRNA and 100% of cases had normal melanocyte development. The transplants in all these animals was done identically to the casper transplants except that the total dose of irradiation was 20 Gy and the cell number was 300,000 cells per animal. The control, WT fish were handled in an identical manner except that they were injected with purified Cas9 protein alone with no gRNAs. The relevant sequences: EDN3b gRNA: 5 0 -GGATAAATGTACTCACTGTG-3 0 F primer for genotyping: 5 0 -TCTCCAGTCTTCTCTGGTGGTT-3 0 R primer for genotyping: 5 0 -GTGTGACAGCGAAAGAGTAACG-3 0 ECE2b gRNA: 5 0 -GGCTCCGGCTATTGAATGGT-3 0 F primer for genotyping: 5 0 -TGACGCATATGTGAAGGTAAAAA-3 0 R primer for genotyping: 5 0 -AACCACAGAGACAGCCATACCT-3 0 Measurement of tumour size in WT and mutants. Adult fish were transplanted subcutaneously as described above, except that fish received either 1 or 2 transplants, since metastases are difficult to visualize in the AB background (as opposed to casper). Each fish was imaged using brightfield and GFP imaging on both the left and right side at 1-40 days post transplant (dpt). WT, EDN3b and ECE2b fish were imaged on a given day using identical exposure times and zoom. Each image was imported from the Zeiss Zen software as TIFFs and analysed in Image J. Each tumour was segmented using the default ImageJ segmentation algorithm and thresholding applied uniformly across all groups. Multiple measurement parameters were extracted, and we found highly consistent results between tumour 'area' and tumour 'perimeter', as expected. The two groups were compared with each other using unpaired ANOVA.
Kaplan-Meier analysis. All animals were followed for a period of 40-45 days and survival measured using the Kaplan-Meier method. The differences between the WT, EDN3b and ECE2b-recipient fish were analysed using the logrank statistic.
Data availability. The RNA-seq data used in this study is available via the NCBI GEO repository with accession at: http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE90143. Any other data is contained in this manuscript or in the Supplementary Information or available from the authors upon request. The human TCGA survival is derived from the publically available cBIO portal (https://cbiologin.mskcc.org).