A regulatory sub-circuit downstream of Wnt signaling controls developmental transitions in neural crest formation

The process of cell fate commitment involves sequential changes in the gene expression profiles of embryonic progenitors. This is exemplified in the development of the neural crest, a migratory stem cell population derived from the ectoderm of vertebrate embryos. During neural crest formation, cells transition through distinct transcriptional states in a stepwise manner. The mechanisms underpinning these shifts in cell identity are still poorly understood. Here we employ enhancer analysis to identify a genetic sub-circuit that controls developmental transitions in the nascent neural crest. This sub-circuit links Wnt target genes in an incoherent feedforward loop that controls the sequential activation of genes in the neural crest lineage. By examining the cis-regulatory apparatus of Wnt effector gene AXUD1, we found that multipotency factor SP5 directly promotes neural plate border identity, while inhibiting premature expression of specification genes. Our results highlight the importance of repressive interactions in the neural crest gene regulatory network and illustrate how genes activated by the same upstream signal become temporally segregated during progressive fate restriction.


Introduction
The neural crest is a migratory stem cell population that gives rise to multiple components of the vertebrate body plan, including the peripheral nervous system and the craniofacial skeleton [1]. The formation of neural crest cells is controlled by a gene regulatory network (GRN) composed of signaling systems, transcription factors, and other regulatory proteins [2][3][4]. Hence, the process of neural crest formation can be subdivided into a series of distinct regulatory modules that flow logically from one to the next. First, extracellular signaling molecules (WNTs, FGFs, BMPs and Notch) present in the early embryo initiate the induction of the neural plate border. The sequential activation of these signal transduction pathways promotes the expression of transcription factors termed neural plate border specifier genes (e.g., PAX7, MSX1, ZIC1, TFAP2A) [5]. These factors cooperate with signaling systems to activate transcription of neural crest specifier genes, which include neural crest markers like SNAI2, FOXD3, and SOX10. Neural crest specifier genes, in turn, activate the effector genes that endow these cells with their unique properties, such as multipotency and migratory capabilities [4].
The neural crest GRN is one of the most comprehensive regulatory networks described for a vertebrate cell type, and has served as a platform for investigating molecular mechanisms involved in cell fate commitment and vertebrate evolution [6][7][8]. Yet, we still have a superficial understanding of the molecular logic encoded in the architecture of the network. Metazoan GRNs are structured as a collection of sub-circuits, in which small groups of genes are linked together by positive or negative interactions [9,10]. GRN sub-circuits are the functional units of the program, performing tasks required to change transcriptional states or to establish spatial domains of gene expression [11]. While past functional and biochemical studies resulted in a substantial expansion of the neural crest GRN, only a few examples of regulatory sub-circuits have been described in this cell type [12]. Delineation of additional regulatory motifs is necessary to understand how this genetic program orchestrates developmental transitions.
Recently, we have identified Wnt target gene AXUD1 (also known as CSRNP-1) as a critical component of the neural crest GRN [13]. This transcription factor forms a protein complex with PAX7 and MSX1 to directly promote the expression of specifier genes. Notably, while PAX7 and MSX1 are present during gastrula stages [14,15], specifier genes including FOXD3 and SOXE factors are only expressed following the start of AXUD1 expression at Hamburger and Hamilton stage 7 (HH7) [13]. This indicates that AXUD1 plays an essential role in activating the neural crest specification program. However, direct activators of this transcription factor have yet to be identified. While AXUD1 is hypothesized to be downstream of canonical Wnts, it does not accurately recapitulate the activity of this signaling system, as Wnt ligands are already present in the early gastrula [16], but AXUD1 transcription only starts during neurulation. This discrepancy suggests a complex mode of regulation that may involve integration of numerous inputs. Identifying how upstream factors control AXUD1 expression may reveal new mechanisms responsible for the emergence of neural crest identity.
Here we describe a regulatory sub-circuit that controls the onset of neural crest formation in vertebrate embryos. By surveying the AXUD1 locus with chromatin conformation capture, we identified a tissue-specific enhancer that interacts with the promoter of this gene in neural crest cells. Dissection of this enhancer allowed for the delineation of a sub-circuit that operates downstream of Wnt signaling to control the temporal expression of neural crest genes. AXUD1 is activated by Wnt signaling and neural plate border genes TFAP2A, MSX1, and ZIC1. Crucially, direct repression by Wnt effector gene SP5 restricts the timing of AXUD1 expression to prevent premature expression of specification genes. We assembled these interactions in a sub-circuit that links Wnt-target genes in an incoherent feedforward circuit, which controls the sequential activation of genes in the neural crest lineage. Our results highlight the importance of repression in the neural crest GRN and demonstrate how the interactions between Wnt target genes drive progressive fate restriction.

Tissue-specific regulatory elements interact with the AXUD1 promoter
AXUD1 is first detected in the neural crest lineage at (HH7), before the onset of specification markers FOXD3, ETS1 and SOX9 [13]. In situ hybridization in avian embryos shows robust expression of this transcription factor at HH8-9 (Fig 1A and 1B). To identify the cis-regulatory elements that control AXUD1 expression, we first employed Chromatin Conformation Capture (3C) in neural crest cells isolated from embryos at these developmental stages. Dorsal neural folds from the cranial region of HH8-9 embryos were micro-dissected with iridectomy scissors and dissociated into a single-cell suspension. Cells were fixed with formaldehyde, and the crosslinked chromatin was digested with restriction enzyme NCOI ( Fig 1C). After proximity ligation, we used RT-PCR to estimate the frequency of interaction between the promoter and sectors of the AXUD1 locus (defined by the presence of NCOI sites, see Methods). For this analysis, we profiled a region spanning~100kb between the two genes that flank AXUD1. The resulting 3C plot (Fig 1D) contained the expected high interaction values near the promoter, but also uncovered a large interaction peak~55kb upstream of the transcription start site (TSS). This interaction profile was consistent in all biological replicates of the experiment (n = 3).
To identify active enhancers of AXUD1, we characterized genomic regions that interact with its promoter using datasets of chromatin accessibility (ATAC-seq), H3K27Ac association and TFAP2A occupancy in neural crest cells [17]. TFAP2A is a pioneer transcription factor that marks active neural crest enhancers [18]. This analysis led to the identification of 11 putative cis-regulatory elements ( Fig 1E). To test if these elements were active, we performed RT-PCR for enhancer RNA (eRNA) in dissected neural folds of HH9 chicken embryos. We detected high levels of transcription both in TFAP2A-occupied putative enhancers adjacent to the 3C peak (elements 1-4) and in elements adjacent to the AXUD1 promoter (elements 9-10) ( Fig 1F). To test if these regions were able to drive gene expression, we conducted transgenic reporter assays in chick embryos. Putative enhancers (ax.1-ax.6, Fig 1D and 1E) were cloned in the pTK-eGFP vector [19] (Fig 1G) and electroporated in the ectoderm of gastrula-stage embryos (HH4). Of these elements, only ax.2 was active in neural crest cells at HH9 (Fig 1H), whereas the other constructs tested (S1 Table) were unable to drive reporter expression ( Fig  1I). Activity of ax.2 was consistent with previous testing of putative neural crest enhancers based on genome accessibility [20]. Histological analysis confirmed that ax.2 was specifically active in neural crest cells (Figs 1J and S1), and we thus refer to it as AXUD1 enhancer 1 (Axud1E1).
In order to establish the temporal pattern of Axud1E1 activity, we next examined transgenic embryos at different stages of neural crest development. We detected specific signal from the Axud1E1 reporter construct in both pre-migratory (Fig 2A) and migratory neural crest ( Fig  2B-2D). Double fluorescent in situ hybridization (FISH) for AXUD1 and eGFP in transgenic embryos showed that prior to migration, the enhancer drives reporter expression only in Chromatin conformation capture identifies a tissue-specific enhancer that interacts with the AXUD1 promoter. (A-B) Whole mount in situ hybridization for AXUD1, depicting the specific mRNA expression during neural crest specification stages HH8 (A) and HH9 (B). (C) Schematic representation of a Chromosome Conformation Capture-qPCR (3C-qPCR) experiment. Crosslinked neural crest cells were incubated with restriction enzyme NCOI and a DNA ligase. These steps allow the formation of hybrid DNA molecules combining restriction fragments that were in close proximity in nuclei of the cells. Primers spanning the AXUD1 locus were paired with a primer anchored in the AXUD1 promoter to amplify hybrid DNA junctions and quantify the interaction frequency with the promoter. (D-F) Identification of active enhancers in the AXUD1 locus. 3C-qPCR interaction map for the AXUD1 locus reveals regions of high interaction frequency with the promoter region (blue dotted line, D). Gray dotted lines highlight the six elements tested in transient transgenesis assays (see below). Error bars represent ± SEM. Purple and blue lines in (D) represent two replicates of the same 3C experiment. ATAC-seq, H3K27ac and TFAP2A CUT&RUN profiles at AXUD1 locus depict regions of accessibility and active chromatin regions (E). (F) eRNA quantification (RT-PCR, normalized to reference gene) for the regions numbered in (E) indicates the level of transcription in the promoter region and putative distal regulatory elements. Error bars represent ± SEM. Statistical significance determined via an unpaired t-test. (G) Reporter vector used in transgenesis reporter assays. The construct consists of the candidate enhancer region cloned upstream of the HSV-tk minimal promoter driving eGFP expression.

PLOS GENETICS
neural crest cells that express AXUD1 (Fig 2E and 2F). Since AXUD1 is rapidly downregulated after neural crest cells delaminate and begin to migrate [13], we examined transcription of the reporter construct at later developmental stages. FISH for eGFP at HH11 showed that late migratory neural crest cells display no transcription of the reporter gene (S2A Fig). This indicates that the labeling of this cell population with eGFP ( Fig 2D) is likely due to the stability of eGFP protein [21]. These results show that Axud1E1 activity accurately recapitulates the spatial and temporal expression of the endogenous gene during early neural crest development.
To understand the regulation of Axud1E1, we performed a series of deletions in the original reporter construct (Fig 2G). We found that the original 2.5kb fragment contained two smaller

PLOS GENETICS
elements capable of driving reporter activity in neural crest cells (Fig 2H and 2I). The first element (Axud1E-500, located in the 3' region of the original construct) spanned 500bp, and was highly active in AXUD1+ cells (Figs 2H and S2B). Chromatin immunoprecipitation (ChIP) confirmed that Axud1E-500 was associated with H3k27Ac and TFAP2A (S2C Fig). The second element (Axud1E1-SE), located at the 5' region of Axud1E1, displayed similar but weaker activity (Figs 2I and S2D). Since both elements are active in the same cells (S2E and S2F Fig) and Axud1E1-SE is located in a region of closed chromatin (Fig 2G), we hypothesized that the latter acts as a shadow enhancer [22]. To test this, we first examined the activity of both cis-regulatory elements in FACS sorted neural crest cells. Quantification of eRNA with qPCR showed specific activity of Axud1E-500, while Axud1E-SE was transcriptionally silent (S2G Fig). Next, we assessed the existence of a compensatory mechanism between the two enhancers by disrupting Axud1E-500 with CRISPR/Cas9. We electroporated the right side of chick embryos with a CAS9 expression vector and gRNAs targeted to the Axud1E-500 [23], whereas the left side was transfected with CAS9 and a control gRNA. eRNA quantification with qPCR showed that targeting of Axud1E1-500 with two different gRNAs results in a significant activation of Axud1E-SE ( Fig 2J) in neural crest cells. These results show that AXUD1 is regulated by a robust tissue-specific enhancer (Axud1E-500) that is in close proximity to a shadow regulatory element.

Wnt effectors LEF1/CTNNB1 directly regulate AXUD1
Our previous studies identified AXUD1 as a downstream target of canonical Wnt signaling, acting as a secondary effector of the pathway during neural crest formation [13]. Thus, we tested whether Wnts directly regulate this transcription factor via Axud1E1-500. To define if the enhancer responds to Wnt signaling, we disrupted the pathway in transgenic avian embryos transfected with the Axud1E1-500 reporter construct. For these experiments, we transfected embryos with a control morpholino on the left side, and morpholinos targeting components of the Wnt signaling pathway on the right side ( Fig 3A). We found that knockdown of either the WNT1/4 ligands ( Fig 3B) or CTNNB1 (Fig 3C) resulted in a robust loss of enhancer activity (Fig 3D), suggesting Wnt signaling directly regulates Axud1E1. To further characterize this interaction, we next identified which Wnt nuclear effector cooperates with CTNNB1 to mediate the pathway's activity in neural crest cells. We examined the expression of the three TCF/LEF paralogs that promote gene expression downstream of Wnt signaling (TCF7, TCF7L2 and LEF1). In situ hybridization showed that LEF1 is strongly enriched in the neural crest during specification, whereas we were unable to detect discernable expression patterns for TCF7 and TCF7L2 in the head of avian embryos (Fig 3E and 3G). Immunohistochemistry for LEF1 was consistent with the in situ hybridization analysis, showing robust and specific expression of the protein in the neural crest ( Fig 3H, S3 Fig).
These results suggest that LEF1 is the main effector of Wnt signaling in neural crest cells. They also indicate that the tissue-specific regulation of Wnt nuclear effectors may play an important role in the activation of target genes. To examine the regulation of Axud1E1-500 by Wnt nuclear effectors, we conducted a series of ChIP experiments in dissected HH9 dorsal neural folds. While pulldown of LEF1 resulted in co-immunoprecipitation of Axud1E1-500, the other TCF/LEF paralogs (TCF7, TCF7L1 or TCF7L2) failed to interact with the enhancer (Fig 3I). Next, we examined whether the interaction between LEF1 and Axud1E1-500 was stage-specific and dependent on the activation of Wnt signaling. We found that ChIP performed on lateral ectodermal tissue at gastrula stages (HH4) or on neural folds (HH9) transfected with a Wnt dominant-negative construct [16] showed no interaction between LEF1 and the enhancer (Fig 3J). Furthermore, the shadow enhancer Axud1E1-SE, which is not accessible  Fig 3J). Finally, we also found a robust association between CTNNB1 and Axud1E1-500 ( Fig 3K). These results indicate that LEF1 regulates AXUD1 by interacting with a tissue-specific enhancer and is the main effector of canonical Wnt signaling in neural crest cells.

Axud1E1 integrates multiple upstream inputs to regulate gene expression
Enhancer analysis is a powerful tool for the identification of novel GRN interactions. To elaborate on the genetic circuit that controls AXUD1 expression, we further dissected the Axud1E1-500 enhancer. We divided this element into five 100bp segments (A-E), which were mutated by replacing the endogenous sequence with a 100bp segment of the coding sequence of eGFP ( Fig 4A) [24]. We found that mutations in segments A and B did not disrupt enhancer activity in neural crest cells (Fig 4B and 4C). In contrast, mutations in segments C, D and E resulted in a strong loss of reporter expression (Fig 4D and 4F). We next quantified the effects of each mutation using flow cytometry. Embryos were co-electroporated with the original 2.5kb Axud1E1 enhancer driving expression of mCherry, and each of the mutant constructs upstream of eGFP (Fig 4G and 4H). We then dissociated transgenic embryos and quantified eGFP expression in mCherry-positive cells. We observed that the A and B mutants were still active in neural crest cells, despite a reduction in eGFP expression. Activity of the C, D and E mutants was indistinguishable from the negative control vector (Fig 4I).
This analysis suggests that segments C,D, and E of Axud1E1-500 are critical for enhancer activity. Consistent with this, a reporter construct containing the simultaneous mutations in A and B is active in neural crest cells (Figs 4J and S4A). We thus generated a 300bp version of the enhancer containing only segments C, D and E. This construct (Axud1E1-300) (Fig 4K and  S4B Fig), was able to drive specific eGFP expression ( Fig 4L). Next, we screened this 300bp enhancer core using the JASPAR database to identify binding sites of potential upstream regulators. We found motifs for genes expressed in the neural crest lineage, like MSX1 and ZIC1, and for the Wnt target gene and multipotency factor SP5 ( Fig 4M). To assess the importance of these motifs, we generated versions of the enhancer carrying mutations in these binding sites (Fig 5A and S2 Table) for transgenic reporter assays. Enhancer variants were transfected on the right side of chicken embryos, whereas the wild-type enhancer was transfected on the left side. While mutation of the MSX1 and ZIC1 motifs resulted in loss of enhancer activity (Fig 5B, 5C and 5E), loss of SP5 binding sites drastically increased enhancer output (Fig 5D  and 5E). These findings suggest that MSX1 and ZIC1 are activators of Axud1E1 while SP5 acts as a repressor of the enhancer.
To confirm the roles of neural plate border genes (MSX1 and ZIC1) and the multipotency factor SP5 in the regulation of AXUD1, we performed loss-of-function experiments for these transcription factors. Gastrula-stage embryos (HH4) were electroporated with control and targeted morpholinos and incubated until HH9. MSX1 and ZIC1 morphants exhibited a reduction in reporter expression (Fig 5F and 5G). In contrast, knockdown of SP5 resulted in an increase in enhancer activity (Fig 5H and 5I). Furthermore, the SP5 morphant side of the embryos displayed a robust upregulation of the endogenous Axud1E1 ( Fig 5J) and premature transcription of AXUD1 at HH6 (Fig 5K). To confirm that SP5 acts as a repressor of AXUD1, we overexpressed SP5 using pCI-H2B-RFP [25]. Analysis of dorsal neural folds from single embryos (left vs. right) with RT-PCR and in situ hybridization showed a significant reduction in the expression of AXUD1 (Fig 5L and 5M) following SP5 overexpression.
Finally, we examined the regulation of Axud1E1 accessibility. Axud1E1 is associated with TFAP2A (Figs 1E and S2C), a pioneer transcription factor involved in chromatin remodeling during neural crest specification [17]. Despite this, we found that mutation of several TFAP2A binding sites in the Axud1E1 did not affect the activity of the enhancer in reporter assays (S5C Fig). This led us to postulate that TFAP2A regulates Axud1E1 at the chromatin level. Consistent with this, we found that knockdown of the pioneer factor caused a significant reduction in accessibility of the endogenous enhancer (S5D Fig), as assayed by ATAC-qPCR. Taken together, these data indicate that AXUD1 is regulated by a combination of activating and inhibitory inputs, where TFAP2A promotes chromatin accessibility, while MSX1 and ZIC1 act as bona fide transcriptional inputs. In contrast, SP5, which is also directly regulated by canonical Wnt signaling [26], represses AXUD1 via Axud1E1.

SP5 directly regulates components of the neural crest GRN
Our dissection of Axud1E1 allowed for the delineation of a regulatory circuit that controls AXUD1 expression, in which Wnt signaling and neural plate border specifier genes act as positive inputs, and SP5 acts as a repressor. To determine how SP5 regulates AXUD1 during neural crest specification, we examined the expression patterns of these genes using in situ hybridization. We found that SP5 was robustly expressed in the neural plate border at HH6 (Figs 6A and S6A), just prior to the onset of AXUD1 expression. At specification stages, SP5 expression is rapidly downregulated in the neural crest lineage (Figs 6B and 6C and S6B-S6D), just as AXUD1 transcription begins. At this developmental stage, SP5 transcripts are detected in the non-neural ectoderm (S6C ' Fig), whereas AXUD1 is restricted to the neural crest ( Fig 6D). Consistent with this, double FISH shows co-localization of SP5 and neural plate border PAX7 (S6E Fig), whereas we found no overlap between the former and neural crest marker TFAP2B at migration stages (S6F and S6G Fig). The temporal segregation in the expression of AXUD1 and SP5 is also evident in their regulation by canonical Wnt signaling. ChIP for CTNNB1 shows that the binding of the nuclear effector to Axud1E1 significantly increases from HH5 to HH9, whereas association to a SP5 enhancer follows the opposite trend (S7A Fig). To confirm that SP5 directly regulates Axud1E1, and to gain further insight into its function, we mapped the occupancy of this transcription factor using Cleavage Under Targets and Release Using Nuclease (CUT&RUN) [27]. This technique relies on a protein A-Mnase fusion protein, which binds to an antibody of choice and cleaves DNA fragments associated with the targeted transcription factor (Fig 6E). CUT&RUN analysis in neural folds dissected from HH7 embryos revealed 1016 SP5-bound regions, the majority of which were located in intergenic regions (87%) (Fig 6F). Motif enrichment analysis of these peaks identified a GC-rich motif common to Sp/Klf zinc-finger transcription factors [28] (pValue = 1e-155) (S7B Fig). SP5 peaks were above background levels, and positive for H3K27Ac (Fig 6G and 6H), indicating that SP5 interacts with active enhancers. Consistent with our enhancer dissection (Fig 5), we also observed robust association of the transcription factor with Axud1E1 ( Fig 6I).
The results from our genomic analysis indicate SP5 regulates a large number of genes in the neural crest lineage. Gene ontology analysis shows SP5 targets are involved in many developmental processes like nervous system development, neurogenesis, and morphogenesis ( Fig 6J). Our results are also consistent with reports that SP5 regulates its own transcription, as we detect strong occupancy in the SP5 gene promoter (Fig 6K). Other neural crest genes, including SOX9 (Fig 6K), TFAP2A, ZIC1, and MSX1 also contained SP5 peaks in active enhancers in their loci (S7D Fig). Taken together, our results allow for the assembly of a regulatory sub-circuit that controls sequential activation of Wnt target genes during neural crest development (Fig 6L).

Discussion
Cis-regulatory analysis has served as a powerful tool to decode the genetic programs that control metazoan development [29]. Dissection of tissue-specific elements can uncover how the MSX1, ZIC1 and SP5 binding motifs defined using JASPAR database. Kb, kilobase; HH, Hamburger and

PLOS GENETICS
multiple inputs are processed to generate specific domains of gene activity. Previous studies on neural crest enhancers have resulted in the identification of new genes involved in the formation of this cell type [24,25,[30][31][32]. Tissue-specific enhancers have also been useful for labeling and isolating pure populations of neural crest cells for genomic analysis [7,33,34]. Whereas previous studies have employed evolutionary conservation, histone marks or accessibility, our results indicate that chromosome conformation capture is a viable strategy for rapid identification of cis-regulatory elements. The subsequent dissection of Axud1E1 revealed an intricate combination of upstream inputs, which included effectors of signaling systems, multiple neural plate border genes and both positive and inhibitory inputs. The complexity of AXUD1 regulation was further underscored by the existence of a shadow enhancer adjacent to the core of Axud1E1. These results support our previous suggestion that AXUD1 occupies an important position in the neural crest GRN and is a direct target of Wnt signaling [13].
Comparative analyses of developmental GRNs show that these genetic programs are composed of collections of sub-circuits, which are responsible for executing well-defined tasks [11]. Here we outline a sub-circuit centered on Wnt signaling and effector genes SP5 and AXUD1, which supports the sequential activation of regulatory modules in the neural crest GRN. The sub-circuit ( Fig 6L) is structured as an incoherent feedforward loop [35], in which Wnts activate both SP5 [26,36] and AXUD1, while SP5 represses AXUD1 by interacting with a tissue-specific enhancer. These types of motifs are termed incoherent since the upstream signal (Wnts) directly activates a target gene (AXUD1) but also promotes the expression of a repressor (SP5) of said gene. We propose that the sub-circuit relies on repressive interactions between Wnt target genes to promote temporally distinct phases of gene expression. The inhibitory interaction between SP5 and AXUD1 prevents these factors and their downstream targets from being expressed at the same time. Since SP5 activates neural plate border genes [37] and AXUD1 promotes the expression of specification genes [13], this regulatory architecture prevents any overlap between these two GRN modules.
Temporal operation of the sub-circuit involves the transition between the two states shown in Fig 6L. The activation of AXUD1 and the specification program depends upon a reduction in the levels of SP5 and on the presence of TFAP2A, MSX1 and ZIC1. Curiously, neural plate border genes are also regulated by SP5 [37], which directly interacts with tissue-specific enhancers located at their loci (S6D Fig). This dual role of SP5, as both an activator (of neural plate border genes) and a repressor (of specification genes) is consistent with studies that have shown that it is able to exert both functions [36,38]. We postulate that once SP5 promotes the neural plate border program, the expression of these genes is stabilized by numerous positive feedback loops [4]. SP5 is subsequently downregulated, which releases repression of AXUD1 and allows specification to occur. Like other multipotency factors, SP5 expression decreases during avian neural crest development (S7D Fig), presumably due to a reduction in levels of Wnt activity [39] (in Xenopus, expression of SP5 can be detected in migratory cells, but there is still a reduction in expression of the gene [37]). Nonetheless, two features of this sub-circuit control the separation of the neural plate border and specification programs. First, SP5 needs to activate TFAP2A, MSX1 and ZIC1 before these genes can promote AXUD1 expression. Second, the direct repression of AXUD1 by SP5 prevents the former from being transcribed prematurely.
These results highlight the importance of repressive interactions in the neural crest GRN. A significant shortcoming of the current version of the network is that it is built primarily of positive regulatory links; only a small number of studies have focused on transcriptional repressors, like SNAI2 [40]. This is likely a consequence of technical limitations (identifying inhibitory interactions is technically challenging) but also because the importance of repression in GRNs has often been underestimated. Genomic analysis has demonstrated that cell fate specification programs rely heavily on inhibitory interactions, which act to silence alternative fates [41], or to promote stepwise activation of gene expression. Additional cis-regulatory analyses or combinations of genomic approaches [42] will be required to delineate the role of repression in the neural crest GRN. Our findings suggest the existence of a set of inhibitory interactions that control the timing of developmental transitions during neural crest development. Elaboration of network circuitry and identification of additional sub-circuits will be necessary to shed further light on the molecular logic of neural crest development.

Ethics statement
All experiments were performed in accordance with Memorandum of Understanding and Agreement 16274-1, which was approved by the Institutional Biosafety Committee of Cornell University.

Collection of chick embryos
Fertilized Leghorn White chicken eggs were obtained from the Department of Animal Science, University of Connecticut. Eggs were incubated at 37˚C until embryos reached the desired developmental stages. Embryos were collected and cultured according to the EC protocol [43], and staged using the Hamburger and Hamilton staging system [44].

In Situ Hybridization
For in situ hybridization, embryos were fixed in phosphate buffer saline (PBS) containing 4% paraformaldehyde (PFA) for 2 hours at RT. Following fixation, embryos were dissected, washed with PBST, dehydrated and stored in methanol at −20˚C. Whole-mount in situ hybridization was performed as previously described [45]. For double in situ hybridization, we used the Tyramide TSA system from Perking Elmer (TSA Plus Cyanine 5 and Fluorescein, NEL754001KT) as previously described [46].

3C-qPCR
Chromosome Conformation Capture (3C) libraries were generated based on previously published protocols [47][48][49] with modifications to account for lower amounts of starting material. Neural folds were dissected from HH9 embryos (n = 100-120 per replicate) in Ringers solution and dissociated in Accumax (Innovative Cell Technologies, #AM105). After crosslinking, cells were washed in cold PBS complemented with protease inhibitors, resuspended in Lysis Buffer and kept on ice for 10 min. Lysed cells were then incubated in 1.2x CutSmart (New England Biolabs) restriction buffer and incubated for 1h at 37˚C with 0.3% SDS. After addition of 2% Triton X-100 lysed cells were re-incubated for 1h at 37˚C. DNA digestion was performed with 400UI NCOI (New England Biolabs, #R3193) at 37˚C overnight. After enzyme inactivation with 1.6% SDS and overnight incubation at 65˚C, cells were transferred to 1.15x ligation buffer (New England Biolabs, #M0202) and incubated with 15% Triton X-100 for 1h at 37˚C. Ligation was performed for 4h at 4˚C with T4 DNA ligase (New England Biolabs, #M0202). Eluted DNA was treated with 10mg/ml Proteinase K overnight at 65˚C. After RNase treatment, samples were purified with phenol-chloroform, concentrated using Amicon Ultra-0.5 ml 30k columns (Millipore, #UFC5030BK) and quantified in Qubit. A control library was assembled by mixing equimolar amounts of PCR products spanning the AXUD1 locus with minimal overlap. Amplification efficiency (slope and intercept) of primers pairs (fragment primers plus constant primer) was verified via qPCR of 10x dilutions of digested/ligated control library [50]. Interaction frequencies were determined by the normalization of ligation products/loading control (GAPDH) as described in [50] (S3 and S4 Tables). The 3CqPCR primers are described in S5 Table.

Quantitative reverse transcription PCR (RT-PCR)
eRNA expression levels were determined through RT-PCR in microdissected neural crest cells. The tissue was lysed in lysis buffer from RNAqueous-Micro Total RNA Isolation Kit (Thermo-Fisher, #AM1931) and the total RNA was isolated following manufacturer's protocol. cDNA was then synthesized using SuperScript III Reverse Transcriptase and Random Primers (Ther-moFisher, #18080051) according to the kit's protocol. RT-PCR was performed using Power Sybr Green PCR master mix (Thermo Fisher, 4368577) in an ABI viia7 RT-PCR machine. eRNA levels at the AXUD1 locus are presented normalized to reference gene HPRT1.
To quantify changes in gene expression caused by SP5 loss-of-function or gain-of-function, we microdissected single neural folds from control and targeted sides of the embryo, which were subsequently lysed in lysis buffer from Power SYBR Green Cells-to-CT Kit. RNA extraction and cDNA preparation were performed according to the kit's protocol). RT-PCR was performed using Power Sybr Green PCR master mix (Thermo Fisher, 4368577) in an ABI viia7 RT-PCR machine. Ct values of all genes were normalized to reference gene HPRT1 and expressed as a fold change compared to the control sample. The qPCR primer sequences are listed in S5 Table. Transient transgenesis Enhancer plasmids, morpholinos and expression vectors were transfected in chick embryos at HH4 by ex ovo electroporation, as previously described [13]. Constructs were injected between the epiblast and vitelline membrane of embryos at a concentration of 1-2ug/ul and electroporated with platinum electrodes (five 50ms pulses of 5.1V, with an interval of 100ms between pulses). In all gene knockdown and overexpression experiments, the embryos were injected bilaterally with the control reagent on the left side and the targeted reagent on the right side. Following electroporation, embryos were cultured in albumin at 37˚C until they reached appropriate developmental stages. Embryo survival was >90% and all embryos were screened to ensure uniform electroporation and proper embryo morphology prior to further downstream analysis.

Enhancer-reporter assays
Putative enhancers (S1 Table) defined by DNA accessibility and TFAP2A binding profiles in neural crest cells were amplified from HH10 chicken genomic DNA and cloned in the pTK-eGFP vector [19]. Single cell measurements with flow cytometry were performed to quantify the intensity of enhancer variants. Heads from three HH9 embryos co-electroporated with the pTK-eGFP construct and Axud1E1-2500bp:mCherry, were dissected and dissociated in Accumax (Innovative Cell Technologies, #AM105). 500 cells per sample were analyzed and plotted in categorical scatter plots using Seaborn Python library. ZIC1, MSX1 and SP5 binding sites in Axud1E1-300 were identified by scanning this sequence using the JASPAR database of transcription factor binding profiles, with a minimal relative profile score threshold of 80% [51]. Mutant constructs (S2 Table) were cloned into pTK-eGFP, and compared to the wild-type enhancer in double-sided electroporations.

Knockdown assays
All loss-of-function assays were performed with double-sided injections of morpholinos and their respective controls to allow single-embryo internal controls. Wnt knockdown was performed by the combined inhibition of WNT 1 and 4 by morpholinos at 1.25uM each, while CTNNB1, SP5, MSX1, ZIC1 and TFAP2A morpholinos were used at 1.5uM (S5 Table). Embryos were incubated until the desired stage and imaged to evaluate enhancer reporter activity.

SP5 overexpression
For overexpression assays the SP5 expression constructs were assembled by insertion of the full-length cDNA sequence of avian SP5 in a pCI-H2B-RFP backbone [52]. The coding sequences were PCR amplified from an HH8 cDNA library. The SP5 expression construct was electroporated as described above, and paired single neural folds (left side = control, right side = knockdown) were dissected and processed individually for RT-PCR.

CRISPR-Cas9 mediated enhancer loss-of-function
We employed a CRISPR-Cas9 system optimized for chick embryos to disrupt the activity of Axud1E1-500 [23]. Two gRNAs were designed using online resources (crispor.tefor.net) and cloned downstream of the U6 promoter in the cU6.3 vector (S5 Table). To assess the effect of endogenous Axud1E1-500 knockdown, gastrula-stage embryos were electroporated with a pCAGG-nls-Cas9-nls-GFP vector and the Axud1E1-500 gRNAs. A control gRNA was used on the left side of the embryo. Embryos were re-incubated at 37˚C. At stage HH9, embryos were screened for robust GFP expression in both sides, and half heads were dissected for control and targeted sides of the embryo. Pools of three half heads were then dissociated in Accumax (Innovative Cell Technologies, #AM105) for 30 min. After dissociation, cells were resuspended in HANKS solution supplemented with 0.5% BSA. Control and target GFP+ cell suspensions (150-300 cells) were sorted into 50 μl of lysis buffer from the Power SYBR Green Cells-to-CT Kit (ThermoFisher, 4402955) using a BD AriaFusion cell sorter. Samples were processed following the manufacturer's protocol. Axud1E1-500 and Axud1E1-SE expression levels were determined using RT-PCR as described above. cU6.3 and pCAGG-nls-Cas9-nls-GFP vectors were a gift from Dr. Marianne Bronner.

ATAC-qPCR
To examine the control of Axud1E1 accessibility, we performed ATAC-qPCR in embryos transfected with a TFAP2A morpholino. HH4 embryos were bilaterally injected with control and targeted morpholinos, as described previously. Embryos were incubated at 37˚C until HH9, when dorsal neural folds were surgically dissected. Paired single neural folds were individually processed for DNA tagmentation following the ATAC protocol described in [17]. For each genomic location, specific enrichment was quantified by RT-PCR (S5 Table). Raw CT values were first normalized to a control region defined by the absence of transcription factor binding, no enrichment of active histone marks and low DNA accessibility in our CUT&RUN and ATAC datasets. Changes in chromatin accessibility in morpholino treated samples were subsequently compared to control samples.

CUT&RUN
Neural folds were dissected from HH7-8 embryos (n = 20 per CUT&RUN experiment). Cells were dissociated in Accumax for 20min at RT under mild agitation. CUT&RUN experiments were carried out as previously described [17]. Briefly, cells were immobilized on BioMag Plus Concanavalin A magnetic beads (Bangs Laboratories, BP531) and incubated with rabbit anti-Sp5 (Abcam, # ab36593) antibody (1:50) overnight at 4˚C. After washing away unbound antibody, protein A-MNase was added to a final concentration of 700ng/mL and incubated for 1h at 4˚C. Cells were cooled to 0˚C and CaCl2 was added to a final concentration of 2mM to activate the MNase enzyme. MNase digestion was performed for 45min and terminated by the addition of 2XSTOP buffer. The protein-DNA complexes were released by centrifugation and digested with proteinase K for 10 min at 70˚C. DNA fragments were isolated via phenol-chloroform extraction and ethanol precipitation. Protein A-MNase was kindly provided by Dr. Steven Henikoff [27].

CUT&RUN library preparation
CUT&RUN libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, #E7645) following the manufacturer's protocol. Fragment analysis was performed with ABI 3730xl DNA Analyzer to perform quality control for the libraries. Equimolar concentrations of the libraries were pooled using the KAPA Library Quantification Kit-ROX Low (Roche, #07960336001) and sequenced with paired-end 37bp reads on an Illumina Next-Seq500 instrument.

CUT&RUN data analysis
Paired-end sequencing reads from the CUT&RUN libraries were trimmed using Cutadapt [53]. Reads were filtered for those with a minimum length of 25bp or longer and aligned to the reference chicken Galgal5 assembly using Bowtie2 [54]. Picard MarkDuplicates tool was used to mark duplicate reads and BAM files were filtered with SAMtools to discard unmapped reads (those that were not the primary alignment, reads failing platform/vendor quality checks, and PCR/optical duplicates (-f 2 -F 1804). Peak calling was performed using MACS version 2.1 with a pValue cutoff of 0.01, skipping the shifting model and extending read sizes to 200bp (nomodel-extsize 200). Representative heatmaps showing the SP5, H3K27Ac and IgG enrichment at SP5 bound peaks were generated using the deepTools2 package [55]. The GO-category analysis was performed with the clusterProfiler R package [56] to assay for over-represented Biological Processes,with a pValue cutoff of 0.05. Motif enrichment analysis was performed using the HOMER findMotifsGenome package [57].
Supporting information S1  Chromatin immunoprecipitation for CTNNB1 shows temporal changes in the regulation of SP5 and AXUD1 by Wnt signaling. Association of CTNNB1 with an SP5 enhancer (Sp5E) decreases during neural crest specification. Conversely, binding of the Wnt effector to Axud1E1 significantly increases from HH5 to HH9. Error bars represent ± SEM. The statistical significance was determined via an unpaired t-test. (B) Fragment size distribution of the two replicates of SP5 CUT&RUN read pairs. Motif enrichment analysis via HOMER for regions occupied by SP5 shows enrichment for GC boxes, similar to other Sp/Klf Zn2+-finger transcription factors. pValue indicates significance of motif occurrence as reported by HOMER.