An ancestral Wnt-Brachyury feedback loop and recruitment of mesoderm-determining target genes revealed by comparative Brachyury target screens

Transcription factors are crucial drivers of cellular differentiation during animal development and often share ancient evolutionary origins. The T-box transcription factor Brachyury plays a pivotal role as an early mesoderm determinant and neural repressor in vertebrates; yet, the ancestral function and key evolutionary transitions of the role of this transcription factor remain obscure. Here, we present a genome-wide target gene screen using ChIP-seq in the sea anemone Nematostella vectensis, an early branching non-bilaterian, and the sea urchin Strongylocentrotus purpuratus, a representative of the sister lineage of chordates. Our analysis reveals an ancestral gene regulatory feedback loop connecting Brachyury, FoxA and canonical Wnt signaling involved in axial patterning that predates the cnidarian-bilaterian split about 700 million years ago. Surprisingly, we also found that part of the gene regulatory network controlling the fate of neuromesodermal progenitors in vertebrates was already present in the common ancestor of cnidarians and bilaterians. However, while several neuronal Brachyury target genes are ancestrally shared, hardly any of the key mesodermal downstream targets in vertebrates are found in the sea anemone or the sea urchin. Our study suggests that a limited number of target genes involved in mesoderm formation were newly acquired in the vertebrate lineage, leading to a dramatic shift in the function of this ancestral developmental regulator.

Introduction T-box genes are a class of transcription factors that evolved in the common ancestor of protists and animals and have been found to play important roles in the development of all animals. The founding member, Brachyury, also called T in mice (and more recently TbxT), is a key determinant of mesoderm in vertebrates (Cunliffe and Smith, 1992;Herrmann et al., 1990;O'Reilly et al., 1995;Smith et al., 1991;Wilson et al., 1995). It is expressed in the presumptive mesoderm and continues to be expressed after gastrulation in the notochord, a de ning structure of all chordates (Kispert et al., 1995;Smith et al., 1991;Stemple, 2005;Wilkinson et al., 1990). Heterozygous mouse mutants of the T gene locus show a shortened tail (hence the Greek name Brachy = short, ury = tail), while homozygous mutants die in utero lacking most of the body posterior to the forelimbs (Beddington et al., 1992;Stott et al., 1993). Knockdown of brachyury in Xenopus also impairs the posterior development of the embryo, while overexpression of brachyury mRNA in explanted animal caps leads to mesoderm differentiation (Conlon and Smith, 1999;Cunliffe and Smith, 1992;O'Reilly et al., 1995;Smith et al., 1991). Similar roles have been found in other vertebrates (Kispert et al., 1995;Schulte-Merker et al., 1994), demonstrating the conserved role for Brachyury as a pivotal mesoderm regulator across vertebrates.
During anterior-posterior axis elongation of the vertebrate embryo, brachyury expression is mainly restricted to the most caudal part of the embryo and the notochord (Smith et al., 1991;Yamaguchi et al., 1999). This caudal end of the developing embryo contains a population of bipotent cells called neuromesodermal progenitor cells (NMPs) (Davis and Kirschner, 2000;Gouti et al., 2014;Henrique et al., 2015;Koch et al., 2017;Selleck and Stern, 1991), which contribute to the neuronal differentiation of the spinal cord, as well as the paraxial mesoderm and notochord development (Gouti et al., 2017;Tzouanacou et al., 2009). Some studies de ne these NMP cells by the expression of brachyury and sox2, a member of the soxB1 gene family (Gouti et al., 2014;Koch et al., 2017), while others argue that brachyury expression along with Wnt/FGF signaling is su cient for the cells to become NMPs even in the absence of Sox2 protein (Mugele et al., 2018). Regardless of the de nition of NMPs, the mesodermde ning role of Brachyury has been established within these cell populations. The notochord-speci c function of Brachyury also appears to be conserved in non-vertebrate chordates. In the ascidian Ciona intestinalis, representing the sister group to the vertebrates, brachyury expression and function is restricted to the notochord and brachyury knockdown abolishes notochord development (Takahashi et al., 1999). In amphioxus, there are two brachyury paralogs. While to date no functional data are available for amphioxus, the expression pattern of bra2 is very similar to that in vertebrates, suggesting that it likely plays a similar role in this cephalochordate (Holland et al., 1995), while bra1 appears to be more restricted to the notochord (Satoh et al., 2021). Thus, Brachyury is key for mesoderm formation and in particular for notochord development in vertebrates, and, possibly, in chordates. However, whether Brachyury has a mesodermal function outside chordates is less clear.
Brachyury has also been studied in non-bilaterian animals, including cnidarians. For instance, in the freshwater polyp Hydra, the brachyury homolog Hybra1 is expressed at the hypostome and acts as an early marker of head formation during budding, regeneration and embryonic development (Technau and Bode, 1999;Technau and Scholz, 2003;Technau et al., 2000). While hydrozoan cnidarians do not form a blastopore during gastrulation, the area giving rise to the hypostome is equivalent to the blastoporal region in other cnidarian clades. In the anthozoan sea anemone Nematostella vectensis, brachyury is expressed at the blastopore margin . Functional studies in Nematostella vectensis and the coral Acropora digitifera suggest that brachyury may be involved in the formation of the pharynx and the endo-ectoderm boundary (Servetnick et al., 2017;Yasuoka et al., 2016). Strikingly, when Hybra1 mRNA is injected into blastomeres of the animal pole of Xenopus embryos, it is capable of inducing mesoderm formation with high e ciency, showing that in the context of the vertebrate, the cnidarian transcription factor behaves like its cognate vertebrate homolog (Bielen et al., 2007a;Marcellini et al., 2003). Thus, the protein is conserved enough to mimic the function of the vertebrate Brachyury and to act as a mesoderm determinant. These observations taken together raise the question: How were the genomic targets and the gene regulatory network downstream of Brachyury altered in the course of animal evolution, in order for this transcription factor to acquire such diverse developmental roles? Therefore, to shed light on the evolution of Brachyury function, genome-wide surveys of its downstream genes in taxa with evolutionarily informative phylogenetic positions are needed (Fig. 1A).
To reconstruct the evolutionary changes in the Brachyury gene regulatory network, we carried out genome-wide ChIP-seq and identi ed Brachyury target genes in two invertebrate species with key phylogenetic placements: the non-bilaterian, diploblastic sea anemone Nematostella vectensis and the deuterostome sea urchin Strongylocentrotus purpuratus, belonging to the sister group of chordates. We revealed an ancestral developmental kernel conserved between the cnidarians, echinoderms and vertebrates that includes a feedback loop between Brachyury, FoxA and canonical Wnt signaling. This circuit is involved in the axial patterning of the primary body axis. By contrast, only very few homologous genes that are crucially involved in vertebrate mesoderm formation and differentiation are also Brachyury targets in Nematostella or Strongylocentrotus. Notably, we did nd a signi cant number of neuronal target genes shared between Nematostella and vertebrates, many of which are repressed as revealed in brachyury knockdowns in both phyla. This suggests that axial patterning and suppression of neuronal targets may be the ancestral function of Brachyury and the acquisition of target genes conveying the role of Brachyury in mesoderm formation emerged within chordates.

Materials And Methods
Animal culture Nematostella polyps were kept as previously described (Fritzenwanker and Technau, 2002;Genikhovich and Technau, 2009a) at 18°C and largely in the dark. Spawning was induced by a combination of light and elevated temperature (Fritzenwanker and Technau, 2002). Embryos were raised at 21°C. Adult Strongylocentrotus individuals were kept in circulating seawater aquaria at Stazione Zoologica Anton Dohrn in Naples. Spawning was induced by vigorous shaking of gravid animals. Embryos were raised at 15°C in ltered Mediterranean Sea water (FSW) diluted 9:1 with deionized water (Andrikou et al., 2013;Anishchenko et al., 2018;Leahy, 1986;Perillo et al., 2018) Antibody generation A 6xHis-tagged full-length clone of Nematostella Brachyury was expressed in E. coli strain BL21. After transformation, bacteria were grown overnight in LB supermedium + Ampicillin. Then, 10 ml LB medium / Ampicillin were inoculated with 500 µl of overnight culture, grown for 3 h (O.D. 600 = 0,5-0,6), IPTG was added to a nal concentration of 1 mM, and bacteria culture was incubated for 21 h at 4°C. Brachyury protein was puri ed using Ni-NTA agarose (Qiagen), followed by polyacrylamide gel electrophoresis. The Brachyury band at ~50 kDa was excised and used to inject two rabbits (PRIMM Biotech). For Strongylocentrotus Brachyury, two polyclonal antibodies from 2 different rabbits were raised against the recombinantly produced C-terminal domain excluding the T-box domain. The antibodies were a nity puri ed by PRIMM Biotech prior to their use.
For both Nematostella and Strongylocentrotus the reactivity of polyclonal antibodies from each rabbit were assessed by both Western blot and immunohistochemistry, con rming their speci city (Supplemental Fig. 1). For the ChIP replicates, antibodies from the two different rabbits was used for each.
Chromatin Immunoprecipitation and Library preparation ChIP-Seq was carried out essentially as previously published (Schwaiger et al., 2014). In brief, embryos at gastrula stage were xed with 2% formaldehyde for 12 min. Nuclei were collected by Dounce homogenizing the embryos. Chromatin was fragmented to an average size of 100-200bp using a Covaris S1 ultrasonicator with the following settings: Duty Cycle: 20%, Intensity: 5, Cycles/Burst: 200, Duration: 240s, Mode: Frequency Sweeping. To carry out the immunoprecipitation, 300 µg protein (as quanti ed by Bradford assay) of the fragmented chromatin was incubated with 5 µg of Nematostella or Strongylocentrotus anti-Brachyury antibodies. Chromatin from ~20 immunoprecipitations was pooled to obtain ~2ng of DNA, as measured by a Pico Analyzer. This was then used for library preparation according to Illumina ChIP-seq DNA Sample Prep Kit instructions (Catalog number IP-102-1001). The quantity and quality (size distribution) of the libraries was con rmed using an Agilent Bioanalyzer. Deep sequencing was performed at the Vienna BioCenter Core Facilities (VBCF) (https://www.viennabiocenter.org/facilities/) with 50 bp SE HiSeq 2500 reads.

Ortholog detection
For ortholog detection we used the OMA database (Altenhoff et al., 2013(Altenhoff et al., , 2018. Before we subjected our in silico translated protein sequences to OMA, we excluded peptide sequences smaller than 100 amino acids. To detect the orthologs over the wide evolutionary distance of our investigated organisms, it was necessary to use a less stringent 'Length tolerance ratio' = 0.2 parameter of OMA. The other OMA parameters used default values.

ChIP-seq Analysis
In addition to our ChIP-seq datasets for Brachyury in Nematostella and Strongylocentrotus we acquired raw reads of published Brachyury ChIP-seq datasets of mouse (Koch et al., 2017) (GEO: GSE93524) and Xenopus (Gentsch et al., 2013) (GEO: GSE48560). Reads from M. musculus, X. tropicalis, S. purpuratus, N. vectensis were mapped to genome versions GRCm38 (mm10), v9.1, v5.0, v1.0, respectively, with BWA using the BWA-MEM algorithm (Li and Durbin, 2010). ChIP-seq peaks were called using Peakzilla (Bardet et al., 2013) using Input sequencing as control. Depending on the bimodal distribution of the reads, the Peakzilla algorithm assigns a score to each peak. We examined the relationship between Peakzilla peaks and their scores, and the coverage of ChIP-seq reads across the genomes, thereby determining the optimal score cutoff to lter high-quality peaks in each dataset, where the lowest scoring peaks that pass the cutoff are still visually detected as peaks in the coverage tracks. Consequently, we removed peaks with a score less than 2 for Mouse and Xenopus, less than 0.5 for Strongylocentrotus and less than 5 for Nematostella.

Target gene selection
To identify target genes, a commonly used procedure is to select the closest gene to an identi ed peak. However, in more compact genomes like Nematostella or Strongylocentrotus (5-10 genes per 100 kb) peaks in the intergenic space are often only marginally less distant to a gene in the other direction.
Therefore, we identi ed the two closest genes to the binding site. Next, using ortholog information obtained earlier with OMA (Altenhoff et al., 2013(Altenhoff et al., , 2018, we determined whether either of the two closest genes has an ortholog that is a target in any other species considered in this study. In about 20% of Nematostella, 24% of Strongylocentrotus, 3% of Xenopus and 28% of mouse peaks, no ortholog information was found, thus, the closest gene to the binding site was assumed to be the target. If both closest genes had orthologous targets, we selected the one with larger number of orthologous target genes. In case an equal number of the orthologous target genes was found, we kept both putative targets and marked them with ++ in the Supplementary Tables 1-4 (Supplemental Fig. 1G).

Motif analysis
To nd and compare the binding motifs of Brachyury we used the MEME-ChIP (4.11.2) suite (Machanick and Bailey, 2011). For the number of motifs to be predicted in a sequence we selected 'anr ' (Any Number of Repetitions) mode and 1e-10 for the e-value cut-off criterion. The discovered motifs were scanned against the non-redundant sequence pro les obtained from the JASPAR database (Khan et al., 2018). The other parameters were set to default.

Co-occupancy of motifs
To investigate what other motifs are present in the vicinity of Brachyury binding motif and their spatial arrangement, we made use of MEME-ChIP with an in-house python script, which nds all the motifs in a peak and their absolute distances from the center of the peak. It also groups the motifs into motif groups based on motif families in JASPAR (Khan et al., 2018). In case of overlapping motifs, it selects the motif with better FDR adjusted p-value assigned by Fimo from the MEME-ChIP (Machanick and Bailey, 2011) suite.
To knockdown Strongylocentrotus brachyury, we injected 2-4 pl of a 200 µM solution of Morpholino previously characterized with the sequence: CGCTCATTGCAGGCATAGTGGCG (Annunziata and Arnone, 2014). We quanti ed the mortality at 24 hr after fertilization and discarded all experiments with <80% survival. Embryos were either xed for WMISH or used for isolation of RNA extraction at 24 hr. PolyAenriched RNA samples were collected and processed for library production using the Lexogen kit and sequencing (50bp single-end HiSeqV4).

Differential gene expression analyses
For Nematostella, we generated 6 biological replicates of each translational and splice morpholino. Principal component analysis (PCA) showed that the replicates cluster together (Supplemental Figure  2D). For Strongylocentrotus, we generated three biological replicates (Supplemental Figure 2E). This was compared to three replicates of Morpholino knockdown in Xenopus (Gentsch et al., 2013) Figure 2F) and two replicates of mouse knockout mutants (Koch et al., 2017) Figure 2G).
Reads aligned to genomic features were counted with featureCounts (Subread v1.6.2) (Liao et al., 2013(Liao et al., , 2014. The counts les generated were used as input to SARTools (v1.6.0) (Varet et al., 2016), which is a wrapper R script around two widely used differential expression analysis tools, namely DESeq2 (Love et al., 2014) and edgeR (McCarthy et al., 2012). Genes with Log2fold change of 1 and FDR corrected p-value (q-value) of 0.1 are considered as differentially expressed. The RNA-seq reads for Strongylocentrotus were quasi-mapped to the Strongylocentrotus purpuratus genome v.5 associated transcriptome using Salmon (v0.11.3) (Patro et al., 2017). The resulting count les were used in DESeq2 (Love et al., 2014) for differential expression analysis, differentially expressed genes with p-adjusted value of less than 0.05 were considered signi cant.
Whole mount in situ hybridisation (WMISH) WMISH in Nematostella was carried out as described before (Genikhovich and Technau, 2009b) with a few minor modi cations. Brie y, embryos were xed in 4% PFA in Phosphate-Buffered Saline (PBS) for 1 h at RT, and then washed 5x with 0.2% Tween in PBS at RT. The embryos were then washed a single time with 50% methanol in PBS with 0.2% Tween and lastly transferred into 100% methanol and stored at -20°C until further processing. Hybridization was performed at 1 ng riboprobe/µl, overnight at 63°C, and detection was used anti-digoxigenin (1:2000) antibody (Roche). If no strong signal and mostly background began to develop in the NBT/BCIP solution, the embryos were washed 5x in PTw and then left to stand in PBS with 2% w/v triton-100x until the background was at acceptable levels. This was repeated until a strong signal was observed.
In Strongylocentrotus, WMISH was done as previously described (Perillo et al., 2021). Brie y, embryos were xed in 4% PFA in MOPS Buffer overnight at 4°C, then gradually dehydrated to 70% ethanol and kept at -20°C until use. Fixed embryos were gradually rehydrated in MOPS buffer, washed several times with the same medium and were pre-hybridized for 3 hours in the hybridisation buffer. Hybridization was carried out with 0.1 ng/μl of probes and incubated for 1 week. The signal for colorimetric in situ hybridisation was developed using anti-digoxigenin AP (Alkaline Phosphatase) conjugated antibody and the AP substrate, while uorescent in situ hybridization was developed with uorophore-conjugated tyramide (1:400 reagent diluents, Perkin Elmer).

Immunohistochemistry
Immunohistochemistry in Nematostella was essentially performed as described in (Genikhovich and Technau, 2009b) with the following changes. Brie y, embryos were xed at 4°C in 4% PFA in PTwTxD (1xPBS with 0.2% Tween, 0.2% Triton X-100 and 0.2% DMSO). All further washing steps were performed with PTwTxD. They were then washed 1x with PTwTxD and ice-cold acetone and put on ice. Once the embryos had settled down the acetone was removed, and they were washed 10x times with the cold PTwTxD. The embryos were then blocked 2 h at RT in the blocking solution containing 20% v/v of heat-inactivated sheep serum and 80% of 1% w/v BSA in PTwTxD. The blocking solution from the embryos was substituted with the rabbit anti-Bra preadsorbed in the bocking solution (1:500) for the time of the blocking and incubated overnight at 4°C on a table rocker. The next day, the embryos were washed 10x with cold PTwTxD and blocked once again for 2 h at RT in the blocking solution. Blocking solution was then replaced with preabsorbed goat anti-rabbit Alexa Flour 568, 1:2000, Invitrogen) mixed with DAPI ( nal concentration 5 µg/ml) and Alexa Fluor 488 phalloidin ( nal concentration 5 U/ml, Invitrogen) and incubated overnight at 4°C. Next day, the embryos were washed 10x in PTwTxD, in ltrated overnight with Vectashield (Vector labs), and imaged.
Immunohistochemical detection of Brachyury in Strongylocentrotus was performed as described in (Perillo et al., 2021). Brie y, the embryos were xed in 4% PFA in ltered sea water for 15 min at room temperature (RT), dehydrated in 100% ice cold methanol for 1 min and washed several times with 1x PBS. PBS was replaced with a blocking solution (4% Sheep Serum, 1% BSA in 1x PBS) and the samples incubated for 1h at RT. The blocking solution was replaced with the solution containing the primary antibodies diluted in blocking solution (Brachyury 1:100, Nkx2.1 1:600, SoxB2 1:500) and incubated either overnight at 4°C or 90 min at 37°C. Embryos were washed several times with PBS, the secondary antibody goat anti-rabbit AlexaFluor 488 (Invitrogen) diluted in PBS was added and then washed again several times.

Phylogenetic analysis
Phylogenetic analysis of T-Box, Sox, Zic, TFAP, NR2, RNF, Rfx family members was performed using maximum likelihood methods. Orthologs and paralogs were obtained from the respective proteomes using BLASTP with the protein sequences from Nematostella and mouse as query. The peptide sequence alignments were generated using MAFFT with the parameters -maxiterate 1000 -localpair (Katoh and Standley, 2013) and cut to the unambiguously aligned core region. From the alignments maximum likelihood trees were reconstructed using IQ-TREE 1.6.12 Substitution models of LG+I+G4 (T-Box and Sox) and LG+F+G4 (Zic, TFAP, NR2, RNF, Rfx) were chosen after running model selection in IQ-TREE. Support values were generated using UFboot (Minh et al., 2013) as implemented in IQ-TREE with 1000 (T-Box and Sox) and 10000 (Zic, TFAP, NR2, RNF, Rfx) bootstrap samples, respectively.

Genome wide detection of Brachyury binding sites in Nematostella and Strongylocentrotus
In order to reconstruct the evolution of Brachyury function, we analyzed its genomic targets in representatives of phylogenetically informative groups (Fig. 1A), the sea anemone Nematostella and the sea urchin Strongylocentrotus. We rst generated antibodies against Nematostella and Strongylocentrotus Brachyury proteins (see Materials and Methods, Supplemental Fig. 1A-D). Antibody stainings of gastrula stage embryos of Nematostella confocal microscopy con rmed that the antibodies detect nuclei at the ectodermal margin of the blastopore, where mRNA expression is also detected by whole mount in situ hybridisation (WMISH) (Fig. 1B). Similarly, in Strongylocentrotus gastrulae, Brachyury protein is detected at the endodermal margin of the blastopore as well as at the future stomodaeum, where the archenteron will break through, re ecting the mRNA expression (Fig. 1C). Loss of function experiments by morpholino knockdowns and localized gain of function experiments con rmed the speci city of the antibodies in Nematostella (Supplemental Fig. 1A-C) and Strongylocentrotus (Supplemental Fig. 1D) (Annunziata and Arnone, 2014).
We next used these antibodies to perform ChIP-seq in early gastrula stage embryos. For Nematostella, we identi ed 2389 putative binding sites in two highly reproducible replicates, which mapped to 1543 putative target genes ( Fig (Schwaiger et al., 2014) as de ned by the combination of p300, H3K27ac, and H3K4me1, suggesting that about 1/3 of the identi ed enhancers are bound by Brachyury at the gastrula stage (Supplemental Fig. 2A,B). Similarly, a large fraction (338/490) of the Brachyury binding sites detected by ChIP-seq in Strongylocentrotus are found in open chromatin as detected by ATAC-seq (Supplemental Fig. 2C). However, we also found an even larger fraction in both species that are not associated with active enhancers or open chromatin (see Discussion). To compare our ndings with other species, we re-analysed previously published ChIP-seq data (Gentsch et al., 2013;Koch et al., 2017) using our peak calling thresholds and found for mouse and Xenopus 4000 and 2497 peaks, respectively. Corroborating the published results (Gentsch et al., 2013;Koch et al., 2017), our called mouse peaks corresponded to 3060 putative target genes, while in Xenopus 1376 target genes were detected (Fig. 2B).
When analysing the distribution of the Brachyury binding sites with respect to target genes, we found that the distribution re ects the genome sizes of the species. In species with larger genome size the peaks are distributed over a larger range, sometimes beyond 100 kb while in species with smaller genome size the peaks are much more in the vicinity of transcription start site (TSS) (Fig. 2H-I). This is particularly obvious in the case of Nematostella, where about 50% of all Brachyury binding sites are located within 1kb upstream or downstream of the TSS, compared to all three deuterostome species (5-12%) (Fig. 2D-F). The smaller genome size also make the gene density in Nematostella higher (about 1 gene / 10 kb) than other species concerned (Putnam et al., 2007). Different transcription factors tend to have strikingly different positional speci city within the enhancer regions (Grossman et al., 2018). From the binding pro le of Brachyury ( Fig. 2D-G), Brachyury seems to be strongly concentrated around the TSS site in all the species under consideration.
Brachyury binding motifs are highly conserved The Brachyury DNA binding motif was originally found to be a palindromic sequence by an in vitro Selex approach Müller and Herrmann, 1997), which was supported by ChIP-seq in Xenopus embryos (Gentsch et al., 2013). We detected 700 and 238 palindromic Bra binding motifs in the sea anemone and sea urchin Brachyury ChIP-seq peaks, respectively. This indicates that the ability for dimerization is an ancestral feature of Brachyury ( Fig. 2B-C). We detected signi cantly higher Brachyury peaks with palindromes compared to peaks with half palindromes (Supplemental Fig. 3A), suggesting a higher binding a nity of the dimer compared to the half palindromes. However, we found no signi cant correlation between the distribution of half-palindromes or palindromes with the expression level or GO category of the target genes (data not shown). Thus, the biological signi cance of the distribution of halfpalindromes versus palindromes remains obscure at this point. However, the detailed comparison of the binding motifs of Brachyury in Nematostella, Strongylocentrotus, Xenopus and mouse showed that the canonical binding motifs are highly conserved throughout metazoans (Fig. 2C). Similar motifs have also been found among the ATAC-seq peaks in the protist Capsaspora (Sebé-Pedrós et al., 2016). This suggests that the DNA-binding domain and the corresponding binding motif have been conserved from protists to humans, suggesting a strong positive selection. This might explain why the overexpression of the cnidarian or even Capsaspora brachyury mRNA in frogs induces the formation of mesoderm (Bielen et al., 2007b;Marcellini et al., 2003;Sebe-Pedros et al., 2013).
In order to detect potential co-factors or competitors of Brachyury in the four species, we searched the peaks for an enrichment of other common motifs, also considering their position with respect to the peak summit. As expected, Brachyury motifs locate centrally, while other transcription factor binding motifs are enriched within ~50bp from the peak center, suggesting possible co-binding (Supplemental Fig. 3B).
When scanning the peak sequences with FIMO for known motifs matching to the enriched motifs identi ed in any species, we found a similar distribution of putative binding sites for homeodomain, bHLH, Pax, HMG (Sox) and Fox proteins in all species (Supplemental Fig. 3B). Especially in Nematostella, a substantial number of putative Sox, Fox and Hox motif sites are found together with Brachyury (Supplemental Fig. 3B). This is of particular interest, as the expression of foxA, foxB, as well as soxB1 and soxB2 is partially overlapping with brachyury in Nematostella (see below).
Brachyury can act as activator or repressor of developmental regulator genes in both the sea anemone and sea urchin Next, we wished to study the role of Nematostella and Strongylocentrotus Brachyury on gene expression by knockdown of Brachyury function. In Nematostella, we generated 4, 5 and 6 transcriptomes of gastrula stage embryos each injected with control, splicing or translation blocking morpholinos, respectively (Supplemental Fig.2D). In Strongylocentrotus, we used a previously published and validated translation blocking morpholino (Annunziata and Arnone, 2014). A principal component analysis showed that splice morpholino and translation morpholino transcriptomes cluster together (Supplemental Fig. 2D, E), distinct from the control Morpholino experiments, without batch effects. Differential gene expression analysis revealed 85 differentially expressed ChIP target genes (corresponding to 23 % of ChIP targets) in bra knockdowns (44 down-regulated and 46 up-regulated) in Strongylocentrotus (Supplemental Fig.3C). In Nematostella, the overlap between the ChIP target and differentially expressed genes is 90 genes (6% of ChIP-targets), of which 38 are down-regulated and 52 are up-regulated (Supplemental Fig. 2C). We compared these data with transcriptomes comparing knockouts or knockdowns from mouse and Xenopus, respectively (Supplemental Fig. 2F-H). Of 3060 mouse putative target genes, 108 (3.5% of ChIPtargets) are differentially expressed in brachyury mutants, among them 81 downregulated and 27 are upregulated (Koch et al., 2017). By comparison, in Xenopus, of 1376 Brachyury target genes 73 (5.3% of ChIP-targets) are differentially expressed in Brachyury knockdowns (40 are down-regulated and 33 upregulated) (Gentsch et al., 2013). Thus, the proportion of genes regulated in relation to ChIP target genes is in line with earlier published reports of 1% to 10% (Farnham, 2009;Scacheri et al., 2006;Yang et al., 2006).

Gene Ontology (GO) analyses of the Brachyury ChIP targets found in Nematostella and
Strongylocentrotus shows an enrichment in categories of transcriptional regulation, Wnt signaling, multicellular organism development, signal transduction and biological regulation (Supplemental Fig 3C), which is similar to the GO categories of Brachyury ChIP targets reported in vertebrates (Gentsch et al., 2013;Koch et al., 2017;Tosic et al., 2019). This suggests that in these species Brachyury is a developmental regulator that mainly acts positively or negatively to control other developmental genes.
Brachyury forms an ancestral gene regulatory loop with Wnt signaling and activates or represses key blastoporal genes When analysing the target genes in Nematostella and in Strongylocentrotus, we found that many members of Wnt, FGF, Notch and BMP signaling pathways as well as several transcription factors were among the direct target genes ( Fig. 3; Fig.6, Supplemental Table 4). Interestingly, in vertebrates Brachyury forms a feedback loop with canonical Wnt3 signaling as well as with FGF signaling in the tailbud and the neuromesodermal progenitors (NMPs) (Garriock et al., 2015;Goto et al., 2017;Martin and Kimelman, 2010). In Nematostella and Strongylocentrotus, Wnt signaling is active in the blastopore region where brachyury is expressed (Davidson et al., 2002;Gross and McClay, 2001b;Kusserow et al., 2005;Logan et al., 1999). Moreover, functional manipulation of the Wnt signaling pathway has demonstrated that brachyury is activated by Wnt/beta-catenin signaling (Lebedeva et al., 2021;Röttinger et al., 2012). To functionally investigate the regulation of key target blastoporal transcription factors and signaling pathway members by Brachyury, we carried out WMISH in wild-type and morphant embryos in Nematostella and Strongylocentrotus. We focused on Brachyury ChIP targets that were regulated in our RNA-seq experiments or were previously shown to have an overlapping or mutually exclusive expression pattern with brachyury. In both sea anemone and sea urchin, brachyury is upregulated upon Morpholino mediated knockdown, suggesting a negative feedback loop of Brachyury on its own expression (Figs. 4A,5A).
In Strongylocentrotus, zygotic brachyury expression initiates at the blastopore just before the onset of gastrulation (18 hpf). Shortly before the archenteron breaks through (24 hpf), the future stomodaeum (oral ectoderm) of the late gastrula also starts to express brachyury (Fig. 1C). Therefore, Brachyury is expected to have different target genes in these two domains. For instance, different genes are coexpressed with brachyury in these two domains (e.g., hox11/13b is co-expressed with brachyury in the blastoporal region, while gsc is co-expressed with brachyury in the oral ectoderm ( Fig. 5A and Supplemental Fig. 5B). We validated several ChIP target genes that showed no change of expression in the differential gene expression analyses in brachyury morphants. Among the brachyury co-expressed genes of the presumptive endoderm, foxA (Tu et al., 2006) and wnt16 (Cui et al., 2014) are downregulated, while hox11/13b (Arenas-Mena et al., 2006) and otx (Peter and Davidson, 2010) are upregulated in Brachyury knockdowns. Moreover, the expression domain of the mesodermal markers ets1 and ese (Rizzo et al., 2006) is expanded, suggesting a repressive function of Brachyury on these genes in wild type embryos ( Fig. 5A and Supplemental Fig. 5B). These results show that, like in vertebrates, Nematostella and Strongylocentrotus Brachyury can act both as direct activator and repressor on different target genes. By contrast, some target genes, which are also co-expressed with brachyury in the blastoporal region (e.g., soxC, wnt8, eve) (McClay et al., 2018;Peter and Davidson, 2010;Wikramanayake et al., 2004) remain unaffected by the brachyury knockdown (Supplemental Fig. 5A), suggesting that other factors than brachyury may play a decisive role in their regulation.
Taken together, in both the sea urchin and the sea anemone Brachyury activates numerous members of the Wnt-beta catenin and the Wnt-PCP pathway in the blastoporal region of the gastrulating embryo, as well as foxA, which in sea urchin is expressed both in the blastopore and the future stomodeum. Since brachyury in turn is downstream of canonical Wnt signaling both in Strongylocentrotus and in Nematostella (as well as other cnidarians), we conclude that in cnidarians, sea urchins and vertebrates, there is a positive feedback loop of Brachyury and Wnt signaling at the blastopore and its derivative tissues.

Nematostella Brachyury represses neuronal genes in the oral domain
One of the surprising recent ndings in vertebrates was that Brachyury has a dual role in the differentiation of "neuromesodermal progenitors": it activates mesodermal genes but directly also represses neuronal genes, thereby antagonizing the function of Sox2 in the promotion of neural fate (Koch et al., 2017). Surprisingly, in Nematostella, Brachyury also regulates many genes that are expressed at the aboral half, thus not overlapping with brachyury expression. During gastrulation, the aboral half is the domain of early neurogenesis (Nakanishi et al., 2012). Notably, many of these aborally expressed target genes are involved in neurogenesis, e.g. the achaete-scute homolog ash-A, islet-1, tbx2/3, masterblind (mbnl), rfx4, noc and lhx-1, which are expressed in single scattered cells, typical for early neurogenic factors (Layden and Martindale, 2014;Layden et al., 2012;Nakanishi et al., 2012;Rentzsch et al., 2017). Interestingly, in Brachyury morphants, the expression domain of many of these genes expands to the oral domain, suggesting that Brachyury might be directly inhibiting early neurogenesis in the oral domain of the blastopore (Fig. 4B, Supplemental Fig. 4). Of note, the putative neuronal markers that show no signi cant change of expression do also not show a single cell pattern, but rather a global endodermal expression pattern, suggesting that the function of Brachyury on the expression of neuronal genes is restricted to the ectoderm (compare Fig. 4 and Supplemental Fig. 4).

Sea urchin Brachyury is involved in proper ectodermal patterning and fate establishment of the different ectodermal domains
We also investigated the effect of Strongylocentrotus Brachyury on genes of the oral ectodermal territory, where brachyury is expressed during gastrulation. The oral ectoderm domain is adjacent to the anterior neuroectoderm (ANE) of the embryo that will later give rise to the apical organ of the pluteus larva. The expression of the oral ectoderm (future stomodaeum) gene goosecoid (Angerer et al., 2001), which shows an overlapping expression with brachyury in this region in wild-type embryos, is severely reduced in Brachyury morphants (Supplemental Fig. 5B, 5D). Similarly, the anterior neuroectoderm marker Nkx2.1 is downregulated in Brachyury morphants as revealed both by differential RNA sequencing and immunohistochemistry (Supplemental Fig. 5C, D). Notably, several other genes that are also expressed in the ANE (e.g. fgf9/16/20, fzd5/8, six3/6, soxB2) (Anishchenko et al., 2018;Range, 2018;Wei et al., 2009) are downregulated in Brachyury morphants (Fig. 5B, Supplemental Fig. 5B, C, D), suggesting that Brachyury might have an indirect role in patterning the anterior neural ectoderm. Moreover, genes involved in the speci cation of the oral-aboral (ventral-dorsal) axis such as nodal and bmp2/4 (Duboc et al., 2010) that partially co-localize with brachyury seem to be severely affected, with the expression of nodal reduced and the expression of bmp2/4 abolished (Fig. 5B). Fgf9/16/20 that is expressed in the oral (Röttinger et al., 2008) is also abolished (Fig. 5B). This suggests that the proper formation of the oralaboral axis is compromised in Brachyury morphants. Notably, the expression of the aboral ectoderm marker spec1 (Otim et al., 2004) is also reduced in Brachyury morphants, con rming the general aberrant patterning of the ectoderm.
To summarize, similar to vertebrates, in Nematostella Brachyury appears to contribute to the inhibition of neuronal differentiation in the blastoporal domain, thereby restricting it to the aboral part. In Strongylocentrotus neuronal markers (e.g., hbn, nkx2.1, nkx3.2; Supplemental Fig. 5C and D) appear to be downregulated upon Bra knockdown, however, mostly indirectly due to the fact that in the absence of Brachyury the anterior neuroectoderm is generally mis-patterned. This ectodermal mispatterning results not only in loss of distinct ectodermal territories (oral-aboral-ANE), but also in the inability of the ANE to promote neurogenesis and, thus, proper neuronal differentiation.

Phylogenetic comparisons reveal ancestral targets of Brachyury
In order to reconstruct the evolutionary changes of Brachyury function that led to its mesoderm determination role, conserved among chordates, we aimed to compare the direct Brachyury target genes in organisms with available information, ranging from protists to vertebrates. To this end, we rst generated a reliable set of orthologs in the species under consideration. We used OMA (Orthologous MAtrix) for its ability to report strict orthologs by verifying pairs of genes (Altenhoff et al., 2013(Altenhoff et al., , 2018Roth et al., 2008). All the translated peptide sequences were ltered against proteins shorter than 100 amino acids. This ltering resulted in 8381 (C. owczarzaki), 25729 (N. vectensis), 28658 (S. purpuratus), 25425 (C. intestinalis), 39662 (X. tropicalis) and 21317 (M. musculus) proteins. Our orthology analysis resulted in 1882 orthologous groups shared by all species (4791 orthologous groups where at least 5 species have orthologs, 7922 orthologous groups where at least 4 species have orthologs, 13149 orthologous groups where at least 3 species have orthologs, and 28918 orthologous groups where at least 2 species have orthologs).
We then compared cnidarian and sea urchin Brachyury targets identi ed by ChIP-seq with the published target gene sets from the urochordate Ciona intestinalis (Katikala et al., 2013;Kubo et al., 2010) and two vertebrates, the mouse Mus musculus (Koch et al., 2017) and the frog Xenopus tropicalis (Gentsch et al., 2013) . We added the putative Brachyury targets from the protist Capsaspora owczarzaki determined by the search of Brachyury binding motifs in ATAC-seq peaks (Sebé-Pedrós et al., 2016) to this comparison. The OMA matrix allowed us to assign pair-wise homologs between these species and thereby reconstruct the ancestral target genes at each of the phylogenetic branching nodes (Figure 6; Supplemental Table 1).
Since the mouse has the best-annotated genome, we used it as a reference for the annotations in the analyses. As in all investigated metazoans, regulation of transcription, multicellular organism development and signal transduction are among the most dominant GO categories of the Brachyury targets, we rst focused our analysis on transcription factors. Notably, brachyury is the only target gene that is shared among all organisms, suggesting that self-regulation is an ancestral feature and under strong selection (Supplemental Fig. 6). However, while in some organisms Brachyury negatively regulates itself (Nematostella, Strongylocentrotus), in other species (mouse) Brachyury positively regulates itself.
There is ample evidence that transcription factor binding sites and hence the corresponding target genes can diverge rapidly (Dermitzakis and Clark, 2002). However, important target genes might be under higher selection pressure and maintained over longer evolutionary distances. Thus, by comparison of lineages of various phylogenetic distances, we should be able to identify ancestral target genes of distinct phylogenetic groups. Therefore, in order to reconstruct how the Brachyury target gene set changed during the course of animal evolution, we sought to determine by pairwise comparisons, which target genes evolved early and have been maintained in several lineages and which genes were recruited only in more recent lineages. A gene was considered an ancestral target for a given ancestor node when a homolog is shared between the earliest branching organism (outgroup) and at least one of the other species (ingroup). The pairwise analysis suggested that at least four transcription factors (TFs) were targets in the last common ancestor of Capsaspora and Metazoa (brachyury, rnf113a1, rfx1/2/3, zic) ( Fig. 6;  Supplemental Fig. 7). In fact, brachyury is the only gene that is a target gene of Brachyury in all investigated organisms, suggesting that autoregulation is an ancestral and highly constrained feature of the role of Brachyury. At node II (Nematostella + deuterostomes), we found 58 homologous target TFs shared between Nematostella and at least one of the deuterostome species (Strongylocentrotus, Ciona, Xenopus, Mus) ( Figure 6). To understand the functional consequences of species-speci c Brachyury target genes, we marked the target TFs by their described and annotated function from Xenopus or mouse using the UniProt description (Bateman, 2019). Using this approach, we nd 32% (19/58) of the ancestral node II genes (sea anemone + deuterostomes) having a neuronal function (e.g., ash, atonal, dmbx, emx2, otp, otx, irx, six3/6, soxB1, foxB, foxP, islet, mbnl, mycN), compared to only 14% (8/58) that have a mesodermal function in vertebrates. Notably, some of these, for example tbx2/3 and hes1, are likely to have a role in neuronal differentiation in Nematostella (Fig. 4B). At node III (deuterostomes) we found 18 TFs shared between Strongylocentrotus and the chordates, of which only one (GATA4/5/6) has a de ned mesodermal function in vertebrates, while 26% (5/19) ancestral target TFs have an assigned neuronal function. When Ciona was compared with the vertebrates (node IV), we detected 19 TFs as shared Brachyury targets. These shared targets are therefore supposedly restricted to the role of Brachyury in notochord development. However, the number of shared target genes between C. intestinalis and vertebrates is relatively modest and many of them are not exclusively, if at all, expressed in the notochord (Reeves et al., 2017). Although Ciona shares with vertebrates the expression of Brachyury in the notochord, only three (Six1/2, Tbx18, Tbx6) have described functions in mesoderm formation in vertebrates. One interesting target gene is tbx6, which has crucial functions downstream and in conjunction with Brachyury in the formation of paraxial mesoderm in vertebrates (Chapman and Papaioannou, 1998;Chapman et al., 1996Chapman et al., , 2003Gri n et al., 1998). There is no tbx6 homolog in Cnidaria, nor in any other non-bilaterian phylum. By contrast, there are putative homologs in the deuterostome Strongylocentrotus and representatives of the protostomes (e.g., the mollusc Lottia gigantea; (Supplemental Fig.7A), however, our data suggest that in the sea urchin tbx6 does not seem to be a target of Brachyury. We conclude that tbx6 likely evolved by a gene duplication event in the bilaterian ancestor (although our phylogenetic analysis (Supplemental Fig. 7A) did not provide strong support for a monophyletic group of vertebrate and invertebrate Tbx6). In Ciona, there are three Ci-tbx6 paralogous genes, which show an expression in the developing paraxial muscles, i.e., complementary to brachyury (Takatori et al., 2004). This suggests that tbx6 is likely negatively regulated by Brachyury in C. intestinalis. Thus, tbx6 was recruited as a novel target only in the chordates.
Finally, at node V we identi ed 25 TFs target of Brachyury that are shared between Xenopus and mouse. Since homologs of these genes are not found as targets in the other organisms, they supposedly have evolved newly or have been recruited in vertebrates. Several of these genes are known as crucial regulators involved in mesoderm development. This includes mesp1/2, mesogenin, twist, mef2c and smad6 (Fig. 6). Thus, these TFs are likely to be involved in conveying the role of Brachyury as a mesoderm determination factor during early gastrulation and a pioneering factor in neuromesodermal progenitor differentiation.
By comparison, only few (4/27) target TFs with a neuronal function are not shared with the other lineages.
In summary, this phylogenetic analysis of the target genes shows that numerous ancestral target genes that predate the split of cnidarians and bilaterians have a neuronal function at least in vertebrates (often in most investigated species), whereas a large fraction of key mesodermal target genes was only acquired at the level of the vertebrates.

Conservation of binding motifs of Brachyury across metazoans
Our genome-wide ChIP-seq experiments in the sea anemone Nematostella and the sea urchin Strongylocentrotus con rmed that the binding motif of Brachyury is deeply conserved throughout the animal kingdom and likely even extending to the protist C. owczarzaki. This is consistent with the extreme conservation of the T-domain of Brachyury over >700 Million years and the induction of mesoderm formation by the ectopic expression of Capsaspora or cnidarian brachyury mRNA in the animal hemisphere of the frog embryo (Bielen et al., 2007a;Marcellini et al., 2003;Sebé-Pedrós et al., 2016). Notably, in all organisms we nd both single binding sites as well as palindromes, which were rst predicted by in vitro selex studies Müller and Herrmann, 1997) . Palindromic motifs bind dimers of Brachyury, and the binding of dimers is expected to be stronger than the binding of monomers to single motifs. This may have a functional consequence on the effect on gene expression.
Indeed, we detected statistically signi cantly higher peaks with palindrome motifs in Nematostella. However, we could not detect any signi cant correlation of the distribution of single and palindromic binding sites with respect to the distance to the TSS nor to the category of target genes, nor to the suppression or activation of the target gene (Supplemental Fig. 3A, data not shown). A large fraction of the detected Brachyury binding sites coincide with promoters and enhancers de ned by speci c combinations of chromatin modi cations in Nematostella and by ATAC-seq de ned open chromatin sites in Strongylocentrotus, suggesting that Brachyury binds primarily to predicted cisregulatory elements. However, it is noteworthy that a certain fraction of the Brachyury peaks is also found in regions of closed chromatin. This raises the possibility that in Nematostella as well as in Strongylocentrotus Brachyury can act as a pioneer transcription factor as it has recently been described in vertebrates (Beisaw et al., 2018;Tosic et al., 2019).
Interestingly, a large fraction (about 23% in Nematostella and 27% in Strongylocentrotus) of the Brachyury binding peaks do not harbor a Brachyury motif (passing signi cance cutoffs) (Supplemental Fig. 3). While some of these could be false positives, we assume that many are events of indirect DNA binding through another transcription factor. This suggests that Brachyury might act in concert with a variety of other transcription factors. Indeed, we do nd motifs of Sox, Fox and homeodomain TFs enriched in the vicinity of the Brachyury motifs, suggesting either competitive or cooperative binding.
Interestingly, in Nematostella brachyury is co-expressed with soxB1, soxB2a (also termed sox1), foxA, and foxB and notably, all of them are also direct targets of Brachyury, suggesting that they are interconnected in a gene regulatory network (see below). Similarly, in Strongylocentrotus, homeodomain genes such as hox11/13b, otx, six3/6, are not only co-expressed, but also direct targets of Brachyury, suggesting conserved network interconnectedness. Future work should focus on investigating the physical interactions of these transcription factors.
A conserved feedback loop of Brachyury, FoxA and Wnt signaling Gene Ontology analyses suggest that in Nematostella as well as in Strongylocentrotus, numerous target genes are related to the control of development and transcription (i.e., transcription factors, members of signaling pathways), suggesting that Brachyury plays an ancestral upstream role in the regulation of developmental processes. In particular, we nd numerous members of the Wnt, FGF, Notch, Hedgehog and BMP pathways as targets of Brachyury in both investigated species (Supplemental Table 4). Several Wnt pathway members like dkk1, wnt5b , wnt8a, wnt9a, wnt11b (PCP) are also targets of Brachyury in vertebrates (Gentsch et al., 2013;Koch et al., 2017), and brachyury is also regulated by Wnt3 and by eFGF (Arnold et al., 2000) in a feedback loop. It was recently shown that Wnt signaling acts upstream of brachyury in Nematostella (Lebedeva et al., 2021;Pukhlyakova et al., 2018;Röttinger et al., 2012). Several TCF binding motifs in the brachyury promoter suggest that this regulation is likely direct (data not shown). In this study, our ChIP-seq data in Nematostella have further revealed that numerous wnt ligands (wnt1,2,3,4,6,11,16,A), all four frizzled receptors (frz1, 4, 5/8, 10), dishevelled, and b-catenin are direct targets of Brachyury, demonstrating the direct regulation of the canonical as well as the putative PCP pathway by Brachyury and the existence of a feedback loop like in vertebrates. Similarly, in Strongylocentrotus, brachyury is activated by the canonical Wnt pathway and, in turn, Brachyury targets several genes of the Wnt pathway (wnt1, wnt9, wnt7, wnt16, fz1/2, fz5/8, fz10), highlighting a deeply conserved feedback loop of Wnt and Brachyury (Croce et al., 2011;Cui et al., 2014;Sethi et al., 2012).
Since, in most bilaterian species, brachyury expression and Wnt signaling is co-localized at the blastopore lip, we postulate that brachyury forms an ancestral gene regulatory feedback loop at the blastopore. Furthermore, because Wnt signaling is known to play a crucial role in establishing the primary body axis in vertebrates, echinoderms (and hemichordates) and cnidarians (Niehrs, 2010), we propose that one of the ancestral roles of Brachyury in metazoans was in axial patterning. This kernel is also conserved in some insects (Lengyel and Iwaki, 2002). The blastopore lip appears to act as an ancestral axial organizer where the Wnt-Brachyury feedback loop is conveying the axis-inducing property (Kraus et al., 2007(Kraus et al., , 2016. In line with this, knockdown of brachyury abolishes the axis induction capacity of the blastopore lip in Nematostella (Lebedeva et al., 2021). While we do not have evidence for a control of brachyury by FGF signaling at present, we also nd FGF ligands (fgf8 in Nematostella and fgf9/16/20 in Strongylocentrotus), as well as foxA (and foxB in Nematostella) as deeply conserved target genes of Brachyury, suggesting that these genes might also belong to the axis patterning kernel (Fig. 7A).

The evolutionary conservation of a blastoporal expression of brachyury
Outside metazoans, brachyury and several other T-box genes are found in the protist Capsaspora owzcarzaki, as well as in certain fungi, indicating that it evolved before the emergence of the metazoans (Sebé-Pedrós and Ruiz-Trillo, 2017;Sebe-Pedros et al., 2011. The function of brachyury in these organisms is not known but the T-domain of Capsaspora is su ciently conserved to induce mesoderm when expressed in the frog embryo (Sebe-Pedros et al., 2013).
Among non-bilaterian animals, brachyury has also been detected in sponges (Leininger et al., 2014;Sebe-Pedros et al., 2013), ctenophores and placozoans Spring, 2003, 2005;Yamada et al., 2007Yamada et al., , 2010. While the correspondence of the body plans of sponges and placozoans to cnidarian or bilaterian body plans appears less clear at rst glance, there is evidence for a conservation of the main body axes (Adamska et al., 2007;Dubuc et al., 2019). Furthermore, functional studies in the ctenophore Mnemiopsis leydii suggests that Brachyury is expressed at the blastopore and is required for the invagination and formation of the stomodeum (Yamada et al., 2010).
The blastoporal expression is also conserved in a broad range of bilaterian species, for instance the polychaete Platynereis dumerilii (Arendt et al., 2001), the mollusc Patella vulgata (Lartillot et al., 2002), the hemichordate Ptychodera ava (Peterson et al., 1999;Tagawa et al., 1998) and the ecdysozoan phylum Priapulida (Martín-Durán et al., 2012). However, in different species of brachiopods (Spiralia, Lophophorata), potentially divergent functions of Brachyury are observed with respect to the formation of the gut openings (Martín-Durán et al., 2017). In Phoronopsis harmeri, a member of the brachiopod sister group Phoronida (Luo et al., 2018), brachyury seems to be involved in the patterning of the ventral endoderm and ventral ectoderm, as well as the formation of the posterior ciliary band (Andrikou et al., 2019). In Drosophila melanogaster, the brachyury ortholog brachyenteron is required for the ectodermal hindgut formation, where it is co-expressed with the foxA ortholog forkhead (Berns et al., 2008;Kusch and Reuter, 1999;Reim et al., 2017;Singer et al., 1996). Similar roles have been reported in short germ insects, e.g. Tribolium castaneum (Berns et al., 2008) and the intermediate germ insect, the cricket Gryllus bimaculatus (Shinmyo et al., 2006). Thus, despite their divergent mechanisms of posterior elongation, these insects have retained a conserved role for brachyury in hindgut development and posterior visceral muscles. Taken together, these comparative data strongly support an ancestral blastoporal expression of brachyury (Swalla, 2006;Technau, 2001) , and where investigated, often in conjunction with foxA expression and Wnt signaling. Functional studies in some of these species support a role of brachyury in axial patterning and invagination as well as fore/hindgut development.

The evolution of the mesoderm inducing role of Brachyury in chordates
Since Brachyury has a conserved function in mesoderm development in vertebrates and is required for the notochord development in urochordates, we looked for conserved downstream mesodermal genes in Nematostella and Strongylocentrotus. Our analysis showed that only few genes with a role in mesoderm formation in vertebrates are also Brachyury target genes in Nematostella (e.g., tbx2/3, msx, gata1/2/3; see Supplemental Table 1, Fig. 6). However, some of these genes appear to have also a neuronal function in Nematostella (e.g., tbx2/3, gata1/2/3), suggesting that they acquired a role in mesoderm formation in the deuterostome lineage. Instead, a surprisingly large number of mostly negatively regulated TFs that are involved in neurogenesis are shared between Nematostella and vertebrates, suggesting that the negative regulation of neuronal regulators is more ancestral. Surprisingly, in Strongylocentrotus, the few mesodermal Brachyury target genes (e.g., ets1, gcm, foxn2/3, see Supplemental Table 1 Some of the key factors driving mesoderm differentiation in chordates, such as tbx6 and myoD, appear to have evolved only in the deuterostome or bilaterian lineage, respectively, since they do not exist in cnidarian genomes. In Strongylocentrotus, two myoD paralogous genes are present (myoD1 and myoD2), both with distinct mesodermal functions in skeletogenesis and myogenesis (Andrikou et al., 2013), but there is no evidence that they are target genes of Brachyury. Quite the contrary, our functional studies suggest that, in the sea urchin, brachyury represses mesodermal fate (through repression of the mesoderm speci cation gene ese) and rather promotes endoderm formation.
Tbx6 is a key downstream gene in vertebrates, conveying the mesoderm determination function of Brachyury in feed-forward loops. It is a crucial factor in the development of the paraxial mesoderm, and mutation of tbx6 leads to the formation of two supernumerary neural tubes instead of somites (Chapman and Papaioannou, 1998;Chapman et al., 2003). Tbx6 is also a Brachyury target gene in the ascidian C. intestinalis (Kubo et al., 2010), expressed in the future paraxial muscle and therefore likely to be negatively regulated by Brachyury. The phylogenetic analysis (Supplemental Fig. 7A) shows that tbx6 genes arose in deuterostomes, although the monophyly of invertebrate and vertebrate tbx6 is not supported, hence the vertebrate tbx6 genes might also have an independent origin. In any case, the evolution of chordate tbx6 (arisen through duplication of a member of the tbx2/3/4/5 families) was probably a crucial step in the recruitment of Brachyury to a mesodermal function.
Besides tbx6, another crucial evolutionary step was the acquisition of other key mesodermal downstream differentiation determinants in vertebrates, such as vegt, mesp1/2, twist, mesogenin, mef2, myf5, myf6, myod as direct target genes of Brachyury, linking early mesoderm formation to muscle differentiation. All these genes have been shown to play decisive roles in the speci cation of the mesoderm in vertebrates.
While the Myogenic Regulatory Factors (MRFs) (myf5, myf6, mrf4, myod) are not present in the cnidarian genomes, mesp and twist exist, but are neither direct targets of Brachyury in Nematostella or Strongylocentrotus, nor do they seem to play a major role in axis formation or germ layer formation (Genikhovich and Technau, 2011;Martindale et al., 2004). Thus, the evolution of several key mesodermal determinants as targets of Brachyury includes both adoption of conserved genes to mesoderm formation as well as the evolution of novel genes in the vertebrate lineage.

Brachyury and the evolution of neuromesodermal progenitors
Recent studies have shown a dual function of Brachyury in vertebrates during axis elongation at the tailbud stage. During this stage, the tailbud contains an undifferentiated cell population, which expresses both brachyury as well as sox2, a marker for pluripotency and early neural fate (Bergsland et al., 2011;Wilson et al., 2009;Wymeersch et al., 2016). This cell population has been termed neuromesodermal progenitors (NMPs) and even though its precise de nition and characterization is still debated, there is evidence that these cells have the potential to differentiate into neural tube or into mesoderm, depending on the concentration of Wnt3a. However, some studies argue that brachyury in combination with WNT and FGF signaling is su cient to drive cells to NMPs state (Mugele et al., 2018). Several studies support the view that Sox2 expression drives cells towards neural fate, and a high level of brachyury expression promotes mesodermal fate (Bouldin et al., 2015;Takemoto et al., 2011). Functional as well as ChIP-seq studies have shown that Brachyury and Sox2 mutually inhibit each other, which eventually leads to a segregation of neuronal and mesodermal cell fates (Koch et al., 2017). Brachyury also represses a number of other neuronal target genes, in line with its role in suppressing neuronal fate (Gentsch et al., 2013). Sox2 belongs to the subfamily of SoxB1, together with its vertebrate-speci c paralogs Sox1 and Sox3. There is no one-to-one ortholog of sox2 in cnidarians and sea urchin, but there is a cnidarian and sea urchin soxB1 homolog. Interestingly, while there is no evidence for soxB1 being a target in Strongylocentrotus (although both soxB1 and soxB2 are differentially expressed in knockdowns, see Supplemental Fig. 5B-C), at least in Nematostella, soxB1 is a direct target of Brachyury (Fig. 3D). SoxB1 is maternally expressed throughout the early embryo of Nematostella, but becomes cleared from the oral pole at the same time when brachyury starts to be expressed (Lebedeva et al., 2021), suggesting a possible suppression of soxB1 by Brachyury. During normal development, the expression of soxB1 retracts further to the ectodermal aboral pole, where neurogenesis starts. After gastrulation, a second expression domain in the pharynx appears at the oral side (Magie et al., 2005). Notably, in brachyury knockdown, the retraction of soxB1 from the oral side occurs normally, but the pharynx expression is missing (Lebedeva et al., 2021). This suggests that the pharynx expression of soxB1 (which overlaps with brachyury) is dependent on Brachyury. Thus, like in vertebrate NMPs, a sox2 homolog soxB1 is a direct target in Nematostella, but to date there is no evidence for mutual inhibition. Furthermore, although soxB1 is expressed in the domain of aboral neurogenesis, there is no evidence for a direct role of soxB1 in promoting neural development in Nematostella (Supplemental Fig. 4C).

Concluding Remarks
Our genome-wide target gene screening in an early branching metazoan, the sea anemone Nematostella and the early branching deuterostome, the sea urchin Strongylocentrotus, revealed an ancestral kernel comprised of a feedback loop of Wnt signaling with Brachyury/ FoxA and likely also FGF signaling (Fig.   7A). As this ancestral loop is implicated in the early speci cation of the primary body axis and the gastrulation site, we therefore conclude that the ancestral function of Brachyury was as a crucial component of axis formation and patterning (Fig. 7B). The comparison with target genes in vertebrates, where Brachyury plays a crucial and early role as a mesodermal determinant and a pioneering factor (Fig.   7B), reveals that the sea anemone Brachyury represses numerous neuronal genes, many of which are also repressed as in vertebrates, but Nematostella Brachyury regulates very few homologs of vertebrate mesodermal genes (Fig. 7B). Even in the deuterostome sea urchin, only few mesodermal targets are shared with vertebrates and in general, mesoderm is rather repressed, while endoderm formation is promoted (Fig. 7B). We hypothesize that the novel acquisition of some key factors promoting mesodermal fate (e.g., tbx6, msgn, mesp1/2, myoD, twist), in a feed-forward loop with Tbx6, have evolved only in the chordate lineage, leading to a dramatic and important shift in function of this ancient transcription factor. Figure 1 Lineages sampled in this study and Brachyury protein localization in relation to brachyury mRNA expression. (A) The tree shows the lineages sampled in this study for the comparative analysis of brachyury function highlighted in red. Genome scale Brachyury binding data from Capsaspora owczarzaki, Nematostella vectensis, Strongylocentrotus purpuratus, Ciona intestinalis, Xenopus tropicalis, and Mus musculus, representing Filasterea, Cnidaria, Echinodermata, Tunicata, Vertebrata   Validation of Brachyury target genes in Nematostella by WMISH after morpholino induced knockdown of Brachyury. Relative spatial expression of (A) oral markers brachyury, wnt1, wnt3, foxA, foxB, fgf8 and (B) aboral markers wnt2, ash1, lhx, tbx2/3, rfx4, ngn). The asterisks indicate the oral side of the embryo. Quantitative data are found in Supplemental table 3. Scale bar corresponds to 50 µm. regulatory input in the endoderm (yellow), ectoderm (red) and ectoderm (blue). WMISH of wnt16, foxA, otx, hox11/13b, ets1 (magenta) relative to Brachyury expression (green). White arrowheads indicate coexpression. (B) WMISH of fgf9/16/20, nodal, bmp2/4, frizzled, six3, spec1 in Brachyury morphants compared to wild type 24h embryos. Pictures are projections of confocal stacks. Nuclei are labeled blue with DAPI. # indicates that the ChIP target is also present in the differential RNAseq dataset. Arrows show the embryonic domain in which we see an effect (red arrow: mesodermally derived; yellow: endodermally derived; blue ectodermally derived). l/v: lateral view; v/v vegetal view; o/v oral view; a/v aboral view. up: upregulated gene; down; down-regulated gene. The scale bar is 20 µm.

Figure 6
Apomorphic and synapomorphic genes target genes of Brachyury coding for transcription factors. To assess the evolutionary origin of vertebrate Brachyury target genes and their function, apomorphic and synapomorphic target genes were assigned to nodes of common ancestors. Transcription factor coding genes targeted by Brachyury that are shared between an outgroup and at least one member of the ingroup are considered ancestral for each node. The function of the target genes was annotated according to their described role in mouse in UniProt (see Supplemental table 2). Genes with a function in mesoderm are highlighted in orange, genes with a neural function are highlighted in blue. Genes with both neural and mesodermal functions are highlighted in purple. Genes with other functions are not highlighted. Vertebrate-speci c paralogs were collapsed to metazoan homologs at nodes I-IV. Note that the common ancestor of cnidarians and bilaterians shared numerous transcription factor target genes involved in neural development.