TFAP2 paralogs facilitate chromatin access for MITF at pigmentation and cell proliferation genes

In developing melanocytes and in melanoma cells, multiple paralogs of the Activating-enhancer-binding Protein 2 family of transcription factors (TFAP2) contribute to expression of genes encoding pigmentation regulators, but their interaction with Microphthalmia transcription factor (MITF), a master regulator of these cells, is unclear. Supporting the model that TFAP2 facilitates MITF’s ability to activate expression of pigmentation genes, single-cell seq analysis of zebrafish embryos revealed that pigmentation genes are only expressed in the subset of mitfa-expressing cells that also express tfap2 paralogs. To test this model in SK-MEL-28 melanoma cells we deleted the two TFAP2 paralogs with highest expression, TFAP2A and TFAP2C, creating TFAP2 knockout (TFAP2-KO) cells. We then assessed gene expression, chromatin accessibility, binding of TFAP2A and of MITF, and the chromatin marks H3K27Ac and H3K27Me3 which are characteristic of active enhancers and silenced chromatin, respectively. Integrated analyses of these datasets indicate TFAP2 paralogs directly activate enhancers near genes enriched for roles in pigmentation and proliferation, and directly repress enhancers near genes enriched for roles in cell adhesion. Consistently, compared to WT cells, TFAP2-KO cells proliferate less and adhere to one another more. TFAP2 paralogs and MITF co-operatively activate a subset of enhancers, with the former necessary for MITF binding and chromatin accessibility. By contrast, TFAP2 paralogs and MITF do not appear to co-operatively inhibit enhancers. These studies reveal a mechanism by which TFAP2 profoundly influences the set of genes activated by MITF, and thereby the phenotype of pigment cells and melanoma cells.


Introduction
Gene expression in developing melanocytes and melanoma, a cancer derived from the melanocyte lineage, is regulated by transcription factors including Microphthalmia-associated transcription factor (MITF) and members of the SOXE, PAX and TFAP2 families [1][2][3][4][5][6][7][8][9][10][11]. MITF is required for differentiation of melanocytes during development, and its activity is regulated at both the transcriptional and post-translational levels [3]. In melanoma cells, high levels of MITF activity promote cell proliferation and pigmentation, while lower levels promote an invasive phenotype [12,13]. Mass spectroscopy revealed that MITF interacts with components of both the PBAF chromatin remodeling complex, including BRG1 and CDH7, and the NURF remodeling complex, including RBBP4 [14,15]. Furthermore, chromatin immunoprecipitation of BRG1 in cells depleted of MITF revealed that MITF recruits BRG1 to the promoters of specific genes, including TYR, which encodes the rate-limiting enzyme in melanin synthesis Tyrosinase [15]. Similar analysis suggested that SOX10 also recruits BRG1 to chromatin, and at some loci it does so in co-operation with MITF [15]. Conversely, there is evidence that PAX3 inhibits the activity of MITF at the DCT promoter [16]. Furthermore, low MITF activity is associated with an invasive phenotype, and deletion or knockdown of MITF results in upregulation of genes that promote migration and invasion [17]. MITF CUT&RUN peaks are found near some genes whose expression is upregulated in MITF mutant cells, implying that MITF directly represses their expression [17]. This set of MITF peaks is enriched for the binding site of FOXC1, a transcriptional repressor [18], suggesting MITF has co-factors in its repressive function as well as its activating one.
The activating enhancer-binding family of transcription factors, comprising five members, TFAP2A-E, regulate development of many cell types and organs including neural crest, placodes, epidermis, trophectoderm, heart, kidney, and brain [19][20][21][22][23][24][25][26][27][28]. In several contexts, including melanocyte differentiation, TFAP2 paralogs function redundantly [10,[29][30][31]. For instance, in zebrafish tfap2a loss-of-function mutant embryos the number of melanocytes is lower than normal and pigmentation is delayed relative to in wild-type embryos; this phenotype is exacerbated if tfap2a mutant embryos are also depleted of tfap2e expression with antisense morpholinos [10]. In zebrafish melanoma Tfap2a and Tfap2e also appear to act redundantly to promote proliferation and, interestingly, to suppress cell adhesion and cell migration [32]. Consistent with this, in the skin of mouse embryos with neural-crest specific knockout of Tfap2a and Tfap2b, the two paralogs with highest expression, fewer-than-normal cells express markers of melanocytes [33]. In addition, evidence exists for sub-functionalization among Tfap2 paralogs [34]. Of note, while tfap2e is the paralog with highest expression in zebrafish embryonic melanophores, its clear ortholog in mice, Tfap2e, is not expressed in embryonic melanocytes [35]. This suggests that during evolution there has been some shuffling of the function of individual TFAP2 paralogs among species.
TFAP2 paralogs and MITF appear to co-activate certain genes. For instance, in a human melanoma cell line, the in vitro enhancer activity of an element within an IRF4 intron depended on the simultaneous binding of MITF and TFAP2 [36]. Further, in zebrafish tfap2a and mitfa double mutant embryos there is a greater-than-additive reduction in the number, and level of pigmentation, of melanocytes in comparison to single mutants [33]. Supporting parallel function of Tfap2 paralogs and Mitfa, the promoters of MITF target genes are enriched for TFAP2 consensus binding sites [15,37], and ChIP-seq experiments in primary melanocytes indicate that TFAP2A and MITF bind overlapping regions of chromatin near genes encoding regulators of pigmentation [33]. Collectively, these observations indicate that TFAP2 paralogs co-activate a subset of MITF target genes by binding at the same enhancers. Still unclear, however, is whether TFAP2 paralogs and MITF act cooperatively or independently at enhancers they co-regulate. It is also unclear whether they co-repress enhancers.
TFAP2 paralogs may serve as pioneer factors for MITF, although not all evidence supports this possibility. Pioneer or initiating TFs can bind nucleosome-bound DNA and recruit other TFs that lack this property called settler TFs (reviewed in [38,39,40]). Evidence that TFAP2 paralogs are pioneer factors includes, first, that the TFAP2 binding site is over-represented within DNase1-protected "footprints" in mouse embryonic stem cells induced to differentiate [41]. Second, TFAP2A catalyzes assisted loading of androgen receptor (AR) in epididymis cells [42] and estrogen receptor in MCF-7 cells [43]. Third, the TFAP2 binding site is enriched for at the center of ATAC-seq peaks, implying it has a strong effect on chromatin accessibility [44]. Fourth, ATAC-seq peaks in naïve-stated human ESC showed reduced openness in TFAP2C-KO cells [45], and forcing expression of TFAP2C in human ESC is sufficient to open chromatin at loci where it binds [46]. Finally, TFAP2A, TFAP2B and TFAP2C can bind nucleosomes [47]. Together these findings support the possibility that TFAP2 displaces nucleosomes and thereby facilitates chromatin binding by MITF. However, it is not clear that MITF needs a pioneer factor to bind chromatin. In the dynamic-assisted-loading model, all classes of TFs have short residency on chromatin (reviewed in [38]). In this model, initiating TFs can recruit ATP-dependent chromatin remodelers (nBAF, SWI/SNF, INO80, ISWI, NURD) and thereby make chromatin accessible to other TFs, i.e., the assisted TFs [48]. Unlike the pioneer factor model, the assisted-loading-model predicts that initiation factors are interchangeable, and that the recruitment of chromatin remodeling complexes by initiation factors is dependent on the local chromatin structure and accessibility [48]. As mentioned above, MITF binds various components of the SWI/SNF complex [14,49,50] and the chromatin remodeler CHD7 [15] and so meets the criteria for an initiating factor. If the dynamic-assisted-loading model holds in this situation, MITF would have no need for a pioneer factor like TFAP2 to assist its binding to chromatin.
To address these questions, we used single cell RNA-sequencing (scRNA-seq) to investigate the temporal order of gene expression of mitfa, tfap2 paralogs, and pigmentation regulators in cells of the melanocyte lineage which were included among GFP-expressing cells isolated from Tg(mitfa-GFP) zebrafish embryos at 28 hours post fertilization (hpf). Furthermore, to test how TFAP2 paralogs affect enhancer activity and MITF binding, we deleted TFAP2 paralogs from SK-MEL-28 melanoma cells and assessed: nucleosome positioning, using the assay for transposase-accessible chromatin with next generation sequencing (ATAC-seq); enhancer activity, using cleavage under targets and release using nuclease (CUT&RUN) with anti-H3K27Ac, anti-H3K4Me3, and anti-H3K27me3; and binding of MITF, using CUT&RUN with anti-MITF. We similarly assessed binding of TFAP2A in wild-type cells and in those harboring loss of function mutations in MITF. Our results support the notion that TFAP2 factors behave like the canonical pioneer factor FOXA1: at many chromatin elements bound by TFAP2A, loss of TFAP2 led to loss of enhancer activity, and in a large subset, it also resulted in chromatin becoming condensed. At both subsets of TFAP2-activated enhancers, MITF binding was TFAP2-dependent. In addition, we find evidence that TFAP2 paralogs can also inhibit enhancers, and at a subset of those, they exclude binding of MITF. Finally, the analyses suggest that TFAP2 directly inhibits many of the same genes that MITF inhibits, but we do not find evidence that TFAP2 and MITF co-repress the same enhancers. Together these findings illuminate the mechanisms by which TFAP2 and MITF coordinately regulate differentiation of melanocytes and the phenotype of melanoma cells.

Expression of tfap2 paralogs in the melanocyte lineage precedes activation of Mitfa transcriptional programs in-vivo
If Tfap2 paralogs cooperate with Mitf to promote expression of melanocyte differentiation genes in-vivo, then high level expression of known Mitfa-target genes, including dct, pmel, mreg (zgc:91968) and trpm1b, should only occur in cells that express both mitfa and tfap2 paralogs. To test this prediction, we performed scRNA-seq on GFP + cells sorted from Tg(mitfa: GFP) transgenic zebrafish embryos at 28 hours post fertilization (hpf) using the 10x Chromium system. We sequenced 11,217 cellular transcriptomes and visualized the data in twodimensional space using Seurat [51] and uniform and manifold approximation and projection (UMAP) [52] (Fig 1A). We assigned cell-type annotations using previously identified marker genes [53][54][55][56] and information from the Zebrafish Information Network (ZFIN) [57]. The 28 clusters comprised 11 cell types (Fig 1A); neural crest and basal cell lineages were the most strongly represented (S1A-S1C Fig). Interestingly, only a minority of the clusters expressed mitfa. Sorted cells that express GFP but not mitfa presumably reflect either non-specific activity of the mitfa promoter used in the transgene or earlier transient expression of mitfa revealed by the long half-life of GFP; importantly, for the upcoming analysis, such cells did not express pigmentation genes. We focused on sox10-positive clusters 6-12 as these comprised all stages of the melanocyte lineage, including all mitfa-expressing cells (S1D Fig). Re-clustering of sox10-positive cells identified ten clusters (Fig 1B). Two of these were neural crest (NC) cell populations which shared expression of bona fide NC markers (foxd3 and sox10), but were distinguished by expression in one cluster of members of the Notch signaling pathway (her6, her4 and her12); this cluster may represent migrating cardiac or cranial NC populations [58] (Figs  1C, S2A and S2D). In addition, two clusters expressed markers of a tripotent precursor of melanoblasts, iridoblasts, and xanthoblasts (cdkn1ca, slc15a2, ino80e, id3, mycn, tfec), called MIX cells [59][60][61][62] (S2B Fig). Interestingly, while both MIX clusters express tfap2a, the clusters differed by one and not the other expressing tfap2e and tfap2b -indeed these genes were ranked first and third, respectively, among those differentially expressed between the two MIX clusters (Figs 1G, S3A and Table A in S1 Table). Therefore, we will refer to the two MIX clusters as tfap2-low and tfap2-high (Fig 1B). An additional cluster expressed high levels of melanoblast/ xanthoblast markers (mitfa, impdh1b, gch2, id3) suggesting it corresponds to the previously described MX cluster [59][60][61][62]. Another cluster corresponding to differentiated   (Fig 1B). Within this trajectory, despite expression of mitfa in both tfap2-low and tfap2-high MIX clusters, expression of select pigmentation genes (dct, pmel, mreg and trpm1b) is first detected in tfap2high MIX cells, and becomes higher in MX cells concomitant with higher level expression of mitfa and of tfap2e (Figs 1C-1G and S4A and S4B). These data support a cooperative role for Mitfa and Tfap2 paralogs in activation of pigmentation genes.
Importantly, the phenotypes in single or double mutants thus far examined may be suppressed by compensatory upregulation of other Tfap2 paralogs. Indeed, we detected upregulated expression of tfap2c in tfap2a mutants and in tfap2a/tfap2e double mutants, presumably driven by non-sense mediated decay (S6 Fig) [91]. In conclusion, the in-vivo consequence of loss of Tfap2 function upon melanophore differentiation is difficult to assess because of the potential for redundant function among Tfap2 paralogs and because of pleiotropy (see embryos at 28 hours post fertilization (hpf). Annotated cell clusters as labelled. (B) UMAP obtained after re-clustering sox10-expressing clusters 6-12 (n = 1918 cells). Black line-Monocle pseudotime trajectory analysis starting at foxd3+ neural crest cells showing the progression through different pigment cell clusters as shown. (C-F) Violin plots showing expression of select genes foxd3, mitfa, tfap2e and dct for each cell cluster represented in B. Blue arrows point to the MIX-tfap2-low cluster, and show that mitfa, but not pigmentation genes, are expressed in this cluster (tfap2a is also expressed in this cluster). Red arrows point to the MIX-tfap2-high cluster, and show that mitfa, tfap2e, and pigmentation genes are expressed in this cluster. (G) Dot plot representing the expression of sox10, tfap2 paralogs, mitfa and Mitfa-target genes, and the percentage of cells expressing these genes in neural crest, MIX, MX and melanophore clusters. Expression of tfap2 paralogs are highlighted in MIX tfap2-low and MIX tfap2-high clusters (red box). Size of dots represents percentage of cells expressing the gene, and blue-red scale (low-high) represents relative average expression among cells. (H-K) Lateral views of head and trunk of live embryos at 29 hpf, anterior to the left and dorsal to the top. Genotype as shown. Boxes, regions magnified in accompanying panels H'-K.' (H-H') A wild-type or heterozygous mutant (sibling) embryo with normal melanophores (white arrowheads). (I-I') A tfap2e ui157/ui157 embryo, with melanophores that are normal in terms of number, differentiation and pigmentation (white arrowheads). (J-J') A tfap2a low/low homozygous mutant embryo, with fewer melanophores than tfap2e ui157ui/157 and WT sibling embryos. (K-K') A tfap2a low/low ; tfap2e ui157ui/157 double-mutant embryo, with fewer and paler melanophores than in tfap2a low/low siblings. (L) Box plot illustrating the number of pigmented melanophores in the dorsum of tfap2a low/low , tfap2a low/low ; tfap2e +/ ui157 , and tfap2a low/low ; tfap2e ui157/ui157 double mutant embryos at 36 hpf. Center line, mean; box limits, upper and lower quartiles; whiskers, minimum and maximum values; black dots, number of melanocytes per individual embryo (tfap2a low/low ; n = 9, tfap2a low/low ; tfap2e +/ui157 ; n = 32, tfap2a low/low ; tfap2e ui157/ui157 , n = 10). P-value according to the Student's t-test.

TFAP2A binds open and closed chromatin
To test how TFAP2 paralogs interact with MITF we switched to a cell line model. We chose the SK-MEL-28 melanoma cell line because we have deleted all alleles of MITF from this cell line and have evaluated MITF binding genome-wide using CUT&RUN [17]. Here we carried out (1) CUT&RUN using antibodies to TFAP2A (i.e., TFAP2A peaks), (2) CUT&RUN using antibodies to chromatin marks indicative of active regulatory elements (H3K27Ac and H3K4Me3) [67,68] and indicative of inactive chromatin (H3K27Me3) [69], and (3) ATAC-seq to distinguish between open and closed chromatin [70]. We used IgG as a background control and the MACS2 software to call peaks in each CUT&RUN dataset (S7A and S7B Fig). Based on proximity to transcriptional start sites (TSS), about one-third of TFAP2A peaks appeared proximal to promoters (within 3kb of a TSS). As expected, these elements had strong H3K4Me3 signal (S7B Fig). At promoter-proximal TFAP2A peaks, the H3K27Ac signal in WT cells was relatively consistent, whereas at promoter-distal TFAP2A peaks (greater than 3 kb but within 100 kb of a TSS) the H3K27Ac signal ranged from high to background level (S7B and S7C Fig). About two-thirds of TFAP2A peaks overlapped ATAC-seq peaks, indicating that they were in open chromatin (S7D and S7E Fig). Of note, the read depth (height) of a peak approximates the number of chromosome molecules where TFAP2A binds. The average read depth of the TFAP2A peaks in closed chromatin was only about 50% of that in open chromatin but was nonetheless 80-fold higher than the IgG background read depth (S7D, S7E, S8A and S8B Figs for example loci). Importantly, the TFAP2 binding site was strongly enriched for in both TFAP2A-bound elements where the local ATAC-seq signal was called as a peak and in counterparts where it was not (p < 1 x 10 −1785 and p < 1 x 10 −4375 , respectively), supporting the idea that TFAP2A binds DNA directly even when the DNA is occupied by nucleosomes (S8C and S8D Fig). These results indicate that TFAP2A binds at both open and closed chromatin, consistent with it being a pioneer factor, and at enhancers and promoters with a range of activity levels.

At a subset of TFAP2-activated enhancers, TFAP2 is necessary for chromatin accessibility
Having shown that TFAP2A can bind closed chromatin, we next asked whether TFAP2 factors play a role in opening chromatin. We used CRISPR/Cas9 methods to introduce frame-shift mutations into TFAP2A and TFAP2C, the TFAP2 genes with high expression in SK-MEL-28 cells (S9A-S9D Fig). We then carried out RNA-seq, ATAC-seq, and CUT&RUN with antibodies to H3K27Ac, H3K4Me3, and H3K27Me3, in two independent knockout clones (hereafter, TFAP2-KO cells). Western blot analysis showed an absence of immunoreactivity for both proteins (S9E Fig). Control clones (hereafter WT cells) were derived from the parental SK-MEL-28 line transiently transfected with Cas9 but not with guide RNAs. RNA-seq revealed that expression of 532 genes was downregulated, and expression of 609 genes was upregulated, in TFAP2-KO cells (i.e., in both clones) versus in WT cells. We will refer to these sets as TFA-P2-activated genes and TFAP2-inhibited genes, respectively.
We defined TFAP2-activated enhancers as TFAP2A peaks a) that overlap H3K27Ac and ATAC-seq peaks in WT cells (21,745/ 36,948, or 59% of all TFAP2A peaks), b) that are greater than 3 kb from a transcription start site (11,005/ 21,745, or 51% of TFAP2A peaks at H3K27Ac/ATAC peaks), and c) where the H3K27Ac signal was significantly lower in TFAP2-KO cells relative to in WT cells (adj p < 0.05, log2FC <-1; 3,858/11,005 or 35% of TFAP2Abound enhancers). Interestingly, at about half the TFAP2-activated enhancers (2002/3858), the ATAC-seq signal was significantly lower in TFAP2-KO versus in WT cells (adj. p < 0.05, log2FC <-1) (Figs 2 and 3A, 3E-3E'). The open status of these loci depends on TFAP2 paralogs, and we therefore infer that TFAP2 paralogs directly or indirectly recruit the nucleosomedisplacing machinery. We further infer that TFAP2 paralogs bound these loci as closed chromatin, but it is also possible they bound to another pioneer factor present at these loci. In either case, it participates in a pioneering function, and we refer to this subset of enhancers as TFAP2-pioneered-and-activated enhancers (Fig 2 Box 1A). We refer to the subset where ATAC-seq signal is unchanged between TFAP2-KO and WT cells as non-pioneered TFA-P2-activated enhancers (Fig 2 Box 1B and Fig 3B, 3F-3F').
Both TFAP2-pioneered-and-activated enhancers and non-pioneered TFAP2-activated enhancers were associated with TFAP2-activated genes. Interestingly, the association was stronger for the former subset (Fig 3I and 3J and Table 1, compare rows 4 and 5). Moreover, at both subsets, the H3K4Me3 signal, which is associated with enhancer activity [68], was reduced in TFAP2-KO cells relative to WT cells (Fig 3E" and 3F"). While both subsets were strongly enriched for the TFAP2 binding site and certain other binding sites (e.g., RUNX), the subset pioneered by TFAP2 was more strongly enriched for the SOXE and MITF binding sites, while the non-pioneered subset was more strongly enriched for the FRA1, TEAD and the ZFX binding sites (Fig 3M and 3N). Of note, FRA1 is a pioneer factor [71] which could explain why these elements do not depend on TFAP2 to be free of nucleosomes.

TFAP2 paralogs inhibit enhancers by blocking the opening of chromatin
Because of evidence that TFAP2A directly represses gene expression [72][73][74] we next sought to identify enhancers directly inhibited by TFAP2 paralogs. To this end we filtered promoter-distal TFAP2A peaks for those where the local H3K27Ac signal was higher in TFAP2-KO cells than in WT cells (adj. p <0.05, log2FC>1). Analogously to TFAP2-activated enhancers, TFA-P2-inhibited enhancers were split between a subset where the ATAC-seq signal was higher in TFAP2-KO cells than in WT cells (as illustrated in Figs 2 Box 2A, 3C and 3G-3G'), implying TFAP2 paralogs maintain condensed chromatin at these sites, and a subset where it was unchanged (Figs 2 Box 2B, 3D and 3H-3H'). The first but not the second subset was significantly associated with TFAP2-inhibited genes (Fig 3K and Table 1), and the average H3K4Me3 signal at the first but not the second subset was higher in TFAP2-KO cells than in WT cells (Fig 3G" and 3H"). We define these enhancers as TFAP2-pioneered-and-inhibited. The binding site for TFAP2 was strongly enriched for in these enhancers, as were those for ETS1 and CTCF (Fig 3O), both transcriptional repressors [75,76]. By contrast, the subset of candidate TFAP2-inhibited enhancers where the ATAC-and average H3K4Me3 signals were unchanged between TFAP2-KO and WT cells was not associated with TFAP2-inhibited genes. We infer these elements are not, in fact, TFAP2-inhibited enhancers. In conclusion, at TFA-P2-inhibited enhancers TFAP2 recruits and/or retains a machinery that condenses chromatin and inhibits enhancer activity; the canonical pioneer factor FOXA1 also has this potential [77,78]. Alternatively, TFAP2 itself may instead play a distinct role in keeping enhancers closed rather than actively closing them from a presumed prior open state.
We analyzed whether TFAP2 paralogs directly activate or directly inhibit promoters, and how they do so. Although TFAP2A peaks are frequently found at promoters (8277 genes have a TFAP2A peak within 3 kb of the TSS), it was uncommon for the underlying H3K27Ac and H3K4Me3 signals to be elevated (or reduced) in TFAP2-KO cells relative to WT cells. Based on this quality, we found 119 candidates for directly TFAP2-activated promoters, and 31 TFAP2-pioneered-and-activated enhancers show reduced nucleosome accessibility (ATAC-seq) and reduced levels of active chromatin marks (H3K27Ac and H3K4Me3) in TFAP2-KO cells compared to in WT cells. We infer that TFAP2 paralogs pioneer chromatin access for transcriptional co-activators, like MITF and SOX10 (purple box), that in turn recruit chromatin remodelling enzymes and histone modifying enzymes. (1B) Non-pioneered TFAP2-activated enhancers show loss of active candidates for directly TFAP2-inhibited promoters. Nonetheless, similar to the trends for TFAP2 regulated enhancers, the pioneered subset of directly TFAP2-activated promoters was more strongly associated with TFAP2-activated genes than the non-pioneered subset, and only the pioneered subset of TFAP2-inhibited promoters was associated with TFAP2-inhibited genes (S10A-S10L Fig and additional examples in S11A and S11B Fig).
Interestingly, at the majority of TFAP2A peaks, whether at open or closed chromatin, neither the ATAC-seq nor H3K27Ac signals were TFAP2-dependent. This large subset of quiescent TFAP2A peaks was associated with neither TFAP2-activated nor TFAP2-inhibited genes ( Table 1, row 10). In contrast to the quiescent subset of peaks, the set of all TFAP2A peaks was modestly associated with TFAP2-activated genes ( Table 1, row 9). Nonetheless, we infer that the common practice of using ChIP-seq or CUT&RUN-seq data to identify genes directly regulated by a transcription factor is subject to false positives.

At a subset of MITF/TFAP2 co-bound peaks, TFAP2 paralogs facilitate MITF binding
A prediction of the TFAP2-as-pioneer-factor model is that binding of transcription factors, like MITF, will depend on TFAP2. Among 36,621 MITF peaks in WT SK-MEL-28 cells that we previously identified by CUT&RUN [17], we found that 15,752 (43%) overlap a TFAP2A peak. Of these, 9,413 (60%) were within open and active (i.e., ATAC+ and H3K27Ac+) chromatin (S12A Fig). To assess MITF binding in the absence of TFAP2, we carried out anti-MITF CUT&RUN in TFAP2-KO cell lines. Of note, as MITF RNA levels in TFAP2-KO cells are only about 40-50% of those in WT cells, an across-the-board decrease in the average height (read depth) of MITF peaks was possible. Instead, we observed that the average height of MITF peaks not overlapping TFAP2A peaks was equivalent in TFAP2-KO cells and in WT (S12B and S12D Fig for example at MLANA). By contrast, the height of about 35% (5,443 /15,752) of MITF peaks overlapping TFAP2A peaks was significantly lower in TFAP2-KO cells than in WT cells (adj. p <0.05, log2FC < -1). (Figs 4A-4D, S12C-S12D and all replicates shown in S13 Fig). We refer to these as directly TFAP2-activated MITF peaks.
We reasoned that TFAP2 paralogs could facilitate MITF binding by displacing nucleosomes (i.e., in pioneer factor mode) or alternatively by elevating MITF's affinity for open DNA. Consistent with both models, we observed that directly TFAP2-activated MITF peaks fall in three subsets with respect to the TFAP2-dependence of the underlying ATAC-seq signal. Thus, at about 57% (3,083/5,443) of directly TFAP2-activated MITF peaks the ATAC-seq signal was significantly lower (Fig 4B, 4C and 4F-4F'), at 37% (2,022/5,443) it was unchanged (Fig 4D  and 4H-4H'), and at 6% it was higher (S14A and S14A' Fig) in TFAP2-KO cells compared to chromatin marks but unchanged nucleosome accessibility in TFAP2-KO cells compared to in WT cells. At these enhancers, we infer that TFAP2 paralogs recruit the binding of transcription factors that, in turn, recruit histone modifying enzymes. TFAP2 paralogs also may recruit such enzymes. It is possible that these elements are stably pioneered by TFAP2 paralogs [94]. (1C) At TFAP2-independent elements, neither the nucleosome accessibility nor active histone marks are altered in TFAP2-KO cells relative to in WT cells. Both types of TFAP2-activated enhancer are significantly enriched near genes whose expression is reduced in TFAP2-KO cells relative to in WT cells (i.e., TFAP2-activated genes). Such genes are associated with the gene ontology (GO) terms cell proliferation and pigmentation. TFAP2-independent elements are associated with neither TFAP2-activated nor TFAP2-inhibited genes. Box 2: TFAP2A peaks at closed and inactive chromatin. (2A) TFAP2-pioneered-and-inhibited enhancers show increased nucleosome accessibility and increased levels active chromatin marks in TFAP2-KO cells compared to in WT cells. At these sites we infer that TFAP2 paralogs recruit or stabilize the binding of enzymes that condense chromatin and that inhibit the binding of transcriptional activators that are otherwise inclined to bind at them. These elements are significantly enriched near genes whose expression is elevated in TFAP2-KO cells relative to in WT cells (i.e., TFAP2-inhibited genes). Such genes were associated with the GO terms cell-cell adhesion and cell migration. (2B) At TFAP2-independent elements, neither the nucleosome accessibility nor active histone marks are altered in TFAP2-KO cells relative to in WT cells. These elements were associated with neither TFAP2-activated nor TFAP2-inhibited genes.  Table 2). We infer that at the first subset of TFAP2-dependent MITF peaks, TFAP2 is a pioneer factor (or participates in pioneering function), facilitating access to chromatin for MITF and other transcription factors (illustrated in Fig 4F"). Supporting this prediction, the transcription factor binding sites for MITF, SOX10, RUNX and FRA1 were strongly enriched at such elements (Fig 4F"'). At the second subset, TFAP2 is a , anti-H3K4Me3 CUT&RUN-seq (blue), anti-H3K27Ac CUT&RUN-seq (green) and RNA-seq (magenta) datasets at: (A) an TFAP2-pioneered-and-activated enhancer at the PAX3 (+60 kb) locus; (B) a non-pioneered TFAP2A-activated enhancer at the ENTPD6 (+26kb) locus; (C) an TFAP2-pioneered-and-inhibited enhancer at the ADAM19 (+21 kb) locus; and (D) a non-pioneered TFAP2A-inhibited enhancer at the FGF5 (+40kb) locus. Genotypes as labeled; y-axis in E applies to E-H, etc. (E-H") Violin plots conveying normalized reads of (E-H) anti-H3K27Ac (two independent replicates), (E'-H') ATAC-seq (four independent replicates) and (E"-H") anti-H3K4Me3 (two independent replicates) CUT&RUN-seq at E-E") TFAP2-pioneered-and-activated enhancers, (F-F") non-pioneered TFAP2-activated enhancers, (G-G") TFAP2-pioneered-andinhibited enhancers, and (H-H") non-pioneered TFAP2-inhibited enhancers. The number of peaks in each group is indicated. P-values shown were determined by Student's t-test. (I-L) Hypergeometric analysis of TFAP2 regulated enhancers at TFAP2-activated ( � ) and TFAP2-inhibited ( �� ) genes in WT cells (FDR < 0.05, |log2FC| > 1). ns; not significant. (M-O) Enrichment of transcription factor binding motifs at (M) TFAP2-pioneered-and-activated enhancers, at (N) non-pioneered TFAP2-activated enhancers and at (O) TFAP2-pioneered-and-inhibited enhancers as determined using HOMER motif analysis. P values were calculated using ZOOPS scoring (zero or one occurrence per sequence) coupled with hypergeometric enrichment analysis. TF; transcription factor. transcriptional activator that recruits MITF, also functioning as a transcriptional activator; we presume another protein serves as a pioneer factor at this subset (illustrated in Fig 4H"). Consistent with this notion, the binding site for JUN, a widely deployed pioneer factor [79], site is strongly enriched in these elements (Fig 4H"'). Examples are shown of TFAP2-activated MITF peaks near FRMD4B, CYP7B1, TRPM1, SOX9, EDNRB, MREG, GPR143, SNAI2, MEF2C, MYO5A, PAX3, EN1 and FOXI3 genes (Figs 4B-4D and S12D). At the third subset of TFAP2-dependent MITF peaks, where ATAC-seq signal was higher in TFAP2-KO cells than in WT cells (S14A-S14A ' Fig), TFAP2 may serve as a pioneer factor for MITF in MITF's proposed role as transcriptional repressor [17] (illustrated in S14A " Fig). However, this category of element was not enriched near genes inhibited by either TFAP2 or MITF (Hypergeometric test; p-value = 6.02 x 10 −02 and p-value = 9.12 x 10 −02 respectively). In summary, these results are consistent with TFAP2 facilitating access for MITF as a transcriptional activator to enhancers in both pioneer-factor and non-pioneer factor modes.
To test whether, reciprocally, TFAP2A binding depended on MITF, we carried out anti-TFAP2A CUT&RUN in MITF-KO cells. TFAP2A mRNA and protein levels were equivalent in MITF-KO and WT cells (S9E Fig), and the average TFAP2A peak height was globally equivalent as determined by CUT&RUN. At 13% (717/ 5334) of TFAP2-dependent MITF peaks, the TFAP2A peak was significantly reduced in MITF-KO cells (Fig 4A and 4C). At such loci, the average ATAC-seq signal was reduced in TFAP2-KO cells as compared to WT cells (Fig 4C and  4G-4G'). We termed these peaks mutually-dependent (i.e., mutually-activated peaks) (illustrated in Fig 4G"). Interestingly, mutually-dependent MITF/ TFAP2 peaks were enriched in binding motifs for TFAP2, MITF, BRM2 and TEAD4 but, unlike the other subsets of TFAP2-activated MITF peaks, not for SOXE (Fig 4H"'). SOX10 co-binds many loci together with MITF [15] and if SOX10 is absent from mutually-dependent peaks this may explain the dependence of TFAP2A binding on MITF at these sites. At~40% (288/ 717) of the mutually dependent peaks, and in the gene body surrounding these peaks, the repressive histone mark H3K27Me3, known to be applied by the polycomb complex [69,80], was significantly higher in MITF-KO cells relative to WT cells but, unexpectedly, not in TFAP2-KO cells relative to WT cells, even though MITF binding was lower in TFAP2-KO cells (illustrated in Figs 4G" and S12E-S12G).
In summary, the binding of MITF depends on TFAP2 at about one third of MITF peaks that overlap TFAP2A peaks. Such TFAP2-activated MITF binding occurs both at loci where nucleosome packing depends on TFAP2 (pioneer factor mode) and where it does not (nonpioneer factor mode). At a subset of the former, but none of the latter, TFAP2A binding is, reciprocally, MITF-dependent.

At a subset of MITF/TFAP2A co-bound peaks, TFAP2 paralogs inhibit chromatin access for MITF
In Fig 3 we established that at some TFAP2A peaks, TFAP2 paralogs maintain closed chromatin, and presumably inhibit binding of transcription factors. Consistent with this prediction, among MITF peaks overlapping TFAP2A peaks, 10% (1,605) of the MITF peaks were higher in TFAP2-KO cells than in WT cells (Figs 4A, 4E and S13), implying TFAP2 paralogs inhibit binding of MITF at this subset. At most peaks within this subset (58%, 924/1,605), the ATACseq signal was also significantly higher in TFAP2-KO cells versus WT cells (violin plots, Fig 4I, 4I' and illustrated in Fig 4I"), implying TFAP2 maintains condensed chromatin at this subset. DNA elements underlying this subset were enriched for transcription factor binding sites including those for SP1, NFY, JUN and TFE3 (Fig 4I"'). Moreover, these elements were modestly associated with TFAP2-inhibited genes ( Table 2; row 8).
Of note, at the majority of MITF peaks that overlap TFAP2A peaks (65%, 10,418/ 15,752), the height of the MITF peak was equivalent in TFAP2-KO and WT cells (TFAP2-independent MITF peaks) (S13 and S15A Figs), implying TFAP2 paralogs neither recruit nor repel MITF at these sites. Interestingly, such TFAP2-independent MITF peaks were not strongly enriched for the TFAP2 binding site (S15A Fig), implying that TFAP2 is attracted to many of these sites via other proteins rather than binding directly to the DNA. Such indirect binding may be less avid, as the average height of TFAP2-independent MITF peaks was smaller than that of TFAP2-dependent MITF peaks (S13 Fig, compare MITF signal in first and fourth cluster WT cells). As expected, TFAP2-independent MITF peaks were associated neither with TFAP2-activated nor TFAP2-inhibited genes ( Table 2; row 10).

TFAP2 and MITF co-regulate genes in the melanocyte gene regulatory network
The delayed pigmentation in zebrafish tfap2a/tfap2e double mutants, and the reduced expression of many pigmentation genes in TFAP2-KO cells, is consistent with two scenarios which are not exclusive of one another. In the first, TFAP2 paralogs directly activate MITF expression and thereby indirectly activate expression of pigmentation genes. In the second, TFAP2 paralogs directly activate expression of pigmentation genes. Supporting the first scenario, there is a TFAP2-pioneered-and-activated enhancer in intron 2 of the MITF gene (S16 Fig), and MITF mRNA levels are about 40-50% lower in TFAP2-KO cells than in WT cells. However, arguing against a strong role for this mechanism, many of the genes whose expression was most strongly reduced in MITF-KO cells compared to WT cells were completely TFAP2-independent, or were TFAP2-inhibited ( Table D in S1 Table). To assess the second scenario, we identified the set of genes activated directly by MITF, defined as MITF-activated genes associated with an MITF peak, and the set of genes directly activated by TFAP2, defined as TFAP2-activated genes associated with an TFAP2-activated enhancer (i.e., of pioneered or non-pioneered variety). Supporting the second mechanism, genes activated directly both by TFAP2 and by MITF were enriched for GO terms related to pigmentation and pigment cell differentiation (Figs 5A and 5B) [81]. Of note, genes related to pigmentation were also among those apparently directly regulated solely by MITF or TFAP2 paralogs (Fig 5C). We took a similar approach to identify genes directly inhibited by TFAP2 and/or by MITF (Fig 5D). Genes directly inhibited by both were strongly enriched for GO terms related to cell motility and cell migration (Fig 5D and 5E). Additional categories of genes strongly activated by TFAP2 included genes associated with cell proliferation and cell differentiation (Figs 6A, S17A and S17B) and genes strongly inhibited by TFAP2 included those associated with cell-cell adhesion (Figs 6B and S17C). We next examined the association of TFAP2-activated and TFAP2-inhibited genes with gene expression profiles from melanoma tumors and cell lines with distinct phenotypes [82][83][84][85][86][87][88]. Enrichment analysis showed that melanoma profiles previously found to be associated with high levels of MITF activity [89] were enriched for genes directly activated by TFAP2, including the subset associated with TFAP2-dependent MITF peaks (Fig 5F). Moreover, melanoma profiles associated with low levels of MITF activity were enriched for genes directly inhibited by TFAP2 (Fig 5F). These findings support other observations that the level of TFAP2 expression has a profound effect on the melanoma phenotype [32,90].
Finally, we performed cell proliferation, cell adhesion and cell migration assays on wildtype and TFAP2-depleted SK-MEL-28 cells. Consistent with effects on gene expression, TFAP2-KO cells showed reduced proliferation over 24, 48, 72 and 96 hours compared to WT cells (Fig 6C), and in cluster formation assays, TFAP2-KO cells showed increased cell-cell adhesion compared to WT cells (Fig 6D and 6E). However, in apparent contrast to the effects on gene expression (i.e., higher expression of cell migration genes in TFAP2-KO cells than in WT cells), while WT cells migrated to fill in a scratch within 24 hours, neither TFAP2-KO clone filled it within that time (S18A Fig). This finding also contrasts with the observation that the expression of tfap2e correlates negatively with the migratory capacity of zebrafish models of melanoma [32], but it is consistent with the accumulation of melanocytes in the dorsum of zebrafish tfap2a knockout embryos [23,63,64] and tfap2a/ tfap2e double mutant embryos. It is noteworthy that in MITF-KO cells, similar to in TFAP2-KO cells, cell migration genes are upregulated, but cell migration in-vitro is inhibited [17]. In summary, reduced expression of pigmentation, differentiation and cell growth genes and elevated expression of cell adhesion and cell invasion genes in TFAP2-KO cells as compared to WT cells is largely explained by the direct effects of TFAP2 paralogs on these genes.

DMTN
To test the TFAP2-pioneer-factor model, the two paralogs with highest expression, TFAP2A and TFAP2C, were deleted from the SK-MEL-28 cells. Creation and integration of several genomic datasets permitted us to identify genomic elements that were bound by TFAP2A in WT cells and that either lost or gained H3K27Ac signal in TFAP2-KO cells; we inferred that these elements were enhancers directly activated or inhibited by TFAP2. As expected by the TFAP2-as-pioneer-factor model, at a subset of TFAP2A-activated enhancers, TFAP2 paralogs were required to keep chromatin open, and at all TFAP2-inhibited enhancers they were required to keep it closed. We also showed that TFAP2 paralogs facilitated binding of another transcription factor (MITF), that this occurred at loci where TFAP2 facilitates chromatin access, and that this dependency was not reciprocal (at most loci bound by both transcription factors). These findings imply that TFAP2 paralogs recruit distinct transcription factors to distinct categories of enhancer, that these in turn recruit enzymes that open chromatin or that condense it (TFAP2 paralogs also may directly recruit such enzymes) [93], and that the state of chromatin accessibility determines access of other transcription factors. Therefore, TFAP2 paralogs facilitate chromatin access by other transcription factors, either as a pioneer factor or as a linker that connects pioneer factors to chromatin-remodeling machinery. Notably, the role of MITF at the co-bound enhancers is still largely unexplored.
Of note, there are loci where the pioneer factor FOXA1 appears to recruit proteins like GRG3 (groucho-related gene 3) that condense chromatin [77,78]. At TFAP2-inhibited enhancers, our data do not distinguish between whether TFAP2 binds open loci and closes them or instead binds to closed loci and helps stabilize the closed state. Of note, the majority of TFAP2A peaks did not coincide with enhancers either activated or inhibited by TFAP2 paralogs and were thus quiescent. In some cases, such apparently quiescent TFAP2A peaks may represent loci that TFAP2 paralogs have stably pioneered; a precedent for this scenario is that at a subset of elements pioneered by PAX7, chromatin remains open after the removal of PAX7 [94]. To explore this possibility, TFAP2A and ATAC-seq will need to be monitored in a developmental time series, or in TFAP2-KO cells before and after TFAP2A has been reexpressed.
Like other pioneer factors, TFAP2 can activate enhancers in a non-pioneer factor mode. At a subset of TFAP2A-activated enhancers the H3K27Ac signal, but not the ATAC-seq signal, was lower in TFAP2-KO cells than in WT cells. Evidence that such elements are indeed TFA-P2-activated enhancers is that their average H3K4Me3 signal was also lower in TFAP2-KO cells than in WT cells. We infer that TFAP2 activates these enhancers, but not as a pioneer factor. At such enhancers the continued presence of TFAP2 is necessary for continued acetylation of histone H3 lysine 27 (H3K27Ac), which fits with evidence that TFAP2 binds the histone acetyl transferase p300/CBP [95] and inhibits the NURD histone-deacetylase complex [93]. TFAP2 may attract other transcription factors without affecting nucleosome positioning; indeed, some TFAP2-dependent MITF peaks were found at non-pioneered TFAP2-activated enhancers. The fact that the TFAP2 binding site is not strongly enriched in non-pioneered TFAP2-activated enhancers implies TFAP2 binds these elements indirectly. Although we refer to these elements as non-pioneered TFAP2-activated enhancers, our experimental design could not rule out the possibility that some such enhancers were stably pioneered by TFAP2 as discussed above. However, the observation that these sites are less enriched for the TFAP2 binding site than the TFAP2-pioneered-and-activated enhancers, and that the binding sites of pioneer factors FOS and JUN are enriched [96], supports the alternative model that such elements are simply pioneered by different transcription factors.
Interestingly, at a subset of MITF/TFAP2A overlapping peaks where MITF was lost in TFAP2-KO cells, TFAP2A binding was lost in MITF-KO cells (MITF/TFAP2A mutually dependent peaks). There is precedent for reciprocal dependence of the binding of pioneer factors, in the cases of both FOXA1 and steroid hormone receptors, at subsets of sites where they are cobound [97]. Notably, at many MITF/TFAP2A mutually dependent peaks, the repressive mark H3K27Me3 accumulated in MITF-KO cells. This is consistent with evidence that the SWI/SNF complex, which MITF probably recruits to such loci, competes for access to chromatin against the H3K27Me3-depositing Polycomb repressor complex and that the binding of pioneer factors is impeded by condensed H3K27Me3-positive chromatin [80,98]. We propose that at the majority of TFAP2-dependent MITF peaks, but not at MITF/TFAP2A mutually dependent peaks, a measure of BRG1 binding is retained in MITF-KO cells, possibly recruited by another activator like SOX10, but this is not the case for MITF/TFAP2A mutually dependent peaks.
Finally, by identifying enhancers regulated by TFAP2, whether in pioneer mode or not, we were able to identify with high-confidence genes directly activated and directly inhibited by TFAP2. The genes directly activated by TFAP2, including by the specific mechanism of facilitating MITF access at a nearby enhancer, were enriched in GO terms related to pigmentation and proliferation, and overlapped significantly with expression profiles previously associated with high levels of MITF activity [82][83][84][85][86][87][88]. The MITF "rheostat model" predicts that high MITF activity promotes cell proliferation [2,12]. Proliferation of both MITF-KO [17] and TFAP2-KO cells (shown here) is slower than of the wild-type counterparts. We suggest there are two mechanisms to explain the shared expression profiles and phenotypes between TFAP2-KO cells and MITF-KO cells. First, the expression of MITF mRNA is about 40-50% lower in TFAP2-KO than WT cells. It is unclear how important this effect is because the average read depth of MITF peaks was equivalent in WT and TFAP2-KO cells, in particular at MITF-activated genes that are not also TFAP2-activated (e.g., MLANA). Second, TFAP2 facilitates binding of MITF to enhancers and promoters, leading to higher expression of associated MITFtarget genes. The "rheostat model" also predicts that low MITF activity promotes cell invasion [12], and there was a large overlap of TFAP2-inhibited and MITF-inhibited genes and these genes were enriched for those involved in cell-cell adhesion. We found evidence that TFAP2 directly inhibits enhancers, but this did not involve TFAP2 and MITF cooperativity. Directly-TFAP2-inhibited genes were associated with cell-cell adhesion and TFAP2-KO cells formed large clusters when cultured on low adhesive plates. Interestingly, in zebrafish melanoma, low tfap2e expression correlates with increased metastases [32]. We have previously shown that MITF binds to enhancers near MITF-inhibited genes [17]. However, mechanisms by which MITF inhibits gene expression in melanoma cells and in-vivo remain to be elucidated. Given that an MITF-low is a deadly status for melanoma, an interesting possibility for therapy would be to effectively covert them to MITF-high by manipulating pioneer factors such as TFAP2 paralogs or other MITF-cofactors. In summary, MITF activity in melanoma cells-and thus the phenotypes of these cells-depends in part on the presence of transcription factors that give MITF access to specific regulatory elements.

Ethics statement
D. rerio were maintained in the University of Iowa Animal Care Facility according to a standard protocol (protocol no. 6011616). All zebrafish experiments were performed in compliance with the ethical regulations of the Institutional Animal Care and Use Committee at the University of Iowa and in compliance with NIH guidelines. Zebrafish embryos were maintained at 28.5˚C and staged by hours or days post-fertilization (hpf or dpf).

Single cell suspension and fluorescence activated cell sorting
50 Tg(mitfa:GFP) transgeneic zebrafish embryos [99] were collected at 28 hours post fertilization (hpf) and washed in PBS without Ca2+ or Mg2+ (Life technology). Embryos were manually dechorionated under a dissecting scope (Leica KL300 LED), washed twice in 1X PBS and then cells were dissociated using a pestle and incubated in trypsin (0.25%)-EDTA (Life technology) at 33˚C for 30 minutes (ensuring to pipette mixed every 5 minutes). Reactions were quenched by adding PBS supplemented with 5% fetal bovine serum (FBS, Life technologies). Dissociated cells were centrifuged at 500 g for 5 minutes and cell pellets were washed with PBS-5%FBS before passing thought a Bel Art SP Scienceware Flowmi 40 μm cell strainers (ThermoFisher Scientific). Dissociated cells were re-suspended into single-cell solution and analyzed at the University of Iowa Flow Cytometry Facility, using an Aria Fusion instrument (Becton Dickinson, Franklin Lakes, NJ).

Single-cell RNA-sequencing (10x genomics) and data processing
The sorted cellular suspensions were loaded on a 10x Genomics Chromium instrument to gen- A custom reference genome for the Zebrafish was constructed with GRCz11 primary assembly (Ensembl) using Cell Ranger (Cell Ranger mkref function). 10X Genomics scRNA-seq reads were then processed and aligned to this reference with Cell Ranger (Cell Ranger Count function). Further analysis and visualization was performed using Seurat (v4.1.0) [51], cells with fewer than 200 RNA feature counts and greater than 5% mitochondrial contamination were removed with filtering. RNA counts were normalized, and Find-VariableFeatures was run with the following parameters: selection.method = "vst", nfeatures = 2,000. The cells were originally clustered in a UMAP by RunPCA then RunU-MAP with dims 1:30. FindNeighbors was run with dims 1:30 followed by FindClusters using a resolution of 1.2. Pigment clusters were further analyzed by sub-setting out and reclustering using dims 1:30, which identified 10 clusters. Pseudotime was performed using Monocle3 (v0.2.3.0) [100].

Generation of a zebrafish tfap2e loss-of-function allele
To generate the tfap2e loss-of-function allele, we designed paired (e.g., left and right) zinc finger nucleases (ZFN) targeting exon 2 of the tfap2e locus resulting in non-homologous endjoining and disruption of the open reading frame for Tfap2e. Briefly, the online tool, ZiFiT [101], was used to identify an optimal ZFN target site [utilizing the CoDa approach [102]. Once identified, a custom DNA fragment encoding the entire left or right zinc finger array (ZFA) along with flanking XbaI and BamHI restriction sites was synthesized (Integrated DNA Technologies, Coralville, IA). Subsequently, the ZFA fragment was subcloned into pMLM290 (Addgene, plasmid 21872), which includes a modified FokI nuclease domain [103]. Next, the fully assembled ZFN was PCR amplified, directionally cloned into pENTR-D/TOPO (Thermo-Fisher Scientific), and finally subcloned into pCS2+DEST using Gateway LR Clonase II enzyme mix (ThermoFisher Scientific). Once assembled, the final pCS2+ plasmids were sequence verified, linearized, mRNA synthesized in vitro (mMessage mMachine SP6 Kit, Ambion/ThermoFisher Scientific). Synthesized RNA was cleaned using the Qiagen RNeasy Kit (Qiagen) and both left and right ZFN components were co-injected into 1-cell stage zebrafish embryos. Following injections, embryos were initially screened via PCR and restriction enzyme digest to confirm editing at the target site. Upon confirmation, additional embryos from a similar clutch (F0's) were allowed to develop into adulthood, 'mosaics' identified and out-crossed, and a stable F1 generation isolated.

Cell lines, reagents, and antibodies
The cells referred to as WT throughout the document are the parent SK-MEL-28 (HTB-72) line. They and the derivative line, delta6-MITF knockout cells (referred to as MITF-KO cells in this work), were obtained from the laboratory of Dr. Eirikur Steingrimsson. The cells were grown in RPMI 1640 medium (Gibco #5240025) supplemented with 10% FBS (Gibco #10270106) at 5% CO 2 and 37˚C. Cells were tested for, and determined to be free of, mycoplasma. SK-MEL-28 cells harbor the BRAF V600E and p53 L145R mutations [104].

Generation of TFAP2A; TFAP2C knockout cell lines (TFAP2-KO)
TFAP2-KO clones were generated using the Alt-R CRISPR-Cas9 technology from Integrated DNA Technologies (IDT). Briefly, crRNAs targeting exon 2 of TFAP2A and TFAP2C were designed using the Cas9 guide RNA design checker tool (crRNA sequences below). Equimolar concentrations of crRNA and tracrRNA (IDT, #1072532) were annealed to form gRNA complexes. The ribonucleoprotein (RNP) complex was prepared by mixing gRNAs and Cas9 protein (IDT #1081058). SK-MEL-28 cells were transfected with constructs encoding components of RNP complexes using the Lipofectamine CRISPRmax Cas9 transfection reagent (Thermo-Fisher #CMAX00015) following the manufacturer's protocol. Single-cell colonies were screened by PCR and Sanger sequencing using primers flanking the cut sites (primer sequences below). Mutant clones (clone 2.12 and clone 4.3) were selected and further screened by western blotting, using anti-TFAP2A and anti-TFAP2C antibodies. The control cell lines used in this study were generated following this protocol but without adding gRNA duplexes.

SDS-PAGE and Western blotting
TFAP2-KO and WT cells were washed in ice-cold PBS. RIPA buffer containing protease inhibitors (Roche, cOmplete Mini) was added and cells were lysed on ice for 20 minutes. Cell lysates were centrifuged at 14,000 g for 20 minutes and the quantity of protein in the supernatants was quantified using Bradford assays (Bio-Rad #5000002). Laemmli sample buffer (Bio-Rad #1610747, 5% 2-mercaptoethanol) was added to 20 μg protein and samples were boiled at 95˚C for 5 minutes before being loaded onto a 10% SDS-polyacrylamide gel (Bio-Rad #4568034). Protein was transferred to polyvinylidene fluoride (PVDF) membranes (Thermo Scientific #88520), which were incubated overnight with primary antibody. Membranes were washed 3 times with TBS-T and incubated with horseradish peroxidase-conjugated anti-rabbit or anti-mouse for 1 hour at room temperature, washed, and imaged using an Amersham Imager 600.

ATAC-seq
ATAC-seq was performed according to [105,106] with minor alterations. Briefly, 70,000 TFAP2-KO cells (clone 2.12 and clone 4.3, four replicates each) and WT cells (four replicates) were lysed in ice-cold lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% NP-40: Sigma). Transposition was performed directly on nuclei using 25 μl tagmentation reaction mix (Tagment DNA Buffer #15027866, Tagment DNA Enzyme #15027865 from Illumina Tagment DNA kit #20034210). Tagged DNA was subjected to PCR amplification and library indexing, using the NEBNext High-Fidelity 2x PCR Master Mix (New England Biolabs #M0451S) with Nextera DNA CD Indexes (Illumina #20015882), according to the following program: 72˚C for 5 minutes; 98˚C for 30 seconds; 12 cycles of 98˚C for 10 seconds, 63˚C for 30 seconds, and 72˚C for 1 minute. The PCR product was purified with 1.8 times the volume of Ampure XP beads (Beckman Coulter #A63881). Library quality was assessed using a BioAnalyzer 2100 High Sensitivity DNA Chip (Agilent Technologies). All DNA libraries that exhibited a nucleosome pattern were pooled and processed for 150bp paired-end sequencing.

ATAC-seq peak calling and differential analysis
ATAC-seq was performed using 150 bp paired-end sequencing reads. Raw ATAC-seq reads were trimmed using Trim Galore Version 0.6.3 (Developed by Felix Krueger at the Babraham Institute) and aligned to human genome assembly hg19 (GRCh37) using Bowtie 2 [107,108] with default parameters. Sorting, removal of PCR duplicates, and identification of fragments shorter than 100 bp as the nucleosome-depleted-regions (NDRs), was performed using BAM filter version 0.5.9. DeepTools version 3.3.0 [109] was used to check the reproducibility of the biological replicates and generate bigWig coverage files for visualization. Peaks were called using model-based analysis of ChIP-seq 2 (MACS2, version 2.1.1.20160309.6) [110]. NDRs for which accessibility differed between TFAP2-KO and WT cells were identified using DiffBind version 2.10 [111] with log2 fold-change threshold of >1 and a false discovery rate (FDR) < 0.05. NDRs that are directly regulated by TFAP2 were identified by overlapping differentially accessible NDRs with anti-TFAP2A CUT&RUN peaks; a 1-bp window was used to define overlap. Peaks were assigned to genes using GREAT with a peak-to-gene association rule of the nearest-gene-within-100 kb [112]. Both the raw ATAC-seq files and processed sequencing data presented in this manuscript have been deposited in the Gene Expression Omnibus (GEO) repository (GSE number pending).

CUT&RUN library preparation and data analysis
CUT&RUN libraries were prepared using the KAPA Hyper Prep Kit (Roche). Quality control post-library amplification consisted of fragment analysis using the 2100 Bioanalyzer (Agilent).
Libraries were pooled, brought to equimolar concentrations, and sequenced with 150 bp paired-end reads on an Illumina HiSeq X platform (Novogene, Sacramento, CA). For quality control, paired-end FASTQ files were processed using FastQC (Babraham Bioinformatics). Reads were trimmed using Trim Galore Version 0.6.3 (Developed by Felix Krueger at the Babraham Institute) and then mapped against the hg19 genome assembly using Bowtie2 version 2.1.0 [107,108]. The mapping parameters and peak calling of MACS2 peaks [110] were performed as previously described [113,114] against their matching control IgG samples. Differential analysis of H3K27Ac and of H3K27Me3 signal in WT and TFAP2-KO cells was preformed using MACS2 with a Log2 fold-change threshold of 1, and p-value < 1 x 10 −5 . Differential H3K4Me3, MITF and TFAP2A signal in WT, TFAP2-KO and when mentioned MITF-KO cells was determined using DiffBind version 2.10.0 [111] with a Log2 fold-change threshold of 1, and FDR < 0.05. DeepTools version 3.3.0 was used to check the reproducibility of the biological replicates, to generate bigwig normalized (RPKM) coverage files for visualization and to plot average CUT&RUN-seq and ATAC-seq profiles (-plotProfile) and generate heatmaps (-plotHeatmap) of normalized reads [109]. MultiBigwigSummary was used to extract read counts (-outRawCounts) [109] and Prism was used to generate Violin and Box plots. Peaks were assigned to genes using GREAT with a peak-to-gene association rule of the nearest-gene-within-100 kb [112]

RNA-seq
Four replicate RNA-seq experiments were performed on TFAP2-KO cells (clone 2.12 and clone 4.3) and WT cells. Total RNA was extracted by direct cell lysis using the RNeasy Plus Mini Kit with QiaShredder (Qiagen #47134). RNA samples with an RNA integrity number (RIN) above nine were used for library generation and 150 bp paired-end sequencing on the Illumina HiSeq2500 platform (Novogene, Sacramento, CA). FASTQ sequence files were processed using FastQC (Babraham Bioinformatics) for quality control, and reads were trimmed using Trim Galore Version 0.6.3 (Developed by Felix Krueger at the Babraham Institute) and subsequently aligned to human genome assembly hg19 (GRCh37) using STAR [115]. The output of the-quantMode GeneCounts function of STAR was used for the calculation of differential transcript expression using DESeq2 [116].The rlog function was used to generate log2-transformed normalized counts. Adjusted p-value < 0.05 was used as the threshold for statistically significant differences. Functional enrichment analyses was performed using PAN-THER [117]. A full list of genes differentially expressed between TFAP2-KO and WT cells is provided (Tables B-C in S1 Table).

Motif analyses
Both de novo and known motifs were identified within 200 bp of TFAP2-activated and TFA-P2-inhibited enhancer and promoter peak summits using HOMER (-findMotifsGenome).

Statistical analysis
Fisher's Exact Test was used to assess association between TFAP2-regulated elements (enhancers and promoters) with TFAP2-regulated gene expression. Fisher's Exact Test is a statistical significance test based on the hypergeometric distribution and is used to analyze contingency tables to determine if there is an association between the classes in the columns and rows of the contingency table. Use of this test is relevant for the analysis of differential expression because we classify enhancers as either TFAP2-activated (downregulated in TFAP2-KO cells) or TFAP2-inhbited (upregulated in TFAP2-KO cells) as well as near a target of interest (i.e., TFAP2-activated or TFAP2-inhibited) or not. This test determines if there is an association between genes that are near an TFAP2-regulated enhancer and genes whose expression varies between TFAP2-KO cells and WT cells. Data processing and analysis was performed in R and the code can be found at https:// GitHub.com/ahelv/Differential_Expression. GraphPad-Prism was used to perform Student's t-test as indicated in the figure legends.

Proliferation assay
Cell proliferation was monitored at 24, 48, 72 and 96 hours using the CellTiter 96 Aqueous One Solution reagent kit according to the manufacturer's protocol. Cells were serum depleted for 24 hours before transferring 5,000 cells per well of cell culture treated 96-well plate (corning) containing RPMI media supplemented with 10% FBS and 1% PS. On the day of the experiment, CellTiter 96 Aqueous One Solution was added to cells at a ratio of 6:1 (Media: solution). After 1 hour at 37˚C in a humidified, 5% CO 2 atmosphere, the absorbance at 450nm was recorded using an ELISA plate reader. Background absorbance was measured from cell free wells containing media and CellTiter 96 Aqueous One Solution. The experiment was performed in triplicate and repeated 3 times.

Cluster formation assay
Cluster formation assays were performed as previously described [32]. In brief, cells were trypsinized and resuspended in RPMI media supplemented with 10% FBS and 1% PS. A flat bottom Low Attachment Surface-well plate (Corning) was seeded with 5,000 cells/well in a final volume of 100uL. Clusters were allowed to form over the course of 3 days in a humidified 37˚C incubator at which time the number of clusters were counted per well using Image J software. The experiment was performed in triplicate and repeated twelve times.

Wound scratch assay
A total of 500K cells were seeded per well of 6-well plate (Thermo Scientific, # 1483210) to reach confluent monolayer. Cells were incubated in serum free media for 6 hours before wounding with a 200 μL pipette tip. Scratches were manually imaged on an inverted light microscope (Leica #10445930) every 6 hours over a 24-hour time period. The distance of scratch closure between WT and TFAP2-KO cells were analyzed with Image J software.