A novel mouse model demonstrates that oncogenic melanocyte stem cells engender melanoma resembling human disease

Melanoma, the deadliest skin cancer, remains largely incurable at advanced stages. Currently, there is a lack of animal models that resemble human melanoma initiation and progression. Recent studies using a Tyr-CreER driven mouse model have drawn contradictory conclusions about the potential of melanocyte stem cells (McSCs) to form melanoma. Here, we employ a c-Kit-CreER-driven model that specifically targets McSCs to show that oncogenic McSCs are a bona fide source of melanoma that expand in the niche, and then establish epidermal melanomas that invade into the underlying dermis. Further, normal Wnt and Endothelin niche signals during hair anagen onset are hijacked to promote McSC malignant transformation during melanoma induction. Finally, molecular profiling reveals strong resemblance of murine McSC-derived melanoma to human melanoma in heterogeneity and gene signatures. These findings provide experimental validation of the human melanoma progression model and key insights into the transformation and heterogeneity of McSC-derived melanoma.

M elanoma is the deadliest form of skin cancer. The current model of human melanoma progression proposes that melanomas first arise in skin epidermis during the radial growth phase, and then invade into the dermis for continued expansion during the vertical growth phase 1 . In this latter phase, melanoma cells can establish a heterogeneous tumor composed of subpopulations with distinct proliferative, biochemical and metastatic signatures 2 . Prognosis of melanoma becomes significantly worse once the vertical growth phase is initiated, and subsequent tumor heterogeneity makes the design of effective therapeutic approaches, particularly targeted therapies, increasingly difficult 3,4 . Therefore, it is important to understand the earliest events in melanoma initiation and progression to develop strategies for early detection and intervention. Mice serve as the primary pre-clinical models in cancer research. However, to date, the validity of mouse models in studying melanoma has been questioned due to the fact that mouse melanoma arises in the dermis, while human melanoma arises specifically in the epidermis.
The capacity of melanocyte stem cells (McSCs) to create melanomas remains controversial 5,6 . To date, the study of Moon and co-workers has demonstrated the potential of hair follicle McSCs to form tumors in the well-known Tyr-CreER; Braf CA/+ (Braf V600E/+ ); Pten fl/fl (Tyr-CreER:Braf:Pten) murine melanoma model 5,7 , whereas the study by Kohler et al. 6 , using the same mouse, demonstrated their lack of tumor-forming capacity. Because Tyr-CreER can target both McSCs located in the hair follicle and melanocytes (Mcs) in the dermis 8,9 and melanoma forms primarily in the dermis of these mice 7 , it has proven difficult to conclusively establish the origin of melanoma using this model. Another melanoma mouse model, constitutively expressing hepatocyte growth factor/scatter factor (HGF/SF) for the migration of melanocytes to the epidermis, develops melanoma at the dermo-epidermal junction upon ultraviolet (UV) irradiation [10][11][12][13] . Although this model is thought to share more histopathologic features with human melanoma, it also cannot distinguish between epidermal and dermal melanocytes as a source for melanoma formation. Investigation for a putative vertical growth phase from epidermal melanoma in mouse melanoma studies has also been stymied using these models.
A major difficulty in the treatment of melanoma derives from the multiple levels of heterogeneity of this disease 14 . Complex phenotypic heterogeneity even within a single melanoma is common, in part because melanoma cells can dynamically and reversibly switch between differentiated and undifferentiated states, exhibiting distinct proliferative, invasive and tumorinitiating characteristics [15][16][17][18] . Without a precise understanding of the cell of origin, it remains impossible to delineate how a defined population of normal cells can initiate a transformation process that ultimately gives rise to a heterogeneous tumor. It has long been proposed that cancer cells can recapitulate embryogenesis, thus differentiated cells may acquire the multipotency of their embryonic ancestors to create heterogeneous tumors 19 . Without understanding a cellular origin of a particular melanoma, it remains impossible to test if and how this occurs after normal melanocytes acquire oncogenic mutations.
While oncogene activation and tumor-suppressor gene inactivation are thought to be the main driving events for the transformation of normal somatic cells into malignant tumor cells, the microenvironment has also been considered an active player in tumor initiation and niche signals have been shown to influence transformation in other types of cancer. For example, Wnt signal activation, driven by paracrine ligands, are required for maintenance and renewal of intestinal stem cells, but also promote their transformation during tumorigenesis 20,21 . Notch signaling, required for the proper renewal and differentiation of intestinal epithelium, is also a requisite for intestinal cancer initiation [22][23][24] . However, potential regenerative niche signals that synergize with oncogenic mutations to promote the transformation of normal melanocytes into melanoma remain unknown.
In this study, we generate a c-Kit promoter-driven model for melanoma induction 25 . We show c-Kit expression defines McSCs in the hair follicle (HF) and c-Kit + follicular McSCs give rise to epidermal melanoma upon combined oncogenic Braf V600E induction and Pten loss with response to normal niche signals for anagen onset. Using this model, we also clearly demonstrate that epidermal McSC-derived oncogenic Mcs invade the underlying dermis in a vertical growth phase to establish a heterogeneous melanoma remarkably similar to human melanoma.

Results
The c-Kit promoter defines follicular McSCs. To test the ability of the c-Kit promoter to target McSCs from the hair follicles away from the dermal melanocytes in the skin, we generated c-Kit-CreER; R26R-mTmG (c-Kit-CreER: R26R-GFP) mice in which membrane-bound GFP is expressed by c-Kit + cells and their progeny following tamoxifen (TAM) treatment 26 . 7-week old mice were injected with TAM for 3 days, during a period when hair follicles in the back skin are in telogen phase and only harbor McSCs 27,28 (Fig. 1a). McSCs are defined as dopachrome tautomerase (Dct) + cells in the bulge/secondary hair germ (sHG) of the telogen hair follicle 27,28 . C-Kit and Sox10 have also been reported as McSC markers [29][30][31] . Immunohistochemistry revealed GFP expression in about 70% Dct + McSCs in the bulge/sHG niche of the HF (Fig. 1b, c). These labeled cells persisted for >1 year ( >500 days) and were able to give rise to mature melanocytes ( Supplementary Fig. 1a-c), verifying the ability of the c-Kit promoter to target long-lived McSCs. Immunohistochemistry revealed that GFP + cells in the HF also expressed c-Kit and Sox10 (Fig. 1b). Although GFP expression was also occasionally detected in the dermis, none of the GFP + dermal cells expressed melanocyte and/or melanoma markers, including Sox10, S100b, and Nestin ( Fig. 1b, d, e) [32][33][34] . Rarely, GFP + CD45 + cells were observed in the interfollicular epidermis and dermis, consistent with the known expression of c-Kit in cells of hematopoietic lineage, however, the work of others has shown that this c-Kit-CreER line is not suitable for targeting hematopoietic stem cells (HSCs) because of low expression ( Supplementary Fig. 1d, e) 35,36 . GFP expression was also occasionally detected in Keratin14 + keratinocytes in the interfollicular epidermis ( Supplementary  Fig. 1e). None of the GFP + epidermal cells expressed Dct, consistent with the previous observations that epidermal melanocytes do not reside in the back skin of mice 28 . To further confirm that c-Kit-CreER does not target dermal melanocytes, we crossed c-Kit-CreER; R26R-Tomato reporter mice to Dct-rtTA; tetO-H2B-GFP mice, to GFP tag Dct + cells upon doxycycline treatment (Fig. 1f) 37 . As expected, we observed no Tomato + GFP + cells in the dermis of these mice (Fig. 1g, h), confirming that the c-Kit promoter targets only follicular McSCs.
To further verify this observation, we also examined Tyr-CreER;Braf CA/+ ;Pten fl/fl ;R26R-Tomato (Tyr-CreER:Braf:Pten: Tomato) mice following TAM treatment. We experimentally induced anagen onset in telogen mice and then treated the mice with TAM for 3 days (Fig. 2a). We observed pAKT expression and increased proliferation in follicular McSCs ( Supplementary  Fig. 3g-i) similar to c-Kit-CreER:Braf:Pten McSCs (Fig. 2b, c), with progression of epidermal oncogenic Mcs into epidermal melanoma cells that was never observed in control mice ( Fig. 2h; Supplementary Fig. 3j). It is noteworthy that the emergence of epidermal melanoma cells was associated with the robust and concomitant expansion of Tomato + dermal melanoma cells in the papillary dermis, as previously observed 5,7 . Both epidermal and dermal Tomato + cells expressed melanocyte/melanoma markers S100b, Dct, MITF and Sox10 ( Fig. 2h; Supplementary  Fig. 3e). The different phenotypes of Tyr-CreER to c-Kit-CreER driven mice suggested that dermal melanocytic subpopulation targeted by Tyr promoter may be responsible for the formation of dermal melanoma in Tyr-CreER:Braf:Pten mice.
Nevertheless, to ensure that Tyr + McSCs were responsible for interfollicular epidermal expansion, we injected BrdU during anagen onset to specifically label cycling McSCs but not quiescent dermal Tyr + cells and then followed with TAMinduction ( Supplementary Fig. 3k) 27,41 . BrdU + /Tomato + label retaining epidermal melanocytes were observed, verifying the origin of melanoma Mcs as follicular McSCs (Supplementary Fig. 3l-n).  Fig. 4h). Constitutive Wnt activation in McSCs in the absence of Braf:Pten mutations had no effect, suggesting that synergy with other pathway is needed for Wnt to activate McSCs ( Supplementary Fig. 4i).
Intriguingly, whereas constitutive Wnt signaling was sufficient to induce the expansion of Braf:Pten mutant McSCs during stem cell quiescence, this treatment did not promote their upward migration. However, overexpression of an additional niche factor, Edn1, in c-Kit-CreER:Braf:Pten:Tomato:β-cat-STA crossed to K14-rtTA; tetO-Edn1-lacZ mice 46,47 , promoted both expansion and upward migration of Ki67 + pAKT + McSCs in telogen ( Fig. 3f-j), which is normally observed only during anagen. When we exchanged c-Kit-CreER for Tyr-CreER, epidermal melanoma formation was also observed concurrently with expansion of Tomato + dermal cells (Fig. 3k-n). Of note, overexpression of epithelial Edn1 alone during telogen did not induce proliferation, AKT pathway activation or migration to the epidermis in Braf:Pten McSCs within the quiescent follicular microenvironment ( Supplementary Fig. 4j, k). These data show that Wnt and Edn signals, which are normally produced by epithelial cells during hair growth, provide the niche stimuli to promote follicular Braf:Pten McSCs to transform into melanomas. Thus, normal extrinsic signals can synergize with oncogenic mechanisms to promote McSC-derived melanoma initiation.
Epidermal melanoma cells invade the dermis. The current model of melanoma initiation in human skin posits that melanoma arises and expands first in the epidermis (radial growth phase) before it invades into the dermis (vertical growth phase) 1 . However, the biology of dermal invasion by melanoma Mcs remains poorly understood, in part because prior mouse lineage tracing models could not distinguish between epidermal and dermal Mcs as the source of melanoma 8 . Zebrafish scale models, which can undergo a form of epidermal melanoma, also provide an imperfect likeness to the human condition because they lack strong resemblance to human skin 48 . In contrast, the c-Kit-CreER:Braf:Pten melanoma model can overcome these limitations by targeting only McSCs in the epithelium.
To determine if and how epidermal melanoma cells invade and expand in mouse dermis, we examined skin from induced c-Kit-CreER:Braf:Pten:GFP mice after establishment of epidermal melanoma. At later time points (26 days after TAM), we occasionally found foci of dermal melanoma expansion beneath the epidermis, suggestive of dermal invasion by epidermal melanoma cells ( Supplementary Fig. 5a, b). Although we found very little GFP labeling in CD45 + cells in the skin of these mice ( Supplementary Fig. 5c), most animals died at this point, likely due to defects in c-Kit + cells in other organs, such as interstitial cells of Cajal in the gastrointestinal tract (Sun et al., unpublished).
To reliably follow the phenotypic evolution of epidermal melanoma cells in mice for longer time periods, we performed full-thickness skin grafting experiments. To do this, we first treated c-Kit-CreER:Braf:Pten:GFP mice with TAM for 3 days at 3 weeks of age to genetically induce Braf:Pten mutations and GFP expression in McSCs. Skin from these animals was then grafted onto the backs of immune-deficient J/NU mice (Fig. 4a). After 10 days of engraftment, we observed a horizontal expansion of epidermal melanoma cells (Fig. 4b), similar to that seen in the c-Kit-CreER:Braf:Pten skin in situ (Fig. 2d). By 17 days, GFP + melanoma cells were detected invading the graft dermis (Fig. 4c). This was observed in the interfollicular as well as follicular epithelium. Some melanoma cells, expanding within the hair follicle outer root sheath appeared to invade the underlying dermal compartment ( Supplementary Fig. 6a). The GFP + cells in the dermis were found to downregulate E-cadherin expression and upregulate MCAM expression (Fig. 4d), which is consistent with our and previous observations of human melanoma at the invasive radial growth phase ( Fig. 4e; Supplementary Fig. 6b) 49,50 .
Further investigation revealed that GFP + melanoma cells that transited into the dermis maintained their pigmentation and expression of melanocytic genes, such as Dct, a melanocytic marker ubiquitously expressed in the McSC lineage in normal mice and MITF, the master transcription factor of melanocytic genes ( Fig. 4f) 42,51 . These melanoma cells lacked expression of neuronal markers Nestin, GFAP and Tubb3 (Fig. 4f). By 30 days after engraftment, we found that dermal melanoma cells had significantly downregulated the expression of Dct and MITF as well as pigmentation (Fig. 4g, h). In addition, most GFP + melanoma cells now expressed Nestin, GFAP and Tubb3 (Fig. 4g, h), reminiscent of human melanomas known to frequently express these neuronal markers 32,52,53 . These results indicate that melanoma cells lose their melanocytic phenotype during tumor progression and they begin to express genes including common neuronal markers (Fig. 4i).
Because nude mice lack a proper immune system, which may affect melanoma progression, we also implanted skin from c-Kit-CreER:Braf:Pten:Tomato mice onto wildtype syngeneic littermates ( Supplementary Fig. 7a). These littermates formed dark melanoma by 26 days after implantation ( Supplementary Fig. 7b). and dermis in the implanted area ( Supplementary Fig. 7c). Similar to the results in nude mice, Tomato + cells in these tumors were positive for melanoma markers such as Sox10 and heterogeneously expressing melanocytic markers Dct and pigment, as well as neuronal markers Tubb3 and Nestin (Supplementary Fig. 7c). Moreover, we detected metastasis of Tomato + melanoma cells to the lymph node ( Supplementary Fig. 7d). Some of these cells expressed Dct and pigment, demonstrating their melanocytic signature. No metastasis to the lung or brain was detected at this time point.
To further demonstrate that c-Kit + McSCs can produce invasive melanoma, we topically treated c-Kit-CreER:Braf:Pten: Tomato mice with 4-hydroxytamoxifen (4HT-TAM) to induce transgene expression (Fig. 4j). By 47 days, these mice also formed heterogeneous primary melanoma tumors and lymph node metastases similar to those observed in skin transplants (Fig. 4k, l). We found very few reporter labeled CD45 + hematopoietic lineage cells in the dermal melanoma tumor ( Supplementary Fig. 8). These topically induced mice survived at least 80 days after TAM treatment, similar to Tyr-CreER driven mice (Fig. 4m).
Our results suggest that this switch from a melanocytic to a neuronal-like phenotype in these mice is largely a spatially and temporally regulated process and that it coincides with the transition from the radial to the vertical growth phase.
McSC-derived and human melanomas share molecular signatures. Our lineage tracing and phenotypic analyses suggest that the initiation and progression of murine McSC-driven melanoma resemble models of human melanoma proposed from clinical observations. To begin to understand the molecular mechanisms responsible for the transformation of McSCs to melanoma, we compared the transcriptional profiles of McSCs and their derived melanoma.
To obtain tumors derived from McSCs, we isolated cells from uninduced skin epidermis including hair follicles from either c-Kit or Tyr oncogenic models and transplanted them with wildtype support cells into immunocompromised nude mice (Fig. 5a) 54 . Immediately following transplantation, we induced Braf:Pten mutations and Tomato reporter expression by TAM treatment. Transplanted McSCs from both models formed visible tumors typically by 47 days after transplantation (Fig. 5b). Consistent with a previous report, transplanted McSCs from control mice can enter the hair follicle to become both quiescent cells in the bulge and pigmented melanocytes in the bulb, but never formed tumors ( Fig. 5c; Supplementary Fig. 9) 54 . Tumors from both models were indistinguishable, both containing high numbers of Tomato + tumor cells expressing both S100b and Sox10 (Fig. 5d) 55,56 (Fig. 5e, f). Gene set enrichment (GSEA) analysis also revealed a positive correlation between tumorenriched genes and a human neural crest gene set 57 (Fig. 5f), similar to that described in a zebrafish melanoma model 48 .
Consistent with our 30 day in vivo phenotypic analysis (Fig. 4g), neuronal markers including Nestin, Tubb3, Gfap and Ngfr were upregulated in tumors, again a noted feature of human melanomas ( Fig. 5e-g). This was coupled to the downregulation of melanocytic genes (Fig. 5g). Immunohistochemical analyses with selected targets confirmed our RNA-seq results (Fig. 5h). These analyses confirmed that mouse McSC transformation involved downregulation of melanocytic genes, upregulation of ECM components, melanoma markers and neuronal/neural crest genes, all changes reminiscent of those observed in human melanoma.  (Fig. 5a). Cells segregated into 7 clusters with unsupervised clustering using Seurat analysis (Supplementary Fig. 11a-c). Two very small clusters include contaminating macrophages and keratinocytes and were excluded from the analyses (Supplementary Fig. 11d, e). Among the 5 major clusters, we found that wildtype McSCs (red) clustered into one population, suggesting that McSCs are relatively homogeneous compared with McSC-derived melanoma  (Fig. 6a, b). Little overlap was observed between normal McSC and tumor groups indicating significant transcriptional changes. Overall, all 4 melanoma clusters showed reduced expression of melanocytic markers like Dct, MITF and Pmel, along with elevated expression of melanoma markers (S100b, MCAM and Sox10) and neuronal and neural crest markers (Nestin, Tubb3 and Ngfr) (Supplementary Fig. 12a). Our single cell RNA-seq analyses also revealed that enrichment of signatures identified by bulk analysis, including ECM-interaction pathway, neural crest signatures and neuronal  (Table 1). We also found that Clusters 2 and 3 were most enriched in a larger set of neural crest/neuronal signatures, such as Foxd3 and L1cam ( Fig. 6c; Supplementary Data 1). Cells in Cluster 3 also highly expressing cell cycle genes ( Fig. 6d; Supplementary Data 1). By assigning cell cycle scores to each cell based on its expression of S and G2/M phase markers, we found that only Cluster 3 was enriched in cells with high S and G2/M scores ( Supplementary  Fig. 12b, c), suggesting that this small population of tumor cells could be responsible for tumor growth, presumably giving rise to cells in at least some of the other clusters. Cluster 5, with the lowest expression of neuronal genes, was most enriched in mesenchymal gene signatures ( Fig. 6e; Supplementary Data 1). Cluster 4 appeared to be an intermediate population, sharing 124 genes with cluster5 and 73 genes with cluster2 of the top 400 enriched genes in each cluster ( Fig. 6f; Supplementary Data 1). Immunohistochemical examination of 10 distinct tumors showed that in most of the tumors, Tomato + tumor cells express melanocytic markers Dct and MITF, neuronal markers Nestin, Tubb3 and Gfap, mesenchymal marker Pdgfra and proliferation marker Ki67 (Supplementary Fig. 12d).
Interestingly, compared to clusters 2 and 3, clusters 4 and 5 melanoma cells were enriched in a gene set previously identified in metastatic melanoma cells from patients ( Supplementary  Fig. 13a) 58 . This set includes genes related to stress response, invasion and resistance to treatment with MAPK inhibitors 58 . Of note, one gene in this cluster, SerpinE2 has recently been characterized as a marker of slow cycling cells in human melanoma that is most metastatic 59 . We found that most melanoma cells expressed Axl, a key marker gene for resistance to targeted therapy in human melanoma, whose expression is mostly exclusive from Mitf + cells, consistent with previous human melanoma report ( Supplementary Fig. 13b) 58,60 . Wnt and Edn pathway related genes were expressed by very few melanoma cells, suggesting that these pathways may not be required for tumor maintenance (Supplementary Fig. 13c).
To understand relatedness between clusters, we reconstructed pseudo-temporal trajectories of tumors and McSCs using Monocle. We found three main branches (A-C) (Fig. 6g-i). Branch A composed primarily of McSCs which we anticipate as the initiation site (Fig. 6g, h). Cluster4 was observed at the branching point leading to either clusters 2 & 3 with high neural crest/neuronal (Path1) gene expression or cluster5 with a strong mesenchymal (Path2) signature, suggesting that cluster4 may represent a transition state (Fig. 6i). Branched heatmap (trajectory initiation is at central location) plotted selected gene expression along pseudotime (Fig. 6j). Melanocytic genes was observed along branch A and lost in B and C. At the initiation point of branch B and C, neural crest transcription factors (TFs), such as Sox2, Foxd3, Klf4 and Pou3f1 and signaling molecules were upregulated. Some of these TFs and neuronal genes continue to upregulate along branch B, while in branch C, they were downregulated, followed by the upregulation of EMT transcription factors, such as Twist1 and Snai1 and mesenchymal genes. These data show that the transcriptional heterogeneity inherent to human melanoma, along with many of its features, is found in McSC-derived murine melanoma (Fig. 6k).

Discussion
Despite enormous progress in our understanding of melanoma, the cellular origin of this disease remains poorly understood. Two obstacles have hampered our progress in this area: (1) the lack of specific CreER lines that can unequivocally mark melanocyte stem cells; and (2) a seemingly poor resemblance of mouse melanoma models to human melanoma. In using the c-Kit-CreER model for McSC-specific Braf:Pten activation, we were able to develop a melanoma mouse model resembling some of the key features of human melanoma. In agreement with Moon et al., we observed that tumorigenic McSCs give rise to epidermal melanomas 5 . However, Moon et al. implicated that a second hit, in this case UVB irradiation, was a requisite for formation of epidermal melanoma whereas we have demonstrated that McSCs can give rise to epidermal melanoma during normal onset of hair regeneration ( Supplementary Fig. 14). Further, we have identified the niche factors involved in this process. Wnt and Edn ligands, paracrine niche factors that drive normal melanocyte differentiation during anagen 42,43 , can also fuel the transformation of McSCs into melanoma.
The microenvironment is known to play important roles in tumor initiation, progression and metastasis. Most studies have focused on how niche signals regulate behaviors of tumor cells, and the role of the microenvironment in stem cell to tumor cell transformation remains under-evaluated. Clinical observations revealed that normal donor hematopoietic cells transplanted into leukemia patients can become leukemic, indicating that an abnormal HSC niche could initiate hematopoietic malignancy 61 . This was supported by experiments in which the transplantation of normal hematopoietic cells into retinotic acid receptor gamma (Rary) deficient mice developed myeloproliferative-like disease, whereas transplantation of Rary-deficient hematopoietic cells into wildtype mice did not 62 . A recent study has further demonstrated that carcinogen-driven Bmp2 upregulation by the breast SC niche can induce the transformation of a nonmalignant mammary epithelial cell line into luminal breast tumors 63 . Although oncogenic SCs may also hijack normal niche signals, this has proven difficult to test in vivo in part due to the cellular complexity of the niche. The McSC system is unique in that McSC renewal and Fig. 4 Epidermal melanoma cells invade the dermis and undergo phenotypic switch to a neuronal-like state. a Schematic showing TAM treatment and skin implantation onto nude mice for b-h. b, c Immunofluorescence for GFP at 10 days b and 17 days c after skin implantation in c-Kit-CreER:Braf:Pten:GFP mice. Arrowheads point to GFP + melanoma cells that invade into the dermis. d Immunofluorescence for GFP (green) and E-cadherin, MCAM (red) at 17 days after skin implantation in c-Kit-CreER:Braf:Pten:GFP mice. e Immunofluorescence for Sox10 (red) and E-cadherin, MCAM (green) in human melanoma specimen during the invasive radial growth phase. f, g Immunofluorescence images for GFP (green) and Dct, MITF, Nestin, GFAP, Tubb3 (red) and brightfield images in areas of dermal melanoma at 22 days f and 30 days g after skin implantation in c-Kit-CreER:Braf:Pten:GFP mice. h Dot plot showing percentage of pigment, Dct, Nestin, GFAP and Tubb3 in GFP + dermal melanoma cells (mean ± s.d.; 3 independent tumors from 3 mice were analyzed in each group, n = the number of randomly selected tumor areas analyzed. At least 1500 GFP + tumor cells were analyzed in each group). i Schematic model illustration showing that upon oncogenic induction, McSCs first give rise to epidermal melanoma. Then melanoma cells in the epithelial compartments invade into the dermis. Finally, dermal melanoma cells undergo a phenotypic switch to acquire neuronal signatures. j Schematic showing 4hydroxytamoxifen (4HT-TAM) treatment and analysis regimen. k Immunofluorescence images for Tomato (red) and Sox10, Dct, Tubb3, Nestin (green) and brightfield images in c-Kit-CreER:Braf:Pten:Tomato mouse skin at 47 days. l Immunofluorescence images for Tomato (red), Dct (green) and brightfield images in lymph node of c-Kit-CreER:Braf:Pten:Tomato mouse at 47 days. m Survival graph showing the days of survival of indicated mice (n = the number of mice). Dashed line outlines the boundary of epithelium and dermis. Scale bars, 50 μm. Source data are provided as a Source Data file differentiation is driven by surrounding epithelial SCs, which can be genetically targeted for analysis 42,44,64 . Additionally, unlike most SC systems, the regenerative events of McSCs and their niche cells are synchronized and therefore easier to modulate.
A traditional view in cancer biology is that tumorigenesis recapitulates embryonic development 19   precede the expansion of melanoma 48 . In agreement with this concept, our study showed that transformed McSCs may overcome the restrictions of their adult lineage and re-acquire the plasticity of their embryonic neural crest lineage origin. Nevertheless, our model shows that such a phenotypic change occurs after epidermal melanoma expansion and subsequent invasion into the dermis, instead of being the initial step for transformation. We speculate that the observed difference may be due to differing microenvironment and structural features in zebrafish scales and mammalian skin. Phenotypic switching is a well-known phenomenon in human melanoma and is thought to underlie the phenotypic intra-tumor heterogeneity and resultant drug resistance of an established melanoma 4,17 . Human melanoma heterogeneity has been studied at the transcriptional level in 19 patient-derived melanoma tumors analyzed by single cell RNA-seq 58 . In this work, the intra-tumoral heterogeneity correlated with 1) segregation of proliferative and dormant cells; 2) heterogeneity of Mitf + /Axl − and Axl + /Mitf − populations resistant to MAPK inhibitor therapy; 3) Identification of a cluster of enriched genes related to stress response, therapy resistance and migration. In our single cell analysis of McSC-derived tumors, we observed all 3 conditions, suggesting that our SC-derived tumors phenocopy heterogeneous human melanoma populations. In the clinic, lack of response to therapies frequently correlates with tumor heterogeneity in which cell subsets respond unequally to therapy. We predict that the melanoma tumor in our model is resistant to Braf inhibitors, such as Vemurafenib/Dabrafenib given the fact that most of melanoma cells in our tumor is Axl + /Mitf −58 .
The comparison of our c-Kit-CreER and Tyr-CreER driven melanoma models suggest that dermal melanocytic cells may also be a source of melanoma. The exact distribution, composition and characteristics of dermal melanocytes are not well studied. Yet, it is noteworthy that the presence of dermal stem cells with the ability to regenerate melanocytes have been identified 65,66 , and their melanoma-forming capacity has been postulated [67][68][69] . Future studies identifying specific tool to target dermal melanocytes/stem cells is needed to determine their ability to form melanoma. Another intriguing question is whether fully differentiated melanocytes can produce melanoma 5,6 . A previous important observation suggests that pigmented epidermal melanocytes produce dermal melanoma in the tail skin of Tyr-CreER:Braf:Pten mice. The development of definitive genetic tools to specifically trace differentiated melanocytes during tumorigenic transformation and dermal invasion may effectively validate such observations in the future. To induce Cre recombination, Tamoxifen (TAM) (Sigma-Aldrich) treatment was performed by intraperitoneal (i.p.) injection (0.1 mg/g body weight) of a 20 mg/ml solution in corn oil per day 70 except for one experiment (Fig. 4j-m) in which mice were topically treated with 4-hydroxytamoxifen (4HT-TAM) (Sigma-Aldrich). Mice driven by c-Kit-CreER were injected with TAM for 3 days, mice driven by Tyr-CreER were injected with TAM for 3 or 7 days (see experiment scheme for each experiment), immunodeficient nude mice transplanted with McSCs were injected TAM for 7 days. Topical administration of 4-HT was conducted by preparing a 50 mg/ml solution of 4-HT in dimethylsulfoxide (DMSO). In total 10 μl of the 4-HT solution was topically applied on a defined area of 2 × 2 cm on the depilated back skin daily for 3 days.
Mice that contain tetracycline-inducible transgenes were given doxycycline (1 g/kg; BioServ) containing diet. For BrdU administration, mice were injected daily intraperitoneally during indicated periods with 10 mg ml −1 BrdU (Sigma) solution in PBS at 50 μg/g body weight 54 .
To experimentally induce the synchronized hair follicle cycle, hair depilation was performed during telogen phase (3 weeks or 7 weeks old) to induce anagen phase. We applied hair removing wax (Sally Hansen) to dorsal skin and gently peeled it off or physically plucking hair by hand. Skin biopsies were obtained from back skin of mice at various time points after treatments. Mice are anesthetized by inhalation of isoflurane during both procedures.
Skin grafting. Full-thickness skin grafting was performed as published with slight modifications 71 . Donor mice were injected TAM for 3 days at 3 weeks old. At 1 day after last TAM treatment, both donor and recipient mice were anesthetized by inhalation of isoflurane. A piece of 1.5 cm × 1.5 cm full-thickness back skin was excised with surgical scissors from donor mice and spread onto a sterile gauze soaked with PBS. A piece of 1.5 cm × 1.5 cm skin was excised from the back of recipient mice. The donor skin was then placed and sutured onto the recipient mice by 8 interrupted sutures of 6-0 silk (Ethilon). Skin biopsies were obtained from grafted skin of mice at least 10 days after grafting.
Immunofluorescence. For paraffin sections, tissue (skin, tumor or lymph node) excised from mice was fixed overnight in 4% paraformaldehyde (PFA) at 4°C. After sequential dehydration in increasing concentrations of ethanol and xylene, the tissue was embedded in paraffin. Paraffin sections were cut at 6 μm, deparaffinized and microwaved in 10 mM Tris and 1 mM EDTA (pH 8.0) for antigen retrieval. Tissue sections were then incubated with the primary antibodies listed below for 2 h at room temperature or overnight at 4°C in PBT (PBS + 0.  Immunofluorescence was also performed on formalin-fixed, paraffinembedded human melanoma samples. The samples were obtained from Interdisciplinary Melanoma Cooperative Group (IMCG) biospecimen database of New York University (NYU) Langone Medical Center and evaluated by an attending pathologist as superficial-spreading melanoma. The samples that contain a few dermal invasions (invasive radial growth phase) were used for analysis. Sections were counterstained with 4′,6-diamidine-2′-phenylindole dihydrochloride (DAPI). All immunofluorescence were analyzed and photographed with a Nikon Eclipse Ti microscope. The images were processed using NIS-Elements software and Adobe Photoshop.
X-gal staining. Skin tissue was fixed in 4% PFA for 30 min at 4°C and processed for β-galactosidase activity overnight at 37°C and samples imaged using a Nikon Eclipse Ti microscope.
McSCs isolation. Mice were sacrificed and their fur clipped. Mice were then rinsed in betadine, followed by 70% ethanol. Scalpel blades were used to remove subcutaneous fat and skin was rinsed in PBS and cut into 0.5 cm×0.5 cm pieces, followed by incubation in 0.25% Trypsin for 2 hr at 37°C. Epidermis was separated from the dermis using forceps and scalpel blades and the epidermis was chopped finely and transferred into Media A (DMEM, 10% FBS, 1x penicillin/streptomycin). In total 2.5 × 10 4 melanocytes were combined with the neonatal epidermal (1 × 10 6 cells) and dermal (1 × 10 7 cells) cells from wildtype albino neonates that were prepared according to previous report 72 . The mixture of cells was centrifuged at 200 × g and pellet was re-suspended in 30 µL of Media A and the cellular suspension was injected subcutaneously into nude mice 54 . Nude mice were injected TAM for 7 days immediately after McSC transplantation. The tumors were harvested at 47 days after transplantation.
Tumor preparation for RNA extraction and bulk RNA-seq. Tumors from transplantation assays were collected in Media A and minced finely with a scalpel. Tumors were dissociated by agitation in Media A with 0.35% Collagenase I (Worthington) for 1 h at 37°C, followed by filtering through 100 μm cell strainer to obtain single cell suspension. Tomato + cells were isolated by cell sorting on Beckman Coulter MoFlo cell sorter or Sony SY3200 cell sorter. RNA was extracted from the sorted cells using RNeasy plus micro kit (Qiagen). Purified RNA samples were provided to Genome Technology Center at NYU Langone Medical Center for quality control and sequencing. Paired-end 50 bps raw reads were obtained from Illumina Hiseq-2500 sequencer, demultiplexed using Illumina bcl2fastq conversion software (1.8.4) and mapped to mouse genome (NCBI37/mm9) by using Bowtie aligner (0.12.9) 73 . Mapped reads were assigned to Ensemble gene model (Mus_musculus.NCBIM37.67.gtf) with feature count package 74 . To perform statistical analysis for significant differential expressed genes, edgeR (3.4.2) 75 was used for the data normalization and comparison between experimental and control groups. The p value was adjusted with FDR method. Significant genes were selected as FDR < 0.02 and fold change > 4. The Principal Component Analysis (PCA) plot was generated using DESeq2 (1.22.2) 76 . Quality control data of bulk RNA-seq was shown in Supplementary Table 1.
Single-cell RNA-Seq and data analysis. Tumors were dissociated into single cell suspension as mentioned above. McSC suspension from 3 telogen Tyr-CreER: R26R-Tomato mice and tumor cell suspension from 3 tumors derived from transplantation of McSCs from Tyr-CreER:Braf:Pten:Tomato mice and were used to isolate Tomato + cells by cell sorting on Sony SY3200 cell sorter. The sorted suspensions were mixed into one sample of tumor cell suspension and one sample of McSC suspension.
The sorted cellular suspensions were loaded on a 10x Genomics Chromium instrument to generate single-cell gel beads in emulsion (GEMs). Approximately 5000-10,000 cells were loaded per channel.  77 . Libraries were run on an Illumina HiSeq 4000 as 2 × 150 paired-end reads, one full lane per sample. Sequencing results were demultiplexed and converted to FASTQ format using Illumina bcl2fastq software. The Cell Ranger Single-Cell Software Suite was used to perform sample demultiplexing, barcode processing, and single-cell 3′ gene counting. The cDNA insert was aligned to the mm10/GRCm38 reference genome. Only confidently mapped, non-PCR duplicates with valid barcodes and UMIs were used to generate the gene-barcode matrix. Further analysis and visualization was   (2)(3)(4)(5) showed enrichment in the 3 key signatures observed by bulk RNA-seq (Fig. 5). NES normalized enrichment score performed using Seurat (v2.3.4), an R package containing implementations of commonly used single-cell analytical techniques, including the identification of highly variable genes, dimensionality reduction, standard unsupervised clustering algorithms, and the discovery of differentially expressed genes and markers 78 . The Seurat object was generated from digital gene expression matrices. The parameter of Filtercells is nGene (200-1000) for McSC, nGene (200-5000) for tumor and percentage of mitochondria genes (0 to 0.1). In the standard preprocessing workflow of Seurat, we selected 4614 highly variable genes for following PCA. Then we performed cell cluster and t-SNE. Ten principal components were used in cell cluster with the resolution parameter set at 0.2, resulting in 7 clusters. We excluded two small clusters that are enriched in keratinocytes and macrophages and used the remaining cells to perform variable gene identification (4834 genes), PCA, cell clustering and t-SNE. Nine principal components were used in cell cluster with the resolution parameter set at 0.33. For pseudo-temporal analysis, digital gene expression matrices with annotations from Seurat were analyzed by Monocle v2.10.1 79 . The top 1, 500 genes with highest dispersion (variation/mean) were used to construct the pseudo-temporal trajectory.
Gene Set Enrichment Analysis (GSEA). For the bulk RNA-seq GSEA (Fig. 5f), genes differentially expressed between SC-derived tumor cells and non-transformed McSCs were rank-ordered from high to low based on their fold change. For single cell RNA-seq GSEA (Table 1), genes enriched in each of melanoma cell clusters (clusters [2][3][4][5] in comparison compared to McSCs (cluster 1) were rank-ordered based on fold change. For single cell RNA-seq GSEA (Fig. 6c-e), indicated melanoma clusters were compared to each other and their differentially expressed genes were rank-ordered based on fold change. We then queried these pre-ranked gene lists for their enrichment in 4 annotated gene sets acquired from The Molecular Signature Database (MSigDB): GO_CELL_MORPGOGENESIS_INVOLVED_IN_NEURON_DIFFER-ENTIATION, KEGG_CELL_CYCLE, KEGG_ECM_RECEPTOR_INTERACTION, and LEE_NEURAL_CREST_STEM_CELL_UP 57 , using the preranked GSEA analysis tool 80 . FDR q-value < 0.25 was deemed significant.
Quantification and statistical analyses. The measurement of quantifications can be found in y-axis of dot plots in the figure and in figure legends. The statistical details of each dot plot can be found in the figure (n number, P value). The exact meaning of n number is described in the corresponding figure legend. Pairwise comparisons between two groups were performed by two-sided unpaired statistical analysis using Student's t test. Statistical significances were considered significant if P < 0.05. The exact P value is labeled in the figure. Experimental data are demonstrated as the mean ± standard deviation (s.d.). Statistical analysis and plotting was done using Microsoft Excel and Graphpad Prism.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All RNA-seq data reported in this paper are deposited in NCBI Gene Expression Omnibus(GEO) database. The accession numbers are GSE96966 for bulk RNA-seq and GSE113502 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi] for single cell RNA-seq. The source data underlying all quantifications in Figs. 1-5 and Supplementary Figs. 1-5, 8 and 12 are provided as a Source Data file. All the other data supporting the findings of this study are available within the article and its supplementary information files and from the corresponding author upon reasonable request. A reporting summary for this article is available as a Supplementary Information file.