Tfap2b specifies an embryonic melanocyte stem cell that retains adult multifate potential.

Melanocytes, the pigment-producing cells, are replenished from multiple stem cell niches in adult tissue. Although pigmentation traits are known risk factors for melanoma, we know little about melanocyte stem cell (McSC) populations other than hair follicle McSCs and lack key lineage markers with which to identify McSCs and study their function. Here we find that Tfap2b and a select set of target genes specify an McSC population at the dorsal root ganglia in zebrafish. Functionally, Tfap2b is required for only a few late-stage embryonic melanocytes, and is essential for McSC-dependent melanocyte regeneration. Fate mapping data reveal that tfap2b+ McSCs have multifate potential, and are the cells of origin for large patches of adult melanocytes, two other pigment cell types (iridophores and xanthophores), and nerve-associated cells. Hence, Tfap2b confers McSC identity in early development, distinguishing McSCs from other neural crest and pigment cell lineages, and retains multifate potential in the adult zebrafish.

Introduction

Melanocytes produce a black-brown pigment that gives color and patterns to the animal kingdom and protect from harmful UV irradiation (Mort et al., 2015; Nguyen and Fisher, 2019). How melanocytes develop in embryogenesis and maintain a regenerative melanocyte population in adulthood is central to our understanding of pigmentation and pattern formation through evolution as well as in human disease (Baxter and Pavan, 2013; Irion and Nusslein-Volhard, 2019; Kelsh and Barsh, 2011; Kinsler and Larue, 2018; Nusslein-Volhard and Singh, 2017; Parichy and Spiewak, 2015). In melanoma, a deadly cancer of the melanocyte, neural crest (NC) and melanocyte developmental mechanisms become reactivated in pathological processes driving initiation, metastasis, survival, and drug resistance (Diener and Sommer, 2020; Johansson et al., 2020; Kaufman et al., 2016; Konieczkowski et al., 2014; Marie et al., 2020; Rambow et al., 2018; Shakhova et al., 2012; Travnickova et al., 2019; Varum et al., 2019; White et al., 2011). Therefore, although mechanisms that underpin NC and melanocyte development are important for understanding fundamental biology, in the melanoma context, understanding dysregulation and heterogeneity of developmental lineages provides a rich source of drug targets for the next generation of therapies (Patton et al., 2021).

Zebrafish are uniquely amenable to advanced in vivo imaging of NC and pigment cell lineages and have emerged as a powerful system to study the developmental lineages in pattern formation and melanoma (Kelsh et al., 1996; Owen et al., 2020; Travnickova and Patton, 2021). In zebrafish, embryonic trunk melanocytes are directly derived from the NC, but a rare melanocyte subset in the lateral stripe originates from melanocyte stem cells (McSCs; also referred to as MSCs in zebrafish), a postembryonic progenitor population dependent on ErbB-kinase signaling (the epidermal growth factor family of receptor tyrosine kinases). McSCs are essential for replenishing melanocytes following regeneration or during growth, and ablation or genetic loss of McSCs leads to large patches devoid of melanocytes in the adult (Budi et al., 2008, 2011; Dooley et al., 2013; Hultman et al., 2009; Hultman and Johnson, 2010; Taylor et al., 2011). Imaging in zebrafish over time shows that McSCs reside in a niche at the site of the dorsal root ganglia (DRGs) and give rise to melanoblasts that migrate along nerves to the epidermis, where they differentiate into pigmented melanocytes (Budi et al., 2011; Dooley et al., 2013; Singh et al., 2014, 2016).

Although fate mapping directly from McSCs has not yet been performed (because of a lack of specific lineage markers), lineage tracing from the NC transcription factor Sox10 indicated that the adult zebrafish pattern is governed by progenitors called MIX cells, an NC subpopulation that is multipotent for all three pigment cell types—black-brown melanocytes, silver iridophores and yellow xanthophores—and gives rise to neurons and glia (Singh et al., 2016). Adult melanocytes and iridophores are primarily derived from MIX cells, whereas adult xanthophores have a dual origin and originate from MIX cells and an embryonic population that expands at the onset of morphogenesis (McMenamin et al., 2014; Singh et al., 2016). Although not directly identified, this analysis indicated that MIX cells persist beyond embryonic and larval stages and are the source of almost all melanocytes that pattern the adult fish (Singh et al., 2016). Although it is clear that these MIX cells are derived from the NC in zebrafish development, their molecular signature, how they become specified, and how they relate to McSCs remains unknown.

Here, using single-cell RNA sequencing (scRNA-seq), imaging, and lineage tracing, we discover that the transcription factor Tfap2b and a select set of its target genes specify the ErbB-dependent multipotent McSC population. Although a mutation in its paralog tfap2a is associated with premature hair graying in humans (Praetorius et al., 2013) and results in pigmentation defects in mice and zebrafish (Seberg et al., 2017), the role of tfap2b in melanocyte biology is largely unexplored. In melanoma, TFAP2B marks a critical node in residual disease (Rambow et al., 2018), although its precise function has yet to be defined. We show that tfap2b+ cells reside at the McSC niche and are progenitors for all three pigment cell types and some nerve-associated cells. We conclude that a Tfap2b transcriptional program specifies a subset of ErbB-dependent cells within a larger embryonic MIX+ cell population to become McSCs at the DRGs. These cells are responsible for generating melanocytes and other pigment cells in adult zebrafish, bridging the gap between embryonic and adult patterning.

Results

McSCs maintain an NC identity at the niche

We performed live imaging of zebrafish embryos to visualize McSCs as they migrate into the DRG niche. We used the transgenic reporter lines Tg(mitfa:GFP), which marks melanoblasts, and Tg(nbt:dsRed), which marks the neural tube and peripheral nerves (Dooley et al., 2013). Prior to melanization, during the first 22–27 h post fertilization (hpf), we observed a highly migratory GFP+ subpopulation that traversed the McSC niche and migrated along nerves to the skin (Figures 1A and 1B). We consider these to be melanocytes that develop directly from the NC and use nerves to migrate before generating the embryonic stripes. Subsequently, we observed a distinct, smaller subpopulation of round melanoblasts that was maintained at the niche and a population of stationary columnar cells that coated the nerves (Figures 1A and 1C). These cells have been described previously as McSCs and progenitors (Budi et al., 2011; Dooley et al., 2013; Johansson et al., 2020; Kelsh and Barsh, 2011; Singh et al., 2014, 2016). We confirmed that the small, round mitfa cells were McSCs by treatment with an ErbB inhibitor, which eliminated this population at the niche site as well as stationary nerve-associated progenitors (Figures 1D and 1E; Budi et al., 2008, 2011; Dooley et al., 2013; Hultman et al., 2009).

Figure 1

McSCs maintain an NC identity at the niche

(A–C) Schematic of developing McSCs and melanocytes in the zebrafish embryo. mitfa-expressing melanoblasts (Mbs; green) develop directly from the NC and travel to the skin dorsolaterally (data not shown) or ventrally along the neural tube (NT) and peripheral nerves (red; B). A subset of those cells establishes at the site of the of the perspective dorsal root ganglia (DRGs) and become McSCs (blue; C). McSC establishment is sensitive to ErbB-kinase inhibitors (ErbBi).

(B and C) Melanoblasts migrating along the axons (B) and a newly established McSC at the site of a perspective DRG (C). Shown are confocal stacks (20 μm) of Tg(mitfa:GFP; nbt:dsRed) embryos imaged laterally at 22 hpf. Scale bars, 20 μm.

(D) ErbB kinase activity is required for McSC establishment at the niche (white arrow). Shown are confocal stacks (20 μm) of Tg(mitfa:GFP; nbt:dsRed) embryos treated with DMSO or ErbBi. Standard deviation (STD) projection. Scale bars, 20 μm.

(E) Quantification of McSC niche occupancy. Tukey honestly significant difference (HSD) test; ∗∗∗p < 0.0001 (3 replicates, 5 embryos/condition/replicates). Lines in boxes indicate the medians, and whiskers indicate data within 1.5 interquartile range of the upper and lower quartiles.

(F) McSCs maintain NC identity at the niche. Shown are confocal stacks (30 μm) of Tg(crestin:mCherry; mitfa:GFP) embryos treated with DMSO or ErbBi. McSCs (white arrows) and nerve-associated cells (yellow arrows) are dependent on ErbB kinase. STD projection. Scale bars, 20 μm.

(G) McSCs and nerve-associated precursors express mitfa:GFP and crestin:mCherry, but expression is lost by 120 hpf. Shown are confocal stacks (30 μm) of Tg(crestin:mCherry; mitfa:GFP) embryos. The lower edge of the NT is indicated (white dotted line) on fluorescence and corresponding bright-field images (top). Bright-field, average intensity (AVG) projection; fluorescence, STD projection. Scale bars, 50 μm.

Tg(mitfa:GFP) is expressed by McSCs, but it is not an exclusive McSC marker, nor does mitfa mutation or knockdown prevent acquisition of McSC identity or establishment at the niche (Budi et al., 2011; Dooley et al., 2013; Johnson et al., 2011), indicating that unknown mechanisms confer McSC identity from the NC. To resolve such mechanisms and identify specific McSC markers, we generated a Tg(crestin:mCherry; mitfa:GFP) double-transgenic line that enabled visualization of the mitfa-expressing melanocyte lineage (GFP+) as well as its origin in the NC (mCherry+) (Kaufman et al., 2016) during McSC establishment. We observed newly established McSCs that retained crestin expression (GFP+ mCherry+) and cell columns along the peripheral nerves that were also dependent on ErbB signaling (Figure 1F). Notably, the expression levels for GFP and mCherry fluorescence were heterogeneous for both cell populations (Figures 1F and 1G). By 48 hpf, GFP and mCherry fluorescence was restricted specifically to the McSCs, whereas differentiating melanoblasts maintained GFP+ fluorescence (Figure 1G). In McSCs, expression of both transgenes was maintained through organogenesis before eventually subsiding by 120 hpf (Figure 1G). Thus, unlike differentiating skin melanoblasts, McSCs maintain an NC identity during the establishment and specification phase.

scRNA-seq identifies six distinct NC pigment cell lineage populations

To identify a molecular signature specific for McSCs, we needed to first understand how pigment cell populations arise from the NC. To this end, we performed droplet-based scRNA-seq on GFP+ and/or mCherry+ cells sorted from the Tg(crestin:mCherry; mitfa:GFP) transgenic line using the 10x Chromium system (Figures 2A and S1A), sequencing 1,022 cellular transcriptomes, 996 of which passed our quality control (Figures S1B–S1J). When we visualized the data in two-dimensional space by applying the Louvain clustering algorithm (Butler et al., 2018) and uniform manifold approximation and projection (UMAP) (McInnes et al., 2018; Figure 2B), we found that the number of expressed genes was distributed uniformly across the cells (Figure S1E; Table S1). To assign cluster identities, we employed a combination of known markers and projections to previously published datasets (Farnsworth et al., 2020; Kiselev et al., 2018; Saunders et al., 2019; Wagner et al., 2018; Figures S2C–S2F; Table S2).

Figure 2

scRNA-seq identifies six distinct NC pigment cell lineage populations

(A) Schematic of the experimental protocol. GFP+ and/or mCherry+ cells were isolated from Tg(crestin:mCherry; mitfa:GFP) embryos at 24 hpf.

(B) UMAP of GFP+, mCherry+, and double-positive cells (n = 996 cells) obtained after Louvain clustering (dimensions, dims = 20, resolution = 0.5).

(C) Heatmap showing the average log2 fold change expression of five selected genes per cluster identified in (B). The average log2 fold change expression across the 6 clusters of sox10, crestin, mitfa, mCherry, and GFP expression levels are also presented for comparison.

(D) UMAP representations of (B) with color change from blue (negative) to red (positive) based on log2 mRNA expression of pcna, twist1a, foxd3, aox5, tfec, mitfa, and tyrp1b.

(E) Pseudotime ordering of the cells in (B). Top: pseudocoloring based on pseudotime scores. Bottom: pseudocoloring based on cluster identity.

(F) RNA velocity analysis of the UMAP represented in (B) (top) and simplified representation (bottom).

See also Figures S1–S3 and Tables S1 and S2.

We identified two proliferative NC populations expressing classical markers. One was characterized by the almost exclusive expression of the basic-helix-loop-helix (bHLH) transcription factor and epithelial-to-mesenchimal transition (EMT) gene twist1a (and other twist genes; possibly cranial NC cells) and the second by expression of foxd3, a known Wnt-regulated NC gene (Figures 2B–2D; Table S2). In addition, we identified four distinct clusters that expressed a combination of different chromatophore markers (Figures 2B–2D). Although expression of chromatophore marker genes is indicative of developmental potential, they are not predictive of cell fate; for example, mitfa is expressed in cells that differentiate into xanothophores (Parichy et al., 2000). However, our scRNA-seq dataset captured cells that are consistent with hypothesized cell populations described through sox10+ lineage tracing (Singh et al., 2016), gene network studies (Petratou et al., 2018, 2021), and early imaging studies (Bagnara et al., 1979). One cluster expressed NC genes, id genes, and chromatophore markers from all three pigment cell types, including melanocyte (mitfa, dct, and tyrp1b), iridophore (pnp4a), and xanthoblast (aox5) markers. For simplicity, we refer to this as a MIX+ cluster because it is positive for expression of markers from all three pigment cells: melanocytes, iridophores, and xanthophores. We were intrigued to see that the cells belonging to this cluster expressed less than half of the genes overall. This molecular phenotype is congruent with a stem cell identity because self-renewing hematopoietic stem cells at the top of the differentiation hierarchy express the lowest number of genes and total mRNA, with expression gradually increasing in differentiated cells (Nestorowa et al., 2016; Figure S2B). Cells in the immediately adjacent clusters expressed relatively high levels of melanoblast markers concomitant with markers for two different chromatophores; these included cells that expressed mitfa and aox5 (MX+ cells) and cells that express mitfa and pnp4a (MI+ cells) (Figures 2C and 2D). Melanoblasts were very similar to MI+ cells in gene expression profile but had higher expression of melanocyte differentiation genes. We suggest that the melanoblasts are the skin-associated mitfa cells that start to melanize by 28 hpf and are responsible for the embryonic melanocyte pattern. We did not find evidence of cells that express iridophore and xanthophore markers without mitfa (IX+ cells).

To understand lineage relationships between clusters, we performed a pseudotime analysis (Figure 2E). We found the cells to be part of a developmental continuum that originates in the NC and transitions through a MIX+ stage before differentiating into MI+ cells and MX+ cells or melanoblasts, consistent with a common origin for pigment cells (Bagnara et al., 1979; Petratou et al., 2018, 2021). A branchpoint emerged between melanoblasts and the MI+ and MX+ cells, indicative of two distinct melanoblast populations. This is consistent with our imaging (Figure 1), which shows mitfa melanoblast populations in the skin (before these become pigmented) and lining the nerves (relatively undifferentiated at that stage). Through RNA velocity analysis (La Manno et al., 2018), the developmental lineage relationships were found to be consistent with the pseudotime ordering (Figures 2F and S3). Importantly, despite similar transcriptomes (Figure 2D), our analysis indicate that the MI+ cells are not progenitors differentiating into the embryonic melanoblasts (black dots, Figures 2E and 2F), suggesting two distinct routes of melanoblast formation in zebrafish.

Identification of ErbB-dependent McSCs by scRNA-seq

Given that the MIX+ cells were enriched for mitfa crestin cells (Figures 2C and 2D), we predicted that they could function as a possible source for McSCs. Because McSCs are ErbB kinase dependent, we designed a scRNA-seq experiment in ErbB kinase inhibitor (ErbBi)-treated Tg(crestin:mCherry; mitfa:GFP) embryos (Figures 3A and 3B). Clustering 347 cells derived from ErbBi-treated embryos revealed a loss of MI+ cells compared with our untreated embryos. Further, relative to untreated embryos, the melanoblast population in ErbBi-treated embryos had reduced differentiation markers, and we therefore called these cells early melanoblasts (Figures 3B and S4A). Pseudotime analysis showed that cells from the double-transgenic embryos form a developmental continuum from NC through MIX+ states to MX+ and early melanoblasts (Figures 3C and S4). Interestingly, in contrast to a loss of MI+ cells, twist1a cells were proportionally increased in ErbBi-treated embryos, and the early melanoblasts and twist1a cells were enriched for mesenchymal and proliferative markers (Figures 3D, 3E, S4B, and S4C).

Figure 3

Identification of ErbB-dependent McSCs by scRNA-seq

(A) Schematic of the scRNA-seq experimental protocol for ErbBi-treated zebrafish embryos (24 hpf). ErbBi treatment: 4–24 hpf.

(B) UMAP of GFP+, mCherry+, and GFP+ mCherry+ cells (n = 346 cells) from ErbBi-treated embryos after Louvain clustering (dims = 10, resolution = 0.5).

(C) Pseudotime ordering of the cells in (B). Left: pseudocoloring based on pseudotime scores. Right: pseudocoloring based on cluster identity.

(D) UMAP of GFP+, mCherry+, and GFP+ mCherry+ cells (n = 1,343 cells) from untreated and ErbBi-treated embryos after Louvain clustering (dims = 12, resolution = 1).

(E) UMAP in (D) pseudocolored with the cell origin. Dashed lines highlight clusters enriched with ErbBi-treated cells.

(F) Pseudotime ordering of the cells in (D). Cell states present in the untreated embryos (dashed box) are absent in ErbBi-treated embryos. Top left: pseudocoloring based on pseudotime scores. Bottom left: pseudocoloring based on cluster identity. Right: split views (by treatment). Cell states and their inferred position in the 24-hpf embryo are also indicated.

(G) ErbB kinase-dependent McSCs (red) and MI+ cells (brown) are highlighted on UMAP presented in Figure 2B.

(H) Minimum spanning tree presented in Figure 2E pseudocolored according to the cell states described in (F). The inferred position and the McSC branchpoint are indicated.

See also Figures S1 and S4 and Tables S1 and S2.

We conclude that the ErbB-dependent and nerve-associated progenitors represent MI+ cells (Figures 3C–3E and S4D). However, from our UMAP, we could not easily distinguish a specific ErbB-dependent MIX+ subpopulation that represents McSCs (Figure 3E). To understand the hierarchical dependence of the pigment cells emerging from the NC, we performed an integrated pseudotime analysis between untreated and ErbBi-treated embryos (Figure 3F). In untreated embryos, cells were distributed along a continuum from the early NC to the pigment lineage, with early pigment progenitors as a “pass-through” state. The pigment lineage split into branches ending in defined states for directly developing melanoblasts, axon-associated progenitors, or MX+ cells. When comparing pseudotime ordering of cells between datasets, we observed that melanoblast-derived lineages were reduced overall and that axon-associated progenitors, robustly present in untreated embryos, were almost entirely lost upon ErbBi treatment. In addition, cells obtained from ErbBi-treated embryos were largely deficient for transitioning cell states but enriched for twist1a NCCs (Figure 3F), possibly providing an explanation for the reduced number of foxd3 NCCs in the UMAP. Based on these data, we hypothesize that ErbB-dependent cells capable of transitioning between MIX+ and differentiating pigment cell lineages are McSCs. When we mapped these cells back to the untreated dataset, we found that the McSCs formed a subpopulation within the MIX+ cluster (Figure 3G) and that the McSCs, in fact, represent a branchpoint in the pseudotime analysis, supporting our hypothesis that their identity is a distinct cell state within the MIX+ cell population (Figure 3H).

McSC identity is specified by a Tfap2b transcriptional program

Next we performed a differential expression analysis comparing the transcriptomes of the putative McSCs and the combined pigment cells derived from untreated embryos (Figures 4A and 4B). As anticipated for a MIX+ subgroup, McSCs expressed all pigment progenitor markers (Figure 4A; Table S3). Combined pigment lineages were enriched for expression of pigment synthesis genes, whereas McSCs were enriched for expression of ribosome biogenesis and splicing genes, similar to what has been reported for other stem cell systems (Brombin et al., 2015; Gabut et al., 2020; Gupta and Santoro, 2020; Recher et al., 2013; Table S4). Moreover, McSCs were enriched in genes associated with neurological disabilities in humans compared with the early pigment progenitors (Table S5).

Figure 4

McSC identity is specified by a Tfap2b transcriptional program

(A) Violin plots of pigment progenitor and McSC gene expression levels. McSCs differentially express a subset of genes (McSC genes) and share expression with cells of the pigment cell lineage (early pigment progenitors, MI+, MX+, and Mbs from untreated embryos in Figure 3F).

(B) Rank plot of differential expression analysis between McSCs and all other states of the pigment lineage. The top differentially expressed gene is tfap2b (log2 fold change = 5.63; adjusted p = 4.99e−5).

(C) Clustered heatmap showing the average expression of 36 Tfap2b targets (orthologs of the targets found in chick by Ling and Sauka-Spengler, 2019) enriched in McSCs and 3 non-differentially expressed targets (sox10, mitfa, and tfap2a).

(D) Rank plot of differential expression analysis between MIX+ cells from untreated embryos and ErbBi-treated embryos. Most of the McSC-Tfap2b target genes are depleted from MIX+ cells in treated embryos.

See also Figure S5 and Tables S3, S4, S5, S6, and S7.

Strikingly, we found transcription factor AP-2 beta (tfap2b) to be the top gene specifically enriched within the transcriptome of the putative McSCs (Figures 4A and 4B). Tfap2b is a transcription factor, functionally redundant with Tfap2a, known to regulate NC and melanocyte development in the mouse and postulated to differentially specify melanocyte precursors together with Mitfa in zebrafish (Chong-Morrison and Sauka-Spengler, 2021; Lignell et al., 2017; Rothstein and Simoes-Costa, 2020; Seberg et al., 2017). Although Tfap2a/e are important in melanoma metastasis (Campbell et al., 2021), Tfap2b was identified as a marker of a subpopulation present in residual disease states following BRAF plus MEK inhibitor treatment (Rambow et al., 2018) and is also a marker of residual disease in our zebrafish melanoma models (Travnickova et al., 2019; Figure S5A). To investigate whether Tfap2b plays a functional role in regulating the McSCs, we reproduced the Tfap2b (biotin) chromatin immunoprecipitation sequencing (ChIP-seq) analysis carried out in chicken NC (Ling and Sauka-Spengler, 2019). After mapping the chicken Tfap2b targets to their zebrafish homologs, we examined the expression of these genes in McSCs as well as in the derivative cell populations within our dataset. We found that a select subset of Tfap2b targets are selectively highly expressed in McSCs compared with other cell states (Figure 4C; Table S6). Critically, these Tfap2b+ McSC target genes are enriched in MIX+ cells derived from control embryos compared with MIX+ cells derived from treated embryos (Figure 4D; Table S7). These data indicate that Tfap2b plays a pivotal role in specification of McSCs through activation of a select set of target genes.

Functional requirement for tfap2b in melanocyte regeneration

Next we wanted to determine whether tfap2b is functionally important for McSC-derived melanocytes but not for NC-derived melanocytes. Thus, we used a melanocyte regeneration assay based on a temperature-sensitive mitfa mutation (mitfa) in which embryos are incapable of generating embryonic melanocytes at higher temperatures but capable of regenerating melanocytes from McSCs at lower temperatures (Johnson et al., 2011; Zeng et al., 2015). We found that morpholino (MO)-mediated knockdown of tfap2b did not affect the development of most NC-derived melanocytes but significantly reduced McSC-derived melanocytes in regeneration (Figure 5A). These results were confirmed with a second, splicing-site MO knockdown (Figure 5A). McSCs contribute to a small population of melanocytes in the embryonic lateral stripe and McSC activity can be assessed through this as a second independent assay in zebrafish embryos. Indeed, in tfap2b knockdown embryos, we found McSC-derived late-stage lateral stripe melanocytes to be reduced, confirming that tfap2b is required for McSC-derived melanocytes even in non-regenerating embryos (Figure 5B). We wanted to assess the effect of tfap2b MO knockdown on adult pigmentation; however, although morphants appeared to be healthy within the first 5 days of development, we were unable to raise the fish beyond 15 days of development, likely because of a potential role of tfap2b cells in other lineages (e.g. neural and craniofacial cells).

Figure 5

tfap2b is expressed at the McSC niche and required for regeneration

(A) tfap2b is required for melanocyte regeneration from the McSC. Shown are images of zebrafish embryos and melanocyte quantification following knockdown of tfap2b in a mitfa regeneration assay. Tukey HSD test, ∗∗∗p < 0.0001 (3 replicates, 20 embryos/condition/replicate). Lines in boxes indicate the medians, and whiskers indicate data within 1.5 interquartile range of the upper and lower quartiles. Scale bars, 200 μm. N.I., not injected; Co. MO, control MO; AUG MO, AUG-directed MO; Splic. MO, splicing MO.

(B) tfap2b is required for late-stage melanocytes from the McSC. Shown are Images of zebrafish embryos and melanocyte quantification following knockdown of tfap2b. Only McSC-derived late-developing lateral stripe melanocytes are reduced in tfap2b knockdown embryos. Arrows highlight missing lateral stripe melanocytes. Tukey HSD test; ∗∗∗p < 0.0001 (3 replicates, 20 embryos/condition/replicate). Lines in boxes indicate the medians, and whiskers indicate data within 1.5 interquartile range of the upper and lower quartiles. Scale bars, 200 μm.

(C) tfap2b:eGFP expression in the McSC. Shown is a merged image of a double-transgenic Tg(tfap2b:eGFP; crestin:mCherry) zebrafish (left) and separated channel images (bright-field, GFP, and mCherry channel). White arrows indicate GFP+/mCherry+ McSCs at the DRGs. Scale bars, 50 μm.

(D and E) tfap2b McSCs require ErbB kinase at the niche. Shown are Tg(tfap2b:eGFP;crestin:mCherry) embryos at 24 hpf (D) and 48 hpf (E), untreated or treated with ErbBi. White arrows indicate the McSC niche. Tukey HSD test; ∗p = 0.0172, ∗∗∗p < 0.0001 (3 replicates, 5 embryos/condition/replicate). Lines in boxes indicate the medians, and whiskers indicate data within 1.5 interquartile range of the upper and lower quartiles. Confocal stacks, 30 μm; STD projection. Scale bars, 50 μm.

See also Figure S5.

McSCs are specified by tfap2b at the niche

We next sought to visualize tfap2b McSCs during development. To this end, we isolated a 1kb tfap2b promoter region and cloned it upstream of the GFP coding sequence to generate a transgenic reporter line, Tg(tfap2b:eGFP), and then crossed this with Tg(crestin:mCherry) animals. We identified cells co-expressing GFP and mCherry at the McSC niche and along nerves (Figure 5C). Critically, these cells were absent upon ErbBi treatment, indicating that Tg(tfap2b:eGFP) expression marks the McSC lineage (Figures 5D and 5E).

In extended imaging analysis, we observed that Tg(tfap2b:eGFP) expression peaked during the first 48 hpf but then declined with crestin:mCherry expression. We saw no overt evidence of tfap2b:eGFP expression in adult zebrafish pigment cells, although expression was sustained in the spinal cord, in accordance with a neuronal expression pattern (Knight et al., 2005). We did not find evidence of tfap2b:eGFP expression in other embryonic melanocytes or pigment cells by live imaging, and analysis of our and other datasets confirmed that tfap2b expression is observed principally in the McSC population in the pigment cell lineage (Figures S5B and S5C). These results indicate that Tfap2b functions transiently to establish McSC identity during embryogenesis and is required during the first few days of embryonic development.

tfap2b McSCs have multifate potential for all adult pigment cell lineages

The highly specific tfap2b gene expression in McSCs (Figures S5B and S5C) prompted us to follow the fate of tfap2b cells through development and into the adult pigment pattern. To this end, we cloned the 1 kb promoter upstream of cre (tfap2b:cre) and injected it into the ubi:switch transgenic line so that fish were mosaic for tfap2b:cre integration, facilitating lineage tracing of cells “switched” from GFP+ to mCherry+ (Mosimann et al., 2011; Figure 6A). Consistent with tfap2b neural expression during development (Knight et al., 2005), we found mCherry+ neurons in the dorsal neural tube at 48 hpf (Figures 6B and S5C). Critically, by 6 days post fertilization (dpf), we found mCherry+ McSCs at the niche and in a few melanocytes in the skin. By 13 dpf, we could clearly detect mCherry+ McSCs and nerve-associated cells as well as a few late-stage embryonic melanocytes in the lateral stripe and xanthophores (Figures 6C and 6D).

Figure 6

Individual tfap2b McSCs during development and at the onset of metamorphosis

(A) Experimental protocol overview for tfap2b mosaic lineage tracing. pEXP-GC2Tol2-tfap2b:cre and Tol2 mRNAs were injected into ubi:switch embryos at the zygote stage and imaged at different stages to identify mCherry+ McSCs and progenitors.

(B–D) tfap2b lineage tracing during zebrafish development at (B) 48 hpf, (C) 6 dpf, and (D) 13 dpf (standard length [SL], 4.2 mm). Shown are representative images of 7–8 animals imaged per stage (1 technical replicate). White pseudocoloring is used for mCherry. Maximum intensity (MAX) projection. Scale bars, 50 μm.

(E) Example of a single McSC followed through development of a single fish and imaged at 48 hpf (left), 120 hpf (center), and 10 dpf (SL, 4.1 mm; right). The McSC produces progenitors that populate the ventrally located axonal projection and a lateral line-associated melanocyte. The embryo was injected with 1.5 pg/nL of tfap2b:cre. Scale bars, 50 μm.

(F) Dramatic expansion of progenitors from the McSCs at the onset of the metamorphosis. Shown is a confocal imaging of a clone expanding from a single McSC at 15 dpf (SL, 4.9 mm). Arrows indicate axon-associated cells. White pseudocoloring is used for the GFP channel and red for the mCherry channel in the left panel, and white pseudocoloring is used for the mCherry channel in the right panel. The embryo was injected with 1.5 pg/nL of tfap2b:cre. Scale bars, 20 μm. I, iridophore, M, melanocyte; Mu, muscle (off-target signal); McSC, melanocyte stem cell; N, NT/spinal cord neuron; X, xanthophore.

See also Figure S5.

To follow the fate of tfap2b McSCs over time, we needed to be confident that clones originated from a single or a few recombination events. To this end, we injected limiting dilutions of tfap2b:cre (1.5 and 3.125 pg/nL) and analyzed 70 embryos at 2 dpf. 54 embryos expressed no mCherry at the McSC but only in other cells (e.g. the neural tube), 13 expressed mCherry at a single visible McSC, and three embryos expressed mCherry at a few McSCs (with a lower level of certainty because of overlap with neuronal cells by epifluorescence microscopy). This indicates that these experiments are within the limiting transgene range and that mCherry+ cells or clones at later stages of development are likely the result of a single or a few recombination events in the McSCs. We followed these fish over 15 days and could clearly visualize cells emerging from the McSC along nerves. In Figure 6E, we show an example of an McSC giving rise to cells along nerves and a late-stage lateral line melanocyte by 5 dpf. We also detected an expansion of mCherry+ cells emerging from the McSC site at the onset of metamorphosis (15 dpf), including a newly differentiating melanocyte (Figure 6F).

To understand the fate of tfap2b cells in the adult pigmentation pattern, we analyzed 56 zebrafish at 1 month of age by confocal imaging. By 1 month, mCherry+ melanocytes, xanthophores, and iridophores were particularly prominent in the developing stripes and occasionally in nerve-associated cells (Figures 7A and 7B). We identified MIX clones that spanned the whole dorsoventral axis in hemisegments and were clearly visible in the scales. We observed one fish with a sparse clone made up of only a few iridophores, xanthophores, and melanocytes distributed along the length of the hemisegment and not adjacent to each other (Figure 7C). This may indicate that not all McSCs have equal potential to generate progenitors or may represent a late-stage McSC activation event.

Figure 7

tfap2b McSCs have multifate potential for all adult pigment cell lineages

(A–C) MIX clones in tfap2b lineage tracing analysis in the adult pigment pattern. Shown are the (A) caudal trunk, (B) tail region, and (C) medial trunk. Magnified images of melanocytes are presented in stripes 2D (AI and AII) and 1D (AIII). Orthogonal projections of a melanocyte (AIV) and a xanthophore (AV) show that cells are localized in the scale. Also shown are (B) a MIX clone spanning the dorsoventral axis and (C) a sparse MIX clone, with magnifications of iridophores, xanthophores, and a melanocyte (top to bottom). White pseudocoloring is used for the GFP channel, and red is used for the mCherry channel in (A)–(C). white pseudocoloring is used for the mCherry channel in the magnified panels and AI–AIII and in (C) magnified images. Shown are representative images of more than 20 fish injected with 25 pg/nL of tfap2b:cre. MAX projection. Scale bars, 50 μm.

(D) Frequencies of the different derivatives in clearly defined clones in 12 juvenile ubi:switch zebrafish injected with low doses of the tfap2b:cre plasmid (1.5 pg/nL, 3.25 pg/nL, and 6 pg/nL). Fish SL in millimeters: (1) 9.7, (2) 8.7, (3) 7.1, (4) 10.5, (5) 9.3, (6) 11.6, (7) 11.1, (8) 10.6, (9) 12.4, (10) 11.2, (11) 13.7, and (12)15.5. Clones were named according to their pigment cell composition. Numbers in the black dots represent mCherry+ melanocytes within the clone.

(E−G) Representative images of the clones analyzed in (D). Shown are (E) a MIX clone with a clone-associated nerve (caudal trunk; the arrow indicates clone extension into the anal fin), (F) an MX clone (medial trunk), and (G) an X clone (rostral trunk). MAX projection. Scale bars, 50 μm. All fish were injected with 1.5 pg/nL of tfap2b:cre. I.s., interstripe; St., stripe.

See also Figure S6.

Our experimental cohort of zebrafish included fish injected with a range of tfap2b:cre concentrations (1.5–25 pg/nL). We observed the same patterns of labeled cells in fish injected at all concentrations of the transgene; however, some very large clones spanned one or more segments along the rostrocaudal axis, indicating that the Cre recombination might have occurred in multiple consecutive McSCs. Therefore, we quantitated the clones derived from 12 fish at 1 month of age that had been injected with low concentrations of tfap2b:cre (1.5 pg/nL, 3.25 pg/nL, and 6 pg/nL) and had clearly definable boundaries. Clone numbers ranged from 0–5 clones/fish but, on average, there were 2 clones per animal. Even at this low tfap2b:cre concentration, many fish had large, clearly visible MIX clones that filled the hemisegment (Figure 7D). Some MIX clones also contained nerves (as seen by (Singh et al., 2016)) and even expanded into the fin regions (Figure 7E). Melanocyte numbers per clone varied from 1 to 60. We also observed large clones composed of just melanocytes and xanthophores (MX clones; Figure 7F) as well as a number of smaller xanthophore-only clones (X clones) within the interstripes (Figure 7G). Again, these data suggest that McSCs may not have equal potential to generate all derivatives. Single X clones could also arise from a tfap2b McSC-derived xanthophore (as seen in Figure 6D at 6 dpf) or, alternatively, from a separate, still unknown tfap2b xanthoblast progenitor at later developmental stages. Finally, we performed cryo-sectioning to analyze tfap2b-derived cells in juvenile zebrafish and found mCherry+ cells migrating along nerve tracks as well in the skin and scales (Figure S6). These findings confirm that tfap2b labels McSCs at the niche and that these cells have a multipotent identity with the potential to give rise to all three adult pigment cell types.

Discussion

Here we show that Tfap2b and a select set of its targets specify McSCs from the NC and that these cells reside at the DRG niche. Molecularly, our data point to a subpopulation of cells with MIX+ identity, some of which activate a Tfap2b transcription program to become ErbB-dependent McSCs. We show that Tfap2b is required for melanocyte regeneration from McSCs, providing a functional role for Tfap2b in stem cell potential. When specified, these cells have the potential to generate all three pigment cell types, as well as nerve-associated cells, in the adult pattern.

scRNA-seq combined with live imaging reveals that the spatiotemporal development of pigment cells is complex. Molecularly, our study resolves two developmental pathways for melanocytes: melanoblasts that develop directly from the NC, and MI+ cells that are ErbB dependent. We propose that the melanoblasts contribute to the embryonic pattern in the skin, whereas the MI+ cells are ErbB-dependent crestin mitfa cells that line the peripheral nerves and remain relatively undifferentiated (Figure 1). We also find that the McSC transcriptomes cluster within a larger MIX+ population and propose that not all MIX+ cells in embryogenesis have the potential to contribute to the adult pattern. In contrast, the ErbB-independent MIX+ cells may represent a tripotent precursor cell for the three chromatophore cells in the embryo (Bagnara et al., 1979; Petratou et al., 2018, 2021). Tfap2b and its transcriptional program could therefore distinguish McSCs from other MIX+ cells.

Despite the strong functional evidence for McSC activity at the DRGs, the lack of cell-type-specific markers has prevented investigation of the molecular mechanisms underpinning their biology. Indeed, additional McSC populations may be present in the zebrafish embryo that are dependent on other microenvironment niches, including an ErbB-dependent and blood vessel-associated population dependent on Endothelin factors (Camargo-Sosa et al., 2019). Our lineage tracing analysis indicates that tfap2b McSCs are multipotent and give rise to melanocytes, iridophores, and xanthophores of the adult. These findings provide a foundation for studying how McSCs make fate decisions in growth and replenish pigment cells during tissue regeneration.

Understanding how McSCs are specified and activated in normal development and how they become dysregulated in disease is critical for regenerative medicine and cancer biology. In mammals, multiple McSC populations have been identified at distinct anatomical locations. In the skin, a McSC population residing in the hair follicle is a reservoir for hair and skin melanocytes during the hair cycle and for repigmentation in vitiligo and is a cellular origin of melanoma (Lee and Fisher, 2014; Moon et al., 2017; Sun et al., 2019). In the hair follicle during the resting phase (telogen), McSC populations are found at the bulge region and the secondary hair germ and are functionally heterogeneous (Joshi et al., 2018, 2019). Although both cell populations have melanocyte potential, RNA-seq analysis shows that bulge McSCs express higher levels of NC genes (e.g., Ngfr, Twist1/2, Snai2, and Sox9), whereas secondary hair germ McSCs express higher levels of melanocyte transcriptional network genes (e.g., Sox10, Mitf, Erbb3, Tyr, and Tyrp1), including Tfap2a and Tfap2b. On the palm or sole, the sweat gland serves as a niche for melanocytes and melanoma precursors (Okamoto et al., 2014). In the dermis, a multipotent stem cell that expresses the NC markers NGFR p75 and Nestin is a source of extrafollicular epidermal melanocytes as well as mesenchymal and neuronal cells (Li et al., 2010; Zabierowski et al., 2011). These cell populations may be similar to the multipotent NC-derived Schwann cell precursor (SCP) that resides along the growing nerve, representing a niche for various cell types, including melanocytes (Adameyko et al., 2009; Diener and Sommer, 2020; Ernfors, 2010; Furlan and Adameyko, 2018). The zebrafish McSC anatomical niche suggests that McSCs are functionally analogous to the SCP (Budi et al., 2011; Dooley et al., 2013). Our findings showing that Tfap2b and a select set of target genes specifically marks the McSCs at the DRG will give insight into how nerves provide a niche for melanocyte progenitors (Furlan and Adameyko, 2018). Further, given the enriched Tfap2b expression in hair germ McSCs, our results may have relevance for how hair follicle McSCs selectively contribute to melanocyte regeneration (Joshi et al., 2019).

Melanoma is one of the most aggressive and heterogeneous cancers, and the melanoma transcriptional landscape spans developmental NC and melanocyte lineage signatures, stem cell signatures, and transdifferentiation signatures (Diener and Sommer, 2020; Marine et al., 2020; Patton et al., 2021). We propose that the molecular mechanisms that regulate McSC biology have direct relevance to melanoma pathogenesis. Illustrating this, we recently discovered that the rate of differentiation of McSCs in zebrafish is dependent on a PRL3-DDX21 transcriptional elongation checkpoint and that the same mechanism in melanocyte regeneration portends poor outcomes for individuals with melanoma (Johansson et al., 2020). Importantly, we and others found that Tfap2b is expressed in human and zebrafish melanoma residual disease cell states, a malignant and drug-resistant cell state that contributes to disease recurrence (Marine et al., 2020; Rambow et al., 2018; Shen et al., 2020; Travnickova et al., 2019; Figure S5A). Thus, the developmental Tfap2b mechanism we identify here for zebrafish McSCs could be co-opted in melanoma so that Tfap2b melanoma residual disease cell states may represent a dysregulated McSC developmental lineage.

Limitations of the study

This study uses scRNA-seq coupled with chemical biology to discover the transcriptome of newly established adult McSCs in zebrafish embryos and identifies tfap2b as a functional marker. The current resolution could be improved by performing scRNA-seq at additional time points. We use live imaging and lineage tracing to demonstrate that tfap2b McSCs are the cell of origin for all three pigment cell types and nerve-associated cells. However, very low levels of tfap2b were detected in a small number of cells by scRNA-seq, which may not become McSCs but could contribute to the clonal analysis. An inducible cre recombinase and/or using an integrated promoter approach that uses more than one promoter to target a single cell type could provide further information regarding the fate and function of tfap2b+ McSC derivatives.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, E. Elizabeth Patton (e.patton@ed.ac.uk).

Materials availability

Newly generated materials from this study are available upon request to the Lead Contact, E. Elizabeth Patton (e.patton@ed.ac.uk).

Experimental model and subject details

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

Fish stocks used were: wild-type AB, mitfa (Johnson et al., 2011; Zeng et al., 2015), Tg(nbt:dsRed), Tg(mitfa:GFP) (Dooley et al., 2013), Tg(crestin:mCherry) generated for this study from the plasmid kindly provided from Charles Kaufman (Washington University), Tg(ubi:loxP-GFP-loxP-mCherry; ubi:Switch) (Mosimann et al., 2011), Tg(tfap2b:GFP) (this study). Combined transgenic and mutant lines were generated by crossing. Adult fish were maintained at ∼28.5°C under 14:10 light:dark cycles. Embryos were kept at either 25°C, 28.5°C or 32°C and staged according to the reference table provided by Kimmel and colleagues (Kimmel et al., 1995) or Parichy and colleagues (Parichy et al., 2009). All fish stages are indicated in the figure legends and are as follows (Standard Length-S.L.-expressed in mm): Figure 1: B, 22 hpf; C, 22 hpf, D, 27 hpf; E, 24 or 72 hpf; F, 24 hpf; G, 24, 48 or 120 hpf; Figures 2, 3 and 4: 24hpf, Figure 5: A, 100 hpf; B, 120 hpf, C, 48 hpf; D, 24 hpf; E, 48 hpf; Figure 6: B, 48 hpf; C, 6 dpf; D, 13 dpf (S.L. 4.2); E, 48, 120 hpf or 10 dpf (S.L. 4.1); F, 15 dpf (S.L. 4.9); Figure 7: A, 1 mpf (S.L. 14.2); B, 1 mpf (S.L.11.3); C, 1 mpf (S.L.11.4); D, 1 mpf, Fish 1, S.L. 9.7, Fish 2, S.L. 8.7, Fish 3, S.L. 7.1, Fish 4, S.L. 10.5, Fish 5, S.L. 9.3, Fish 6, S.L. 11.6, Fish 7, S.L. 11.1, Fish 8, S.L. 10.6, Fish 9, S.L. 12.4, Fish 10, S.L. 11.2, Fish 11, S.L. 13.7, Fish 12, S.L. 15.5; E, 1 mpf (S.L. 9.3), F, 1 mpf (S.L. 11.6); G, 1 mpf (10.5 mm); Figure S6: 1 mpf (S.L. 11.6).

Method details

Generation of zebrafish transgenic lines

The promoter region of tfap2b (1kb upstream of the first coding exon) was PCR amplified from zebrafish total genomic DNA with the following set of primers: forward: 5′-GGGGACAACTTTGTATAGAAAAGTTGtacccagagagtcacacatgg-3′; reverse: 5′-GGGGACTGCTTTTTTGTACAAACTTGtGGAATACGCGTGCACTAACAT-3′. Amplicons were cloned into pDONRP4-P1R (Tol2Kit v1.2, plasmid #: 219) to generate the p5E-tfap2b entry clone which were combined with either pME-GFP (Tol2Kit v1.2, plasmid #: 383) (Kwan et al., 2007), and the SV40 polyA sequence from p3E-polyA (Tol2Kit v1.2, plasmid #: 302) into the pDestTol2pA2 (Tol2Kit v1.2, plasmid #: 394) destination vector to generate pEXP(tfap2b:GFP) and expression vectors.

The cre coding sequence was amplified from pMC-CreN plasmid (kind gift Jianguo Shi Check) with the following primers: forward: 5′-ggggacaagtttgtacaaaaaagcaggcttcGCCACCATGCCCAAGAAGAAGAGGAAG-3′; reverse 5′-ggggaccactttgtacaagaaagctgggtcttCTAATCGCCATCTTCCAGCAG-3′. The amplicon was cloned into pDONOR221(Thermo Fisher) producing a middle entry vector, pME-cre that was cloned with the tfap2b promoter from pME-tfap2b, and the SV40 polyA sequence from p3E-polyA into the pDestTol2CG2 destination vector (Tol2Kit v1.2, plasmid #: 395) to generate the pEXPGC2(tfap2b:cre) expression vector.

The pEXP vectors were mixed with Tol2 mRNA (in vitro transcribed with the Ambion mMessage mMachine SP6 Kit, Thermo Fisher, from the Tol2Kit pCS2FA-transposase plasmid–plasmid #: 396); and microinjected into 1-cell stage either AB or ubi:switch embryos, at a final concentration 25 pg/nL and 35 pg/nL respectively. In order to perform lineage trancing with a sufficiently limiting amount of transgene, we also injected pEXPGC2(tfap2b:cre) expression vector at decreasing concentrations (12.5 pg/nL, 6.25 pg/nL, 3.125 pg/nL, or 1.5 pg/nL) with 35 pg/nL of Tol2 transposase.

Zebrafish embryos expressing the tfap2b:GFP transgene were selected and grown to adulthood before crossing with wildtype zebrafish to obtain the F1 generation. Embryos expressing the tfap2b:cre transgene were selected and either grown until the desired developmental stage for confocal imaging or single housed for imaging in epifluorescence.

Zebrafish morpholino oligonucleotides

For the tfap2b morpholino, 1 ng of AUG-directed morpholino (5′-CGTGCACTAACATCTGGGCGGAAAA-3′) or 2 ng of splicing morpholino (5′-GGTGGAAATAATGATAGTCTCACCT-3′) oligonucleotide were injected, as well as the standard control (Gene Tools, LLC). Regenerating melanocytes in tfap2b morpholino-injected and control-injected fish were tested in the mitfa regeneration assays with embryos raised at 32°C for 72 h before down-shifting to 25°C. Imaging and quantification of regenerating melanocytes were performed at 120 hpf. Late-developing melanocytes in tfap2b morpholino-injected and control-injected fish were tested during normal development in AB embryos raised at 28.5°C for 120 h before imaging and quantification of regenerating melanocytes. Representative of 3 biological replicates.

Cryosections

mCherry+ ubi:switch juveniles were cut in the middle of a clearly visible clone before being fixed in 4% PFA overnight at 4°C. The fixed juvenile fish were then soaked in 15% sucrose/PBS overnight at 4°C, then 30% sucrose/PBS for 6 h at 4°C prior to being embedded in O.C.T compound (Sakura) and frozen in isopentane. 40 μm cryosections were cut using a Leica cryostat.

Imaging

Embryos at 4 hpf Tg(mitfa:GFP; nbt:dsRED), Tg(crestin:mCherry; mitfa:GFP) or Tg(tfap2b:GFP; crestin:mCherry) were arrayed in 6-well plates (Corning) containing 0.05% DMSO or 5 μM AG1478 (ErbB-inhibitor, ErbBi, Sigma-Aldrich) in 3 mL of E3 embryo medium and kept at 28°C until imaging time. ubi:switch embryos injected with pEXPGC2(tfap2b:cre) were screened under a fluorescence stereomicroscope for the presence of GFP in the heart at 48 hpf and were then raised as described previously. The same microscope was used to follow 70 embryos injected with either 3.125 or 1.5 pg/nL of pEXPGC2(tfap2b:cre) plasmid. Repeated imaging was performed at 48 hpf, 120 hpf, 10 dpf and 15 dpf. Larvae older than 4 dpf were first soaked in 5 mg/mL epinephrine (Sigma-Adrich) for 5 min, anestethized and quickly mounted in 3% methyl-cellulose prior imaging.

Embryos or fish were selected randomly for confocal imaging as above. Fish older than 5 dpf were imaged while in terminal anesthesia. 1 mpf fish were first soaked in 5 mg/mL epinephrine (Sigma-Adrich) for 10 min prior mounting in low-melting point agarose. Images of randomly picked embryos or selected adult fish were acquired using a 0.4x/0.3, a 10X/0.5 or a 20X/0.75 lens on the multimodal Imaging Platform Dragonfly (Andor technologies, Belfast UK) equipped with 405, 445, 488, 514, 561, 640 and 680 nm lasers built on a Nikon Eclipse Ti-E inverted microscope body with Perfect focus system (Nikon Instruments, Japan). Data were collected in Spinning Disk 40 μm pinhole mode on the Zyla 4.2 sCMOS camera using a Bin of 1 and no frame averaging using Andor Fusion acquisition software. If required: Z stacks were collected using the Nikon TiE focus drive.

Images of cryosections were acquired using 10 and 20X Lenses on a Zeiss Axio-Observer Z1 inverted microscope (Carl Zeiss UK, Cambridge, UK), with a ASI MS-2000 XY stage (Applied Scientific Instrumentation, Eugene, OR). Samples were illuminated using Brightfield or a Lumencor Spectra X LED light source (Lumencor Inc, Beaverton, OR) complete with Chroma #89000ET single excitation and emission filters (Chroma Technology Corp., Rockingham, VT) and acquired on either a Prime BSI Express camera for fluorescence microscopy or a QImaging Retiga R6 Color (Teledyne Photometrics) for color brightfield imaging.

For melanocyte counting, regenerating and normal developing embryos were fixed in 4% PFA/PBS and dehydrated in increasing concentrations of glycerol (Sigma-Adrich). Images were acquired Leica MZFLIII fluorescence stereo microscope with a 1x objective fitted with a Qimaging Retiga Exi CCD camera (Qimaging, Surrey, BC, Canada). Image capture was performed using Micromanager (Version 1.4).

Data were analyzed using Fiji 1.0 and 64bit Java8. Representative of 3 biological repeats.

Single cell suspensions, fluorescence activated cell sorting and library preparation

Tg(crestin:mCherry; mitfa:GFP) were processed in two instances and the following method applied to each treatment separately to obtain two libraries.

300–400 embryos at 4 hpf were divided in two equally sized batches and arrayed in 6-well plates (Corning) containing either E3 or 5 μM AG1478 (ErbBi) in 3 mL of E3 till 24 hpf. A single cell suspension of each batch of embryos was then produced following the method described by Manoli and colleagues (Manoli and Driever, 2012) with minor modifications. Samples were sorted by a FACSAria2 SORP instrument (BD Biosciences UK). Green fluorescence was detected using GFP filters 525/50 BP and 488 nm laser, red fluorescence was detected using RFP filters 582/15 BP and 561 nm laser, and live cells selected with DAPI using DAPI filters 450/20 BP and 405 nm laser. Prior to sorting for fluorescence levels, single cells were isolated by sequentially gating cells according to their SSC-A vs. FSC-A and FSC-H vs FSC-W profiles according to standard flow cytometry practices. Cells with high levels of DAPI staining were excluded as dead or damaged. Cells from wild-type stage matched embryos (without transgenes) were used as negative control to determine gates for detection of mCherry and GFP fluorescence. Then Tg(crestin:mCherry; mitfa:GFP) cells from either untreated or ErbBi-treated zebrafish were purified according to these gates. 10,000 (AB background) fluorescent cells per batch were collected in 100 μL of 0.04% BSA/PBS in LoBind tubes (Fisher Scientific), spun down at 300G at 4°C, resuspended in 34 μL of 0.04% BSA/PBS, and immediately processed using the Chromium platform (10x Genomics) with one lane per sample. Single-cell mRNA libraries were prepared using the single-cell 3′ solution V2 kit (10x Genomics). Quality control and quantification assays were performed using High Sensitivity DNA kits on a Bioanalyzer (Agilent).

Libraries were sequenced on an Illumina NovaSeq platform (1 lane of an S2 flowcell, read 1: 26 cycles, i7 Index: 8 cycles, read 2: 91 cycles). Each sample was sequenced to an average depth of at least 1750 million total reads. This resulted in an average read depth of ∼50,000 reads/cell after read-depth normalisation.

Quantification and statistical analysis

Statistical details of the experiments, n numbers, and dispersion and precision measurements can be found in the figure legends.

Bioinformatics analysis

scRNA-seq data processing and quality check

FastQ files were aligned using the CellRanger (v.2.1.1, 10x Genomics) pipeline to custom zebrafish STAR genome index using gene annotations from Ensembl GRCz11 release 94 with manually annotated entries for GFP, mCherry, mitfa intron 5 and mitfa intron 6 transcripts, filtered for protein-coding genes (with Cell Ranger mkgtf and mkref options). Final cellular barcodes and UMIs were determined using Cell Ranger. Libraries were aggregated (using 10X Cell Ranger pipeline ‘cellranger aggr’ option), with intermediary depth normalization to generate a gene-barcode matrix.

Gene-cell matrices (total: 1519, from untreated embryos: 1022, from ErbBi-treated embryos: 497) were uploaded on R (v.3.6.2) and standard quality control metrics with the Scater package (v.1.12.2) (McCarthy et al., 2017). Only cells with total features >700, log10 total counts >3.0, and mitochondrial gene counts (%) <10 were considered as high quality and kept for further analyses (total:1343, from untreated embryos: 996, from ErbBi-treated embryos: 347). Prediction of the cell cycle phase was performed using the cyclone function in the Scran (v.1.12.1) (Lun et al., 2016).

Clustering, UMAP visualisation and cluster calling

The Louvain clustering of the separated libraries (Figures 2B and 3B) was performed with Seurat (v.3.1.4) (Stuart et al., 2019) using the FindNeighbors and FindClusters functions (cells from untreated embryos: dims = 20, resolution = 0.5; cells from ErbBi-treated embryos: dims = 10, resolution = 0.5) after performing linear dimensionality reduction and checking the dimensionalities of the datasets visualized with elbow plots.

Data were projected onto 2 dimensional spaces using Uniform Manifold Approximation and Projection (UMAP) (McInnes et al., 2018) using the same dimensionality values listed above.

Cluster specific genes were identified using the FindAllMarkers and FindMarkers function in Seurat (v.2.3.4 or v.3.1.4) with default parameters (Wilcoxon Rank-Sum test that compares a single cluster against the others) and then using a Bayesian approach with the scDE package (v.1.99.4) (Kharchenko et al., 2014). See Table S2.

Cluster calling was done after detection of published marker genes for specific cell types and by making unbiased pairwise comparisons based on gene overdispersion against published datasets GEO: GSE112294 (Wagner et al., 2018), GEO: GSE131136 (Saunders et al., 2019), and NCBI SRA: PRNJNA56410 (Farnsworth et al., 2020) using the scMap package (v.1.6.0) (Kiselev et al., 2018) and between the datasets presented in this paper.

The combined clustering for the aggregated libraries used for McSC identification (Figures 3D–3E, S4B–S4D) was performed using the Seurat package (v.2.3.4, dims = 12, resolution = 1) (Butler et al., 2018).

Plots were generated either using Seurat (v.3.1.4) or ggplot2 (v.3.2.1) (Wickham, 2016).

Pseudotime analyses and comparison of the developmental lineages

Differential expression analyses were performed using the identified clusters within each dataset to resolve pseudotemporal trajectories using the setOrderFilter, reduceDimension and orderCells function in Monocle (v.2.12.0) (Trapnell et al., 2014). The minimum spanning tree obtained from cells derived from the untreated or ErbBi embryos was rooted on cluster “twist1a Neural Crest” (Figures 2F, 3C and 3H).

The integrated pseudotime analysis used for the discovery of the McSCs (Figure 3F) was based on the combined clustering (Figures 3D and 3E). The cluster “Other NCC derivatives” was excluded from this analysis because it did not contain pigment cell markers, and we wanted to understand the development of the pigment cell lineage. The top 1000 highly dispersed genes among the untreated embryos dataset were chosen as feature genes to resolve the lineage tree using the setOrderingFilter, reduceDimension, and orderCells functions of Monocle (v.2.12.0). We used the default parameters (except for max_components = 4 and norm_method = log) to generate the 3D trajectory during dimensionality reduction. The same genes were then used to order the cells from the ErbBi-treated embryos and then the combined results were plotted using the PlotComplexTrajectories function and highlight missing states (states 7,8,11 that were collectively called “McSCs” and plotted back in the original UMAPs) in the ErbBi-treated dataset. The transcriptome of the cells belonging to the McSC states from untreated embryos were then compared with the ones from cells of the states composing the pigment lineage (“Axon-associated MIs”, “Directly-developing melanoblasts ”, “MX lineage” and “Pigment progenitors”). The differential expression analysis was performed using a Bayesian approach with the scDE package (v.1.99.4) and the results were plotted using the ggplot2 package (v.3.2.1) and the pathway analysis was performed using the ClusterProfiler R package (v.3.12.0). The same approach was used for the other differential expression analyses presented.

RNA velocity analyses

RNA velocity analyses were performed with the Velocyto R package (v.0.6) (La Manno et al., 2018) using default parameters.

Tfap2b targets

The Tfap2b (Biotin) ChIP-seq data was retrieved from the Gene Expression Omnibus (GEO) under the accession code GEO: GSE125711. The same mapping and peak calling pipeline described by the dataset-linked publication was used, with the setting of FDR <0.01, fold enrichment >2 set to define the final Tfap2b binding element region. R Package “ChIPseeker” (v1.26.2) (Yu et al., 2015) was used to map the peak coordinates to the gene symbols using chicken genome assembly galGal5. Package “homologene” (v1.5.68.21.2.14) was used to identify human homologous for the gene targets identified in Tfap2b (Biotin) ChIP-seq, as well as the differential expressed genes in the zebrafish cluster (with human homolog ZEB2 manually mapped to zebrafish gene zeb2a as we found it was not automatically mapped by the software).

Other statistical analyses

Counts of dorsal melanocytes in the head and trunk region were performed using the Cell Counter plugin on ImageJ Fiji.

The niche occupancy was calculated as the percentage of GFP+ positive cells per number of visible DRG (Figure 1) or as the number of fluorescent cells ventrally to the neural tube (Figure 5). Two images per embryos were acquired (caudally to the first somite, ∼8 somites; rostrally to the urogenital opening, ∼8 somites) and the total number of the niches per embryo were considered.

Statistics for regeneration assays and niche occupancy were performed using running R (v.3.6.2) from RStudio (v.2). For all assays, a normal distribution and equal variance were assumed. For assays with more than two groups, data was analyzed through Analysis of variance (ANOVA), using Tukey-HSD (Honestly Significant Difference) test. Box plots: boxes represent 25th to 75th percentiles, lines are plotted at median. Whiskers represent Min to Max.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
AG1478	Calbiochem/Sigma-Aldrich	Cat No: #658552
FACSMax	AMSbio	Cat No: AMS.T200100
(+/-)–Epinephrine hydrochloride	Sigma-Aldrich	Cat No: E4642-5G
Critical commercial assays
Chromium Single Cell 3′ Library & Gel Bead Kit v2, 16rxns	10x Genomics	Cat No: PN-120237
Chromium Single Cell A Chip Kit, 16 rxns	10x Genomics	Cat No: PN-10000009
Chromium i7 Multiplex Kit, 96 rxns	10x Genomics	Cat No: PN-120292
Tol2kit gateway cloning	(Kwan et al., 2007)http://tol2kit.genetics.utah.edu/index.php/Main_Page	N/A
Deposited data
scRNA-seq data	This paper	GEO: GSE178364
Zebrafish single cell RNA-seq – 1 dpf	(Wagner et al., 2018)https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3067195)	GEO; GSE112294 /GSM3067195
Zebrafish single cell RNA-seq – 1 dpf	(Saunders et al., 2019)https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131136	GEO: GSE131136/GSM3764579
Zebrafish single cell RNA-seq – 5 dpf	(Farnsworth et al., 2020)http://cells.ucsc.edu/?ds=zebrafish-dev	NCBI SRA: PRJNA564810
Zebrafish melanoma datasets	(Travnickova et al., 2019)	GEO: GSE136900
Tfap2b (biotin) ChIP-seq	(Ling and Sauka-Spengler, 2019)	GEO: GSE125711
Experimental models: Organisms/strains
mitfavc7/vc7	(Johnson et al., 2011)	RRID:ZFIN_ZDB-GENO-110330-3
Tg(nbt:dsred) – Tg(Xla.Tubb:dsred)	Professor David LyonsEdinburgh University	ZFIN: ZDB-TGCONSTRCT-081023-2
Tg(mitfa:gfp)	(Dooley et al., 2013)	ZFIN: ZDB-TGCONSTRCT-081203-1
Tg(crestin:mCherry)	This paper from plasmid generated by Kaufman and colleagues(Kaufman et al., 2016)	ZFIN: ZDB-TGCONSTRCT-160208-3
Tg(ubi:loxP-eGFP-loxP-mCherry) - ubi:Switch	(Mosimann et al., 2011)	ZFIN: ZDB-FISH-201123-10
Tg(tfap2b:GFP)	This paper	N/A
Tg(tfap2b:cre)	This paper	N/A
Oligonucleotides
tfap2b promoter F 5′-GGGGACAACTTTGTATAGAAAAGTTGtacccagagagtcacacatgg-3′	Sigma-Aldrich	N/A (custom made)
tfap2b promoter R 5′-GGGGACTGCTTTTTTGTACAAACTTGtGGAATACGCGTGCACTAACAT-3′	Sigma-Aldrich	N/A (custom made)
cre cds F 5′-ggggacaagtttgtacaaaaaagcaggcttcGCCACCATGCCCAAGAAGAAGAGGAAG-3′	Sigma-Aldrich	N/A (custom made)
cre cds R 5′-ggggaccactttgtacaagaaagctgggtcttCTAATCGCCATCTTCCAGCAG-3′	Sigma-Aldrich	N/A (custom made)
tfap2b MO (AUG-directed):5′-CGTGCACTAACATCTGGGCGGAAAA-3′	GeneTools LLC	N/A (custom made)
tfap2b MO (Splicing site-directed):5′-GGTGGAAATAATGATAGTCTCACCT-3′	GeneTools LLC	N/A (custom made)
Standard control morpholino	GeneTools LLC	N/A
Recombinant DNA
pDONOR221	Tol2kit v1.2	plasmid #:218
pDONRP4-P1R	Tol2kit v1.2	plasmid #:219
p5E-actin2	Tol2kit v1.2	plasmid #:299
p3E-polyA	Tol2kit v1.2	plasmid #:302
pME-GFP	Tol2kit v1.2	plasmid #: 383
pDestTol2pA2	Tol2kit v1.2	plasmid #: 394
pDestTol2CG2	Tol2kit v1.2	plasmid #:395
p5E:tfap2b	This paper	N/A
p5E:cre	This paper	N/A
pEXP(tfap2b:GFP)	This paper	N/A
pEXPGC2(tfap2b:cre)	This paper	N/A
Software and algorithms
CellRanger (v.2.1.1)	10x Genomics	RRID:SCR_017344
Scater (v.1.12.2)	(McCarthy et al., 2017)https://bioconductor.org/packages/release/bioc/html/scater.html	RRID:SCR_015954
Scran (v.1.12.1)	(Lun et al., 2016)https://bioconductor.org/packages/release/bioc/html/scran.html	N/A
Seurat (v.2.3.4)	(Butler et al., 2018)https://satijalab.org/seurat/	RRID:SCR_016341
Seurat (v.3.1.4)	(Stuart et al., 2019)https://satijalab.org/seurat/	RRID:SCR_016341
scDE (v.1.99.4)	(Kharchenko et al., 2014)https://hms-dbmi.github.io/scde/	RRID:SCR_015952
Monocle (v.2.12.0)	(Trapnell et al., 2014)http://cole-trapnell-lab.github.io/monocle-release/docs/	RRID:SCR_016339
ClusterProfiler (v.3.12.0)	(Yu et al., 2012) https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html	RRID:SCR_016884
Velocyto R package (v.0.6)	(La Manno et al., 2018)http://velocyto.org/	RRID:SCR_018167
ChIPseeker (v.1.26.2)	(Yu et al., 2015)https://bioconductor.org/packages/release/bioc/html/ChIPseeker.html	N/A
Homologene (v1.5.68.21.2.14)	http://www.ncbi.nlm.nih.gov/homologene	RRID:SCR_002924
ggPlot2 (v.3.2.1)	(Wickham, 2016)https://cran.r-project.org/web/packages/ggplot2/index.html	RRID:SCR_014601
R (v.3.6.2)	http://www.r-project.org/	RRID:SCR_001905
RStudio (v.2)	http://www.rstudio.com/	RRID:SCR_000432
Micromanager (Version 1.4)	http://micro-manager.org	RRID:SCR_00041
Fiji 1.0	http://fiji.sc	RRID:SCR_002285