Identification of in vivo Hox13-binding sites reveals an essential locus controlling zebrafish brachyury expression

All vertebrate embryos form their anterior-posterior axis by adding cells progressively to the posterior end of the presomitic mesoderm and spinal cord from a neuromesodermal progenitor population (NMp) located in the tailbud at the most posterior end of the body (Henrique et al., 2015; Kimelman, 2016b; Mallo, 2020; Steventon and Martinez Arias, 2017; Wymeersch et al., 2021). A positive autoregulatory loop consisting of Wnt signals and the T-box transcription factor Brachyury (Tbxt) plays an essential role in the NMp population by sustaining the mesodermal potential of NMps, such that loss of Brachyury causes a severe reduction in somitic mesoderm and an increase in neural tissue (Gouti et al., 2014; Martin and Kimelman, 2012; Yamaguchi et al., 1999). In contrast, too much Wnt and Brachyury are also detrimental because the mesoderm becomes overpopulated at the expense of neural tissue (Bouldin et al., 2015; Garriock et al., 2015; Gouti et al., 2014; Jurberg et al., 2014; Martin and Kimelman, 2012; Tsakiridis et al., 2014; Turner et al., 2014). Hence, the Wnt-Brachyury autoregulatory loop needs to be controlled carefully in order for the embryo to produce the correct amount of mesoderm and neural tissue as the body extends; however, little is known about the precise mechanisms controlling the levels of these factors.

We recently showed that the Hox13 genes, specifically hoxa13b and hoxd13a, enhance the ability of the zebrafish brachyury ortholog tbxta (also known as no tail) to promote formation of the mesoderm during the somitogenesis stages by helping to maintain expression of tbxta, wnt3a and wnt8 (Ye and Kimelman, 2020). However, the direct targets of these Hox13 factors were unknown. Although the importance of Hox genes in embryonic development has been known for decades, identifying in vivo targets has been a major challenge, particularly in vertebrates, and only a limited number of verified vertebrate Hox targets are known (reviewed by Mann et al., 2009). Although large-scale in vitro screens have been useful for identifying DNA-binding motifs for Hox proteins, both in Drosophila and vertebrates (Berger et al., 2008; Jolma et al., 2013; Noyes et al., 2008), the sequences are relatively nonspecific, making in silico predictions of Hox-binding sites challenging, and essentially impossible in cases where transcription is regulated by distal enhancer elements. However, an important finding from one of these studies, which used protein produced in cell culture, is that the posterior Hox factors, Hox9-13, bind not only the canonical Hox motif shared by all Hox factors, but also a slightly different alternative binding motif (Jolma et al., 2013).

One major problem for identifying in vivo Hox-binding sites in vertebrate embryos is the paucity of high-quality antibodies specific to Hox proteins, which is an even bigger problem for zebrafish than for mammals because of less commercial interest. In addition, approaches, including chromatin immunoprecipitation sequencing (ChIP-seq), are not feasible for tissues such as the tailbud, which contains only a few hundred cells. We solved the first of these problems by creating two transgenic zebrafish lines that allow us to produce controlled levels of Hoxa13b with a small epitope tag, circumventing the need for a Hox13-specific antibody. We have also adapted a relatively new technology, called CUT&RUN, for use on zebrafish embryos, which permits quantitative identification of in vivo binding sites for transcription factors using very limited numbers of cells (Meers et al., 2019; Skene and Henikoff, 2017).

Using these approaches, we identified at least 576 candidate Hoxa13b target genes in the zebrafish tailbud during the middle of the somite-forming stages. Using epigenetic markers in parallel, we further characterized the targets associated with transcriptional activation to identify genes that are potentially activated by Hoxa13b. Importantly, we found two transcriptionally active Hoxa13b-binding sites far upstream of the tbxta transcription start site that have been conserved over hundreds of millions of years of evolution. These enhancer elements also each contain a conserved perfect consensus-binding site for Tcf7/Lef1, which are the Wnt-regulated transcription factors in cells, revealing how Hox factors and Wnt signaling could intersect to promote tbxta expression. We reveal that the proximal of the two elements, which in zebrafish embryos shows higher levels of Hoxa13b binding in the CUT&RUN analysis, is sufficient to drive posterior expression of tbxta in the tailbud when joined to a proximal tbxta promoter, resolving the previously puzzling inability in zebrafish to find a tbxta promoter that sustains expression in NMps during axis extension. Importantly, we show that deletion of the proximal enhancer element by using CRISPR produces shorter embryos because of defects in the most posterior somites, as well as reduced expression of tbxta in NMps. Finally, we reveal an unexpected role for the transcription factor Rbpj, which binds near the Hox13-binding sites in tbxta, in regulating tbxta expression and posterior body formation. In addition to contributing a rich source of data for understanding gene expression within the tailbud, our results provide a potential molecular basis for understanding how Hox13 proteins act to promote axis formation, and reveal a crucial Hox-binding enhancer element regulating the Brachyury/Wnt autoregulatory loop.

In order to solve the lack of high-quality antibodies to Hoxa13b, we produced two different transgenic lines that each add a FLAG epitope tag to Hoxa13b. One of the lines produced a heat shock-inducible Hoxa13b-FLAG, whereas the other placed the Hoxa13b-FLAG tag at the endogenous hoxa13b locus. We used both approaches because the expression level of endogenous hoxa13b in the tailbud is low, as is true for all Hox13 genes (Ye and Kimelman, 2020), and we were concerned that this low expression might cause us to miss some of the authentic Hoxa13b target genes because of a weak signal. By using both approaches, we hoped to maximize our chances of identifying the complete set of Hoxa13b-binding sites. Inducible epitope-tagged Hox factors have recently been used in differentiated mouse embryonic stem cells to identify Hox targets in that system, confirming the utility of this approach (Bulajic et al., 2020).

To produce the inducible transgenic line, 27 amino acids containing three FLAG tags were added to the C terminus of Hoxa13b (Fig. 1A). In order to monitor expression, GFP was placed after the 2A peptide, which allows GFP to be co-expressed with Hoxa13b as a separate protein. This coding region was placed under control of a zebrafish heat-shock promoter to make the line HS:hoxa13b-FLAG-GFP. We previously showed that a similar transgenic line that produces a nonepitope-tagged version of Hoxa13 produces a severe axis formation defect when heat shocked at 40°C (Ye and Kimelman, 2020). This new line produces the identical phenotype when heat shocked under the same conditions, demonstrating that the epitope tag does not impair Hoxa13b function (Fig. S1).

The second transgenic line was produced by taking the genomic hoxa13b sequence, which includes one intron, and placing the FLAG-2A-GFP sequences at the C terminus of the coding region. This was flanked by a unique DNA sequence called Mbait, which is not present in the zebrafish genome, allowing the plasmid to be cut in vivo by a specific Mbait gRNA without creating off-target effects (Kimura et al., 2014). This construct was injected into zebrafish embryos along with Cas9, the Mbait gRNA and gRNAs targeting the 5′ and 3′ untranslated regions (UTRs) of the endogenous hoxa13b locus in order to insert the modified hoxa13b at the endogenous hoxa13b locus. We identified one line that showed the expected expression of the FLAG-tagged protein at the C terminus of the embryo in NMps and mesoderm (Fig. 1C-G, Movie 1). This knock-in (KI) transgenic line is referred to as KI:hoxa13b-FLAG-GFP.

To identify in vivo Hoxa13b-binding sites, we adapted the CUT&RUN procedure (Meers et al., 2019; Skene and Henikoff, 2017) for use in zebrafish. In order to produce Hoxa13b-FLAG in embryos without perturbing development from HS:hoxa13b-FLAG-GFP, we first established that heat shocking at 36°C produced just enough GFP to detect, but did not cause any phenotypic alterations to embryonic development. We then manually isolated tailbuds from embryos of both the inducible and KI transgenic lines at the 15-18 somite stage, which is in the middle of the somite-forming stages, and then dissociated and permeabilized the cells. Based on the expression pattern in the KI line (Fig. 1F, Movie 1), we estimate that ∼50% of the cells in our dissected tailbuds produced Hoxa13b-FLAG protein. Anti-FLAG antibody was added, followed by Protein AG-Micrococcal nuclease (pAG/MNase), which cuts just adjacent to the FLAG-tagged transcription factor (Fig. 1B). Once the genomic DNA is cut, the complex is released into solution and the associated DNA is subjected to DNA sequence analysis. In parallel, we isolated tailbuds from wild-type (WT) embryos and performed CUT&RUN using antibodies to H3K4me1 and H3K27ac in order to identify the regions of the genome in the tailbud associated with transcriptional activation (Calo and Wysocka, 2013). Finally, as a control, we performed CUT&RUN on WT embryos with the anti-FLAG antibody, which was not expected to bind any endogenous zebrafish transcription factors.

Given that CUT&RUN has not been performed previously in zebrafish, we first wanted to be sure that it could be done effectively. We evaluated the success of our CUT&RUN experiment by checking the size distribution of all aligned fragments for the histone modification samples. We observed a strong enrichment of mononucleosome-sized fragments, with a small number of dinucleosome-sized fragments, which is a typical pattern for histone modification CUT&RUN data from a successful experiment (Fig. 2A,B; Meers et al., 2019; Skene and Henikoff, 2017). In contrast, the control antibody samples showed a random distribution of sizes. Correlation analysis of the samples from independent experiments showed good reproducibility among biological replicates and that the H3K27ac and H3K4me1 modifications were often found together, as expected (Fig. 2C). Finally, we used heat maps to examine signal enrichment at specific locations, as used previously (Henikoff et al., 2020; Skene and Henikoff, 2017). Whereas DNA fragments produced from the control DNA showed no enrichment at the called-peak sites, DNA fragments produced using the histone modification antibodies showed central enrichment at these sites (Fig. 2D,E).

Given that the most comparable ChIP-seq data using the histone modifications were produced with 24 h post-fertilization (hpf) whole embryos (Bogdanovic et al., 2012; de la Calle Mustienes et al., 2015) whereas we used posterior explants at 16-17 hpf, we were not able to compare our CUT&RUN results directly to stage-matched ChIP-seq data through bulk analysis. However, comparing a variety of posteriorly expressed genes that remain on at both time points, such as sox2 and cdx4, showed similar results for the histone modifications in both methods (Fig. 2I). The results of the histone modification analysis are provided in Tables S1 and S2, and will be useful for identifying transcriptional regulatory regions for genes involved in NMps, presomitic mesoderm formation and spinal cord neurogenesis. In summary, our results show that CUT&RUN can be used effectively in zebrafish. In addition, whereas the previous ChIP-seq data required 1000 whole embryos per biological replicate (Bogdanovic et al., 2012; de la Calle Mustienes et al., 2015), our results were obtained with 30 tailbuds per replicate, opening up this type of genomic analysis to a variety of studies in zebrafish that were not previously possible, such as analyzing chromatin signatures in small regions of the embryo or tissue-specific fragments.

We next proceeded to use CUT&RUN on our transgenic lines expressing Hoxa13b-FLAG. Results from both lines showed clear central enrichment of signal at the called-peak regions (Fig. 2F,G). Comparing the results using HS:hoxa13b-FLAG-GFP and KI: hoxa13b-FLAG-GFP, we found that the heat-shock line had a better signal-to-noise ratio; thus, we used this line for initial identification of Hoxa13b-binding sites, and then used the data from the KI line as a second level of analysis for peaks called from the inducible line to find sites identified in both transgenic lines. The reduced signal from the KI line resulted not only from the lower level of expression of Hoxa13b (based on relative GFP expression), but also because, at lower levels of expression, the Hoxa13b-FLAG was competing for binding against other Hox13 factors, and potentially all Hoxa9-13 proteins (see Discussion), making it more difficult to detect binding. Given that the KI data are likely to under-represent endogenous Hoxa13b-binding sites, as shown below for one of the enhancer elements in tbxta, we included peaks in Table S3 that were found with both the inducible and KI lines as well as peaks found only with the inducible line. Table S3 also includes the associated histone marks discussed above for each of the peaks.

Using HS:hoxa13b-FLAG-GFP embryos, we identified a total of 1871 Hoxa13b-occupied peaks, 689 of which were overlapping peaks identified using the KI:hoxa13b-FLAG-GFP line. Only 12.6% of the 1871 binding sites were located within 3 kb of the transcription start site of the nearest gene, whereas more than 60% of the sites were in the distal intergenic regions (Fig. 3A, Fig. S2A). The same result was observed in a recent ChIP-seq analysis of the binding sites of multiple Hox factors in differentiated mouse embryonic stem cells (Bulajic et al., 2020). Based on the TxDb.Drerio.UCSC.danRer11.refGene annotation database, we mapped the 1871 binding sites identified from the HS:hoxa13b-FLAG-GFP analysis to the nearest neighboring gene, which identified 1380 genes. The 689 sites overlapping in the KI data mapped to 576 neighboring genes (Table S3).

We analyzed the DNA peaks for common motifs using the programs MEME and DREME, and found that a Hox motif was the most common motif discovered (Fig. 3B, Fig. S2B). The motif was centrally located in the peak region (Fig. 3C), with footprinting analysis showing that it forms a 19-bp footprint (Fig. 3D) because the bound transcription factor protects the DNA from MNase cutting (Liu et al., 2018; Skene and Henikoff, 2017). Intriguingly, at a specific position in the motif (Fig. 3B), approximately half the fragments had a T, which is the base bound by all Hox proteins (Hox1-13), whereas, in the other half, a C was detected, which is bound specifically by posterior Hox proteins (Hox9-13; Jolma et al., 2013). Thus, Hoxa13b binds both in vitro-predicted motifs (TTTAT and TTTAC) in our system. We did not note a clear preferential association of T versus C with other enhancer marks. In addition, apart from the Hox motif, we did not consistently identify any obvious secondary motifs, suggesting that the posterior Hox proteins, unlike the more anterior Hox proteins that bind TALE co-factors (Mann et al., 2009), rely instead on a diverse group of co-factors in order to achieve specific binding.

We next examined the list of genes bound by Hoxa13b for targets that could explain how Hox13 proteins might support the Brachyury/Wnt loop. We were particularly intrigued to find a Hoxa13b-binding site neighboring the tbxta gene, located 24.6 kb upstream of the start of transcription (Fig. 4A), which we called ‘Hox element 1’. This site was associated with the active enhancer marks H3K4me1 and H3K27ac (Fig. 4A), suggesting that this region is an active transcriptional enhancer in the tailbud. This was of particular interest because a major puzzle with tbxta in zebrafish is that a 2.1-kb fragment upstream of the transcription start site drives strong expression throughout NMps during the gastrula stages and later in the notochord, but does not produce expression in NMps during the somitogenesis stages when the axis extends (Harvey et al., 2010), suggesting that a hitherto unidentified distal enhancer is required to promote the second phase of neuromesodermal tbxta expression. Indeed, in our previous work showing the important effects of Hox13 factors on the Brachyury/Wnt loop, we speculated that the Hox proteins might work through a yet unidentified distal enhancer for tbxta (Ye and Kimelman, 2020).

When we used BLAST to compare this tbxta distal enhancer region, including the Hoxa13b-binding site, to other fish, we found that an ∼200 bp region was conserved in all the ray-finned fish we could analyze, including reedfish from the genus Polypterus, which diverged from all other actinopterygians (ray-finned fish), including zebrafish, ∼400 million years ago (Takeuchi et al., 2009), indicating a strong evolutionary pressure to retain this element (Fig. 4B). The Hox-binding motif is conserved completely among all these fish and, intriguingly, all of them contain the C that is specific for binding by posterior Hox proteins (Fig. 4B, arrowhead).

We were also interested to observe that Hox element 1 contains a perfect match to the consensus-binding site for Tcf7 and Lef1 transcription factors, which are the transcription factors activated by the Wnt pathway (Cadigan and Waterman, 2012), because Wnt expression is required for brachyury/tbxta expression throughout early development (Arnold et al., 2000; Galceran et al., 1999; Martin and Kimelman, 2008; Vonica and Gumbiner, 2002; Yamaguchi et al., 1999). This Tcf7/Lef1-binding sequence is exactly the same one used in the Super TOPflash Wnt reporter (Veeman et al., 2003). Therefore, the combination of a Hox13-binding site and a Tcf7/Lef1-binding sequence suggests that Hox element 1 integrates Hox13 binding and activation by the Wnt pathway.

In examining the CUT&RUN data from the HS:hoxa13b-FLAG-GFP fish, we saw a weaker peak of Hoxa13b binding 28.9 kb upstream of the start of transcription start site in the tbxta gene that was not called as a peak by our algorithm nor was a peak present in data from the KI:hoxa13b-FLAG-GFP fish (Fig. 4A). When we examined this region using a BLAST search, we found that it is also highly conserved, with conservation as far as the coelacanth, a lobe-finned fish (Fig. S3). Both the posterior-specific Hox motif and a consensus Tcf7/Lef1 site were found in this element, suggesting that it serves as a second enhancer element for tbxta. This site is also associated with the activating marks H3K4me1 and H3K27ac (Fig. 4A), demonstrating that it is also an active enhancer in the tailbud. Therefore, we called the region at 28.9 kb ‘Hox element 2’.

To detect a role for Hox elements, we inserted them, along with nonconserved flanking sequences in order to allow for DNA looping, into a Tol2 vector that contained GFP and the 2.1 kb proximal promoter (Fig. 5A). We noticed a site that had strong H3K4me1 and H3K27ac occupancy at −5.4 kb, which was also conserved among different fish species, although no recognizable transcription factor motifs were found in this element (Fig. S4). Therefore, we included it in our initial analysis in case it was important for enhancing tbxta transcription. As shown previously (Harvey et al., 2010), transgenic zebrafish lines containing the 2.1 kb proximal promoter produce strong GFP expression during gastrula stages and at later stages in the notochord, but little GFP expression in the tailbud and posterior somites during somitogenesis stages (Fig. 5B, Fig. S5). In contrast, a transgenic line that included the proximal promoter, the two Hox elements and the −5.4 kb element produced robust expression in the somitogenesis-stage tailbud and posterior somites (Fig. 5C, Fig. S5). We then made a new transgenic line that included only Hox element 1 and its flanking sequences, and found that this also produced strong tailbud expression (Fig. 5D). These results demonstrate that Hox element 1 is sufficient to confer somitogenesis-stage tailbud expression upon the proximal tbxta promoter.

Although, in most fish, these Hox elements are located far from the start of tbxta transcription, we noticed that they were located much closer together in channel catfish (Fig. 5A and Table 1). Therefore, we isolated 5.7 kb of the catfish tbxta upstream region and derived transgenic zebrafish lines containing this construct. The catfish promoter was able to drive strong tailbud expression in zebrafish embryos (Fig. 5E), demonstrating functional conservation of the tbxta upstream region despite it being much smaller in catfish. Interestingly, although the three elements shown in Fig. 5 are conserved between zebrafish and catfish, we could not detect any conservation in the 2.1 kb proximal promoter region between these species, as was also seen previously in a comparison between zebrafish and medaka for the same promoter region (Harvey et al., 2010).

Although our data showed that Hox element 1 is sufficient to regulate tbxta expression in the posterior, we wanted to know whether it is necessary. Thus, we used CRISPR with guide RNA (gRNAs) that flanked Hox element 1 (Fig. 6A), isolated two transgenic lines that had this element deleted and bred these through two generations. Fish Ptbxta^Δ312 had a 312 bp deletion and Ptbxta^Δ314 had a 323 bp deletion and a 9 bp insertion; both produced the same phenotype. In crosses of fish heterozygous for the deletion, a quarter of the fish had a shorter body (Fig. 6B,G; 25%, n=92). When a subset of the embryos were genotyped, all of the short bodies were homozygous for the deletion (100%, n=6), whereas the normal-length bodies were heterozygous or WT (100%, n=6). Examining these shorter-bodied fish revealed no change in the number of somites and, instead, we observed aberrant formation of the most posterior somites compared with WT (Fig. 6H, Fig. S6A). We also examined tbxta expression because we expected this to be changed. At 24 hpf, we identified a portion of embryos with greatly reduced tbxta (Fig. 6C,D). Genotyping of these embryos showed that 100% were homozygous for the deletion (n=13). Given that reduced tbxta expression would be expected to lead to reduced mesoderm formation, we examined the early mesodermal marker tbx16 4 h later (28 hpf) and observed embryos with strongly decreased tbx16 (Fig. 6E,F). Genotyping of these embryos showed that 100% were homozygous for the deletion (n=14). Thus, deletion of Hox element 1 causes strongly reduced tbxta and tbx16 expression and a shorter body.

Given that NMps are bipotential, there is a corresponding increase in neural tissue when the amount of mesoderm decreases (Garriock et al., 2015; Gouti et al., 2014; Martin and Kimelman, 2012; Tsakiridis et al., 2014; Turner et al., 2014). To see whether this occurs in fish lacking Hox element 1, we examined the volume of the neural tube using the neural marker sox2 with fluorescent in situ hybridization (FISH). Whereas WT and heterozygous embryos had the same neural tube volume, embryos homozygous for the Hox element 1 had a markedly increased volume (Fig. 6I, Fig. S6B). Thus, when Hox element 1 was deleted from tbxta, the NMps switched from producing mesoderm to neural tissue in the most posterior end of the embryo.

In addition to Hox- and Tcf7/Lef1-binding sites in Hox element 1, we also found a perfect match to the consensus sequence for Rbpj (Fig. 4B; Johnson and Macdonald, 2011). Rbpj, which in Drosophila is known as Suppressor of Hairless, is best known for its transcription-repressing roles in the Notch signaling pathway, but it also has important Notch-independent roles, which include both activation and repression (reviewed by Johnson and Macdonald, 2011). Zebrafish have two duplicated rbpj ohnologs, rbpja and rbpjb, both of which have maternal and ubiquitous expression through the first day of development (Echeverri and Oates, 2007; Sieger et al., 2003). Consistent with their role in the Notch pathway, knockdown of both ohnologs using morpholino oligonucleotides (MOs) causes defects in cyclic gene expression in presomitic mesoderm and somite boundary defects (Echeverri and Oates, 2007; Sieger et al., 2003). Interestingly, MO-knockdown embryos have posterior defects, particularly in tail formation, although this was not examined in detail. In mouse, null mutants of Rbpj die by 10.5 days post-fertilization and form fewer than six somites (Oka et al., 1995), which is reminiscent of the somite defects in brachyury (T) null mutants (Beddington et al., 1992). Unfortunately, the expression of brachyury was not examined in these embryos, but there was a strong reduction in the levels of the somite marker Meox1. These results show that Rbpj has roles that extend beyond regulation of the Notch pathway in vertebrate embryos, and involve both formation of the mesoderm and cell growth and survival.

In order to determine whether Rbpj binds Hox element 1, we made a new transgenic line, HS:rbpja-FLAG-GFP, which allows conditional expression of a FLAG-tagged version of Rbpja with co-expressed GFP. Embryos expressing Rbpja-FLAG after heat shock did not show any defects using the standard temperature of 40°C. Embryos co-expressing GFP were subjected to CUT&RUN as described above, and central enrichment of signal at the called-peak regions was observed (Fig. 2H). In agreement with the proposed role of Rbpj in Notch signaling, we observed peaks near the presomitic cycling genes her1, her4, her7, hes6, hey1 and hey2, and also at the Notch targets dlc and dld (Table S4). In addition, we observed peaks near cyp26a1, which was shown to be strongly reduced in the tailbud when rbpj MOs were used (Echeverri and Oates, 2007). Importantly, we identified binding at the consensus Rbpj site in Hox element 1 (Fig. S7). We also observed Rbpja binding to Hox element 2 (Fig. S7) and, although we could detect an Rbpj site close to Hox element 2 in all fish we examined, the position of the site was not conserved among species. In zebrafish, for example, the Rbpj site was 40 bp beyond the Tcf7/Lef1 site shown in Fig. S3.

Given the well-known issues with the use of morpholinos, we took an alternative dominant-negative approach to examine the role of Rbpj in zebrafish by creating a transgenic line in which the Engrailed repressor domain was expressed at the N terminus of Rbpja (HS:EnR-rbpja-FLAG-GFP). This construct was effective as EnR-Rbpja-expressing embryos heat shocked at 38°C showed a strong posterior defect, whereas HS:rbpja-FLAG-GFP embryos were unaffected (Fig. 7A). Expression of tbxta was examined 5 h after heat shock and was markedly reduced in HS:EnR-rbpja-FLAG-GFP embryos in the tailbud, but not in WT or HS:rbpja-FLAG-GFP embryos (Fig. 7B-D), whereas cdx1a and snai1a, which are genes that do not have Rbpja-binding sites (Table S4), were unaffected (Fig. 7E-J). In addition, notochord expression of tbxta was not inhibited by EnR-Rbpja (Fig. 7B-D). We also saw no increase in apoptosis in HS:EnR-rbpja-FLAG-GFP embryos 5 h post-heat shock, demonstrating that the reduction of tbxta was due to a change in gene expression and not a loss of cells. Thus, these results show that Rbpja also contributes to the activation of tbxta expression.

Whereas our previous work provided genetic evidence for a role of Hox13 proteins in sustaining NMps, the direct targets of the Hox13 proteins were unknown. Using two new transgenic zebrafish lines and CUT&RUN for identifying in vivo binding targets, we identified at least 576 candidate Hox13 direct-target genes in the tailbud, which potentially increases to almost 1400 if all of the candidate Hoxa13b-binding sites from the inducible transgenic line are included. This is the first use of CUT&RUN in zebrafish, and our results demonstrate that this method will be valuable because a high-quality antibody for a particular transcription factor is not required as long as a small epitope tag can be added and small amounts of tissue can be used, eliminating background signals from other regions of the embryo. Although direct targets of specific transcription factors have been identified in zebrafish previously using ChIP-seq (Jahangiri et al., 2012; Morley et al., 2009), this used large numbers of whole embryos (5000 embryos per biological replicate). Even analysis of chromatin signatures required 1000 whole embryos per replicate (Bogdanovic et al., 2012; de la Calle Mustienes et al., 2015). In contrast, we used 30 tailbud explants per biological replicate, demonstrating that transcription factor-binding data can be ascertained from very small embryonic tissue samples, greatly expanding the ability to perform specific genomic analysis on defined regions of the embryo. This same approach should also work on other embryos, particularly where transgenic approaches are feasible.

The data provided here should be of widespread use for researchers studying many aspects of posterior body formation, including transcriptional analysis of somite boundary formation, because we included analysis in WT embryos of two transcription-activating epigenetic marks. Moreover, we analyzed only one of the more than 500 Hoxa13b-bound genes shown in Table S3, leaving many other interesting targets for future studies. Here, we focused on a transcription-activating enhancer of somitogenesis-stage tbxta expression that we call Hox element 1, which was of interest because, in our previous work (Ye and Kimelman, 2020), we speculated that a distal tbxta enhancer could be a Hox13 target based on the failure of a proximal 2.1 kb promoter fragment to drive expression in NMps after the end of gastrulation. This enhancer is the missing regulatory element that had been revealed in a previous study of the tbxta promoter in zebrafish (Harvey et al., 2010), and is of high practical value because it was not possible previously to drive expression of constructs in NMps in zebrafish beyond the gastrula stage, except for the notochord (Row et al., 2016).

The importance of Hox element 1, as well as Hox element 2, is supported by their conservation in all the ray-finned fish examined, spanning 400 million years of evolution. Interestingly, the distance from these elements to the start of transcription is variable (Table 1), demonstrating that much of the upstream region is not essential for tbxta expression. Indeed, we found that most of the tbxta upstream sequences are not conserved among fish species and neither is the proximal tbxta promoter, as previously noted (Harvey et al., 2010). Although we were not able to find sequences obviously matching these Hox-binding elements in other vertebrates, relatively simple DNA rearrangements or changes in the sequences other than the Hox and Tcf7/Lef1 motifs will preclude finding homology using similarity searches. Especially because the Hox-binding motif is so degenerate, our results show the importance of beginning with in vivo identification of Hox-binding sites.

Although our previous loss-of-function studies (Ye and Kimelman, 2020), and overexpression studies in mouse and chick (Aires et al., 2019; Denans et al., 2015; Young et al., 2009), focused on Hox13 factors, our observation that tbxta Hox-binding elements, as well as binding sites in many other genes, contain a motif specific to the posterior Hox factors in vitro (Jolma et al., 2013) suggests that these factors act in combination to regulate posterior body formation. A recent paper examining binding of HOXC9, HOXC10 and HOXC13 in an ex vivo system (motor neurons produced from embryonic stem cells) using induced epitope-tagged HOX proteins and ChIP-seq, revealed that HOXC13 binds many sites not bound by HOXC9 or HOXC10, although ∼10% of the sites bound all three factors (Bulajic et al., 2020). Therefore, it will be interesting to use CUT&RUN to determine which of the Hox13 in vivo targets in the tailbud are also bound by Hox9-12 to determine which genes are likely to be solely Hox13-regulated genes and which are controlled more generally by all of the posterior Hox genes. The posterior Hox proteins have always been somewhat confusing because, unlike the more anterior Hox proteins that bind TALE co-factors such as Pbx and Meis to increase local DNA-binding specificity, specific co-factors for the posterior Hox proteins have not been identified (Mann et al., 2009). We were unable to find a specific motif commonly associated with the Hoxa13b-binding site in our large dataset of Hox13-binding sites. Potentially, this implies that the posterior Hox factors are more promiscuous in their interacting partners.

In our previous work, we proposed that Hox13 factors act as ‘guarantors’ of gene expression in NMps, following an idea first proposed from studies in Caenorhabditis elegans by Chalfie and colleagues (Zheng and Chalfie, 2016; Zheng et al., 2015). In this view, Hox13 factors may not be absolutely essential but instead provide robustness to the Brachyury/Wnt autoregulatory loop. Our proposal was based on the observation that zebrafish hoxa13;hoxd13 mutants showed synergistic defects in formation of the body axis and in tbxta expression when Tbxta activity was reduced using a temperature-sensitive mutation, but did not show an effect when Tbxta was fully active (Ye and Kimelman, 2020). How Hox13 proteins might act as guarantors is unclear. One possibility is that they enhance the binding of other factors, such as Tcf7/Lef1 and Rbpj, through collaborative interactions creating what Mann and colleagues called a ‘Hoxasome’ (Mann et al., 2009). Alternatively, or in addition, the Hox13 proteins, and perhaps all of the posterior Hox proteins, could act as pioneer transcription factors that help to open the chromatin at the Hox elements, a hypothesis based on recent studies demonstrating that HOX13 proteins as well as other posterior HOX factors change the chromatin landscape by opening up inaccessible regions (Amandio et al., 2020; Bulajic et al., 2020; Desanlis et al., 2020). Although this has been well documented in the limb bud and genital tubercle in mice, the precise role of Hox13 factors in regulating tbxta enhancer elements awaits further study.

During the gastrula stages, the proximal tbxta promoter is sufficient to drive expression throughout the NMps through the use of two enhancer elements, one of which (E2) contains a binding site for the Nodal-response factor Foxh1 (Chen et al., 1997; Harvey et al., 2010). This proximal promoter also contains sequences necessary to drive tbxta expression in the notochord beyond the gastrula stages, although the enhancer element for this expression is not yet known (Harvey et al., 2010; Row et al., 2016). By the end of gastrulation, we propose that a switch occurs such that the distal Hox-binding elements are required to promote tbxta expression in the NMps, primarily because of Wnt signaling stabilizing β-catenin, which binds to Tcf7/Lef1 factors and activates transcription (Fig. 8). We propose that Rbpj and the posterior Hox factors aid in this activation, with the posterior Hox factors acting to ensure robust activation of tbxta expression.

Elimination of Hox element 1 causes phenotypic defects only in the most posterior somites. This might be because of partial compensation by Hox element 2, given that we observed tbxta expression in the tailbud of the mutants lacking Hox element 1 until fairly late in somitogenesis. However, previous experiments with a temperature-sensitive Tbxta demonstrated that a strong reduction in Tbxta activity from the start of gastrulation onward causes severe phenotypic changes in posterior body formation, whereas the same reduction in Tbxta activity beginning at the end of gastrulation produced only mild phenotypic changes restricted to the most posterior somites (Kimelman, 2016a). This surprising minor requirement for Tbxta function during somitogenesis stages in zebrafish is potentially explained by a study that showed that, although most of the neural and mesodermal cells in the posterior body are allocated from a NMp pool during the gastrula stages, a second NMp population remains resident in the tailbud and contributes only to the most caudal somites (Attardi et al., 2018). It is this second population that may depend crucially on the post-gastrula-expressed tbxta in zebrafish and, thus, on the requirement for Hox element 1 to maintain robust tbxta transcription until the end of somitogenesis. Interestingly, zebrafish have a relatively small number of somites (∼32), whereas other fish, particularly teleosts such as eel, have many more (up to ∼200; Ward and Brainerd, 2007); thus, in these more elongated body plans, the role of NMps and of post-gastrula tbxta expression driven by the upstream Hox elements may have a larger role. Although these other fish species are less studied than zebrafish, it will be interesting to study NMps and the role of Tbxta in these species (see also Sambasivan and Steventon, 2020).

Why do the distal tbxta enhancer elements become required after the gastrula stages? In fish and frogs, Bmp signaling on the ventral side of the embryo is required for formation of the tailbud and posterior body development (Kimelman and Martin, 2012; Tuazon and Mullins, 2015), and regulates the gastrula-stage ventral expression of tbxta in the cells that will become NMps of the posterior body (Northrop et al., 1995), acting through one of the elements (E1) in the proximal 2.1 kb tbxta promoter (Harvey et al., 2010). However, after the gastrula stage, Bmp signaling does not have a major role in axis formation (Pyati et al., 2005). Thus, we propose that, during the gastrula stages, Nodal, Wnt and Bmp act combinatorially to activate tbxta throughout all NMps in the body, but, as gastrulation progresses and Bmp signaling declines, the distal enhancer elements are required to sustain tbxta throughout the somitogenesis stages. Although this model in detail is specific to zebrafish and potentially frog embryos, it will be interesting to examine the role of Rbpj and posterior Hox factors in the regulation of Brachyury in other vertebrates in light of the observation that Rbpj is essential for posterior body formation in mouse (Oka et al., 1995).

WT fish were an AB/WIK mixture and were used for crosses when between 3 months and 3 years of age. All animal protocols were approved by the University of Washington Institutional Animal Care and Use Committee. Embryos for all heat-shock lines were heat shocked in a circulating water bath for 30 min at 40°C, unless otherwise noted.

Amplification of cDNA or genomic DNA used the primers listed below using either Platinum Pfx DNA Polymerase (ThermoFisher) or Q5 DNA Polymerase (New England Biolabs) according to the manufacturer's conditions, with an annealing temperature of 60°C.

The hoxa13b-coding region was amplified from 15-somite stage zebrafish embryo cDNA using the primers GGATCCGCTATGACAGCGTCTTTACTCCTC and GTCAACAAGTACAAGGGCATCAGTTATCGAT, with added BamH1 and ClaI sites shown in italics, and inserted into a vector such that the stop codons were removed and a viral 2A peptide (Provost et al., 2007) was placed immediately after the coding region. The GFP sequence was placed immediately after the 2A sequence. This sequence was placed in a Tol2-hsp70 vector and the resulting plasmid was used together with Tol2 transposase to create stable transgenic lines, as previously described (Kawakami, 2004). This line was designated ‘w249’.

hoxa13b genomic DNA was amplified from the 5′ UTR to a primer that removed the stop codon using the primers GTTGGATCCTGACGCACGC and GTCAACAAGTACAAGGGCATCAGTTATCGAT, with the naturally occurring BamH1 site underlined and the added ClaI site in italics. This was inserted into a vector such that, following the stop codon, was a sequence encoding three FLAG tags, followed by the viral 2A peptide and then GFP. This was then flanked on both sides by the Mbait sequence (Kimura et al., 2014). The plasmid was injected into one-cell-stage embryos together with Cas9 and gRNAs targeting the Mbait sequence and 5′ and 3′ UTRs of hoxa13b to produce a CRISPR/Cas9-mediated KI fish. Fish were screened for GFP expression in the posterior region. This line was designated ‘w250’.

The rbpja-coding region was amplified from 15 s zebrafish embryo cDNA using the primers GGATCCAAGATGGCGCCTGTTGTGACAG and CCTCTGCCATGTCCGTCTCCTATCGAT, with the added BamH1 and ClaI sites shown in italics, and inserted into the same vector used for hoxa13b as described above that had the viral 2A sequence between the GFP and Rbpja. A transgenic line was made with Tol2, and was designated ‘w251’.

An Engrailed repressor domain was inserted in-frame in front of the rbpja-coding region in the HS:rbpja-FLAG-GFP plasmid. A transgenic line was made with Tol2, and was designated ‘w252’.

A 2.1 kb fragment of the tbxta promoter was cloned from genomic DNA using GAATTCATACAATTCCTTTGTGCTGTTGCAACAC, which places an EcoRl I site at the 5′ end, and CCATGGTTCCGATCAAATAAAGCTTGAGATAAGTCCG, which places a Nco1 site at the ATG, and inserted into a vector containing GFP and Tol2 elements. Cloning of the upstream elements used the following primers: for Hox element 2, GTCGACCGTTGTTTAAATAAAACGGCGAGATACATG and GTCAGATATGGACAAACCCATCCATCTCGAG, using Sal1 and Xho1, respectively; for Hox element 1, GTCGACTTTAACCTTGACAAGTGTGAATAGCTG and ATCGATACGAATAATTAACTAATAATTTTGATTTCAACTGTAC, using Sal1 and Cla1, respectively; and for the element at −5.4 kb ATCGATTCCTAAAGGTCTGACATGTACTGCG and GAATTCGGTTTTTCATTGTAGTTAACTGTGAGAGC, using Cla1 and EcoR1, respectively. This line was designated ‘w253’.

The 5.7 kb tbxta promoter was cloned with GTCGACCGTCAACCTAGACACATTTAAATTTGGC, which places a Sal1 site at the 5′ end, and CCATGGTTCCAGTCTTATTGGGGGAAAAGCC, which places a Nco1 site at the ATG, and inserted into a vector containing GFP and Tol2 elements.

These were produced using the gRNAs listed below together with Cas9. Tbxta^Δ312 has a 312 bp deletion (AAAGCTTGTCCTGTAGGGGGTGCTATΔGATCCTATACTGTGCTCCAGACTCCACAGA), whereas Tbxta^Δ314 has a 323 bp deletion and 9 bp insertion (insertion in italics) (GCTTGTCCTGTAGGGGGTΔCTGCTGGCCΔTCCTATACTGTGCTCCAGA). In both cases, the Δ indicates the region deleted. Embryos with the deletions were screened using the primers GCTGCACCCAAGAAAAGCAA and GTCACTTTGTTCAACTGTAGCGT using standard PCR conditions, and analyzed on 2% agarose gels.

The gRNAs used to make the KI line were GAGCTGCTGGGCTCCATGTA and GGGCTTGATATTGGTGGTAT. The Mbait gRNA was GGCTGCTGCGGTTCCAGAGG. The tbxta gRNAs used to remove Hox element 1 were CCTGTAGGGGGTGCTATCTC and GAGCACAGTATAGGATCTGG.

The basic CUT&RUN procedure followed that described by Skene and Henikoff (2017) except that a pAG/MNase (Meers et al., 2019) was used instead of the original pA/MNase. Briefly, zebrafish tailbuds were dissected at the boundary of the third newly formed somite from groups of 30 embryos at the 15- to 18-somite stage, with the epidermis first removed prior to cutting, as previously described (Manning and Kimelman, 2015). Explants were then dissociated into single cells with wash buffer [20 mM HEPES, pH 7.5, 150 mM NaCl, 0.5 mM spermidine and a protease inhibitor (Roche; one complete tablet per 50 ml added fresh)] by gentle pipetting up and down 50 times using a 1 ml pipette. Tubes were precoated with 1% bovine serum albumin overnight at 4°C to avoid cells sticking to the wall. The dissociated single cells were then washed twice, bound to activated Concanavalin A-coated magnetic beads, and suspended in 250 µl antibody buffer [wash buffer containing 0.03% (wt/vol) digitonin (Dig) and 2 mM EDTA]. Antibody was added immediately to a final concentration of 1:100 and cells were incubated overnight at 4°C on a nutator. After incubation, beads/cells were washed three times with Dig-wash buffer [wash buffer containing 0.03% (wt/vol) Dig] to remove unbound antibody, followed by resuspension and incubation of beads/cells in 150 µl Dig-wash buffer with pAG/MNase (1:100 dilution), at 4°C for 1 h. After two washes with Dig-wash buffer, beads/cells were resuspended in 100 µl Dig-wash buffer and chilled to 0°C. pAG/MNase was activated by adding 2 µl 100 mM CaCl₂ and the cutting reaction was allowed to continue for 1 h at 0°C. An equal volume (100 µl) of 2× stop buffer (340 mM NaCl, 20 mM EDTA, 4 mM EGTA, 0.05% Dig, 100 µg/ml RNase A and 50 µg/ml glycogen) was added to stop the reaction, followed by 30 min incubation in a 37°C water bath to release the cut DNA fragments from the cells. Beads/cells were placed on a magnetic stand and the liquid containing the cut DNA fragments was transferred to a new tube. DNA was extracted using the phenol-chloroform extraction method (Skene and Henikoff, 2017).

The CUT&RUN experiment was performed in duplicate for HS:hoxa13b-FLAG-GFP and histone mark H3K4me1 and H3K27ac, and in triplicate for KI:hoxa13b-FLAG-GFP. WT embryos with the anti-FLAG antibody were included as the control for each replicate CUT&RUN experiment. The antibodies used were: anti-FLAG (Millipore F3165), anti-H3K4me1 (Abcam ab8895) and anti-H3K27ac (Abcam ab4729).

Libraries for next-generation sequencing were prepared using the NEBNext Ultra II DNA Library Prep Kit, NEB (Ipswich, USA) for transcription factors following the protocol of Nan Liu (Harvard University; https://dx.doi.org/10.17504/protocols.io.wvgfe3w), with some modifications for histone marks because of larger CUT&RUN fragment sizes (Meers et al., 2019). Modifications for histone mark samples included: (1) in step 2, the PCR program for End Prep was 20°C for 30 min, 58°C for 45 min and hold at 4°C; (2) in step 8, 1.5×AMPure XP bead size selection was applied after adaptor ligation; and (3) in step 19, 0.65×AMPure XP bead size selection was applied for the first round of size selection to remove PCR products larger than 500 bp.

Libraries were pooled at equimolar concentrations and paired-end (PE) 150 bp sequencing was performed on a Illumina HiSeq 4000 platform by Novogene Corporation. Bioinformatics analysis was facilitated by using the advanced computational, storage and networking infrastructure provided by the Hyak supercomputer system at the University of Washington. Paired reads were quality checked and trimmed using Trim Galore version 0.4.4_dev (www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Only reads longer than 20 bp after trimming were kept. Trimmed PE reads were aligned to the zebrafish genome (danRer11 masked) using Bowtie2 version 2.3.4.1 (Langmead and Salzberg, 2012) with options: local–very-sensitive-local–no-unal–no-mixed–no-discordant–phred33. The resulting bam files containing the alignment information were converted to bed format using Samtools 1.9 (Li et al., 2009) and Bedtools 2.29.2 (Quinlan and Hall, 2010), or bigwig format using a custom R script (R version 3.6.2) and visualized with the Integrated Genomics Viewer (IGV version 2.7.2; Robinson et al., 2011). Assessments of the mapped fragment size distribution, replicate reproducibility and heatmap of reads over the peak regions were plotted following the scripts of Henikoff et al. (2020). Fragments of 24-121 bp and 149-500 bp were used for peak calling and downstream analysis of transcription factors and histone marks, respectively. MACS2 version 2.1.1.20160309 (Zhang et al., 2008) was used for peak calling with the parameters: callpeak -t input_file -q 1e-2 -f BEDPE -keep-dup all -n output_file_name. The broadpeak option was used for histone mark peak calling.

Peaks overlapping between the duplicates of each group or between the triplicates for KI:hoxa13b-FLAG-GFP were kept and used for subsequent analysis. The overlap analysis was performed using the Bedtools merge function with flags -c 1,4 -o count,collapse. Features merged from two or more experiments were retained for further analysis. The resulting peaks of the HS:hoxa13b-FLAG-GFP sample were thought to be the direct targets of Hoxa13b and they were further annotated by overlapping with the KI:hoxa13b-FLAG-GFP peaks to add another level of confidence, because these two lines have different advantages in terms of Hoxa13b target profiling. The potential target genes near the peaks were annotated based on the TxDb.Drerio.UCSC.danRer11.refGene annotation database using a custom R script, and ChIPseeker (Yu et al., 2015) was used for genome ontology analysis to generate Fig. 3A and Fig. S2A.

Sequences of the ±200 bp center region of the Hoxa13b peaks were used for comprehensive motif analysis including de novo motif discovery (MEME and DREME), similarity to known motifs (Tomtom) and central enrichment analysis (CentriMo) using the MEME-ChIP program, which integrates the above-mentioned programs and automatically groups significant motifs by similarity (Machanick and Bailey, 2011).

The top-rated motif discovered was the Hoxa13b motif; therefore, the footprint analysis was performed using this motif. Single base-resolution cutting probability of the MNase around the peak regions containing the Hoxa13b motif was calculated using a custom R script (R version 3.6.2) and plotted using a R package ggplot2 (Wickham et al., 2016) to reveal the general footprint of Hoxa13b.

Standard conditions were used for alkaline phosphate in situ hybridization (https://wiki.zfin.org/display/prot/Thisse+Lab+-+In+Situ+Hybridization+Protocol+-+2010+update), whereas the Lauter protocol was used for FISH (Lauter et al., 2011). Immunofluorescence was performed using the MF20 antibody (Developmental Studies Hybridoma Bank) for somite muscle staining and the anti-FLAG antibody (Millipore) for revealing Hoxa13b-FLAG expression in the KI:hoxa13b-FLAG-GFP embryos. A standard immunofluorescence stain protocol was used. Briefly, fixed embryos stored in methanol were rehydrated into PBST (phosphate-buffered saline with 0.1% Tween-20) and blocked for 1 h at room temperature, then incubated with the first antibody (1:50 dilution for MF20 and 1:200 for anti-FLAG) overnight at 4°C. After three washes with PBST (15 min each time), embryos were incubated with the secondary antibody (Abcam ab150113, 1:1000 dilution) for 2 h at room temperature, then washed with PBST and stepped into 70% glycerol for later imaging.

Brightfield images for in situ hybridization and for body-length measurement were taken using a Nikon AZ100 microscope with white-light illumination. Wholemount fluorescence images (first and second rows in Fig. 5B-E) were taken with the same microscope using a fluorescence illumination. Confocal images (Fig. 1C-F, Figs S5, S6 and the magnified images in the third row of Fig. 5B-E) were obtained using an Olympus Fluoview 1200 scanning confocal microscope with a 40× oil lens for fixed embryos mounted on a slide or a 40× dipping lens for live embryos. The 3D image stacks were taken at a step-size of 2 μm. 3D reconstruction of the stacked images was done using Imaris 9.2 (Oxford Instruments) and the posterior somite and neural tube volume data were generated by Imaris after manual creation of objectives based on the fluorescent signal.

One-way ANOVA was used for the body length comparison, and Tukey's HSD test was used for comparisons of the somite volume and neural tube volume shown in Fig. 6G-I.

We thank Derek Janssens and Steve Hennikoff for kindly supplying Protein AG/MNase and valuable technical advice; Colin Kenny and Hannah Arbach for helpful discussions; Lauren Saunders for help with Tapestation analysis; Alex Chitsazan for bioinformatics help; and Rex Dunham and Baofeng Su for providing catfish genomic DNA.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.199408