An ancestral Wnt–Brachyury feedback loop in axial patterning and recruitment of mesoderm-determining target genes

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 chromatin immunoprecipitation sequencing 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 signalling 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 endodermal and 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. Brachyury is an early mesoderm determinant and neural repressor in vertebrates. Comparative Brachyury target screens between a sea anemone and a sea urchin reveal an ancestral gene regulatory feedback loop involved in axial patterning, with conserved endodermal and neuronal, but not mesodermal, targets.

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 chromatin immunoprecipitation sequencing 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 signalling 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 endodermal and 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.
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 [1][2][3][4][5] . It is expressed in the presumptive mesoderm and continues to be expressed after gastrulation in the notochord, a defining structure of all chordates 1,6-8 .
Heterozygous mouse mutants of the T gene locus show a shortened tail (hence the Greek name Brachy, meaning short, and ury, meaning tail), while homozygous mutants die in utero, lacking most of the body posterior to the forelimbs 9,10 . Knockdown of brachyury in Xenopus also impairs the posterior development of the embryo, while overexpression of brachyury messenger RNA in explanted animal caps leads to mesoderm differentiation 1,2,4,11 . Similar roles have been found in other Article https://doi.org/10.1038/s41559-022-01905-w 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 chromatin immunoprecipitation sequencing (ChIP-seq) and identified 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 signalling. 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 find a substantial number of neuronal target genes shared between Nematostella, Strongylocentrotus and vertebrates, many of which are repressed as revealed by brachyury knockdowns. This suggests that axial patterning, endoderm specification 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.

Genome-wide detection of Brachyury binding sites
To reconstruct the evolution of Brachyury function, we analysed its genomic targets in representatives of phylogenetically informative groups (Fig. 1a), the sea anemone Nematostella and the sea urchin Strongylocentrotus. We generated antibodies against Nematostella and Strongylocentrotus Brachyury proteins 31 and confirmed their specificity (Fig. 1b, 31 . We next used these antibodies to perform ChIP-seq in early gastrula-stage embryos. For Nematostella, we identified 2,389 putative binding sites in two highly reproducible replicates, which mapped to 1,543 putative target genes (Fig. 2a,b; Methods; Extended Data Fig. 1e). For Strongylocentrotus, we identified 490 binding sites, which corresponded to 391 target genes ( Fig. 2b and Extended Data Fig. 1f). The lower number of detected binding sites in the sea urchin might be due to less effective ChIP enrichments or reflect a biological difference. Notably, a large fraction of the Brachyury peaks in Nematostella (850 out of 2,389) overlapped with the previously 2,559 identified enhancers at the gastrula stage 45 as defined by the combination of p300, H3K27ac and H3K4me1, suggesting that about one-third of the identified enhancers are bound by Brachyury at the gastrula stage (Extended Data Fig. 2a,b). Similarly, a large fraction (338 out of 490) of the Brachyury binding sites detected by ChIP-seq in Strongylocentrotus are found in open chromatin as detected by ATAC-seq (Extended Data Fig. 2c). However, we also found a significant fraction in both species that are not associated with active enhancers or open chromatin (Discussion). To compare our findings with other species, we re-analysed previously published ChIP-seq data 21,46 using our peak calling thresholds and found for mouse and Xenopus 4,000 and 2,497 peaks, respectively. Corroborating the published results 21,46 , our called mouse peaks corresponded to 3,060 putative target genes, while in Xenopus 1,376 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 reflects the genome sizes of the species. In species with larger genome size the peaks are distributed over a larger range, sometimes beyond 100 kilobases (kb), while in species with smaller genome size the peaks are much more in the vicinity of the transcription start site (TSS) (Fig. 2h-i). This is particularly obvious in the case of Nematostella, where ~50% of vertebrates 6,12 , demonstrating the conserved role for Brachyury as a pivotal mesoderm regulator across vertebrates.
Brachyury exhibits a dynamic expression pattern in vertebrate embryos. At the gastrula stage it is expressed around the blastopore and during anterior-posterior axis elongation, at the tailbud stage, brachyury expression is mainly restricted to the most caudal part of the embryo and the notochord 1,13 . This caudal part of the developing embryo contains a population of bipotent cells called neuromesodermal progenitor cells (NMPs) 14 , which exhibit temporal differences in gene expression over developmental time 15,16 and contribute to the spinal cord [17][18][19] paraxial mesoderm 17,18 and notochord development 17,20 .
Some studies define these NMPs by the expression of brachyury and sox2, a member of the soxB1 gene family 15,18,21 , while others argue that brachyury expression along with Wnt/Fibroblast Growth Factor (FGF) signalling is sufficient for the cells to become NMPs even in the absence of Sox2 protein 22 . Regardless of the definition of NMPs, the mesoderm-defining role of Brachyury has been established within these cell populations. The notochord-specific 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 23 . In amphioxus, there are two brachyury paralogues. 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 probably plays a similar role in this cephalochordate 24 , while bra1 appears to be more restricted to the notochord 25 . 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.
In numerous investigated invertebrate species, brachyury is expressed around the blastopore, similar to the early brachyury expression in vertebrates 26 . For instance, in Echinodermata (for example, sea urchins and sea stars), members of the non-chordate deuterostome group Ambulacraria, brachyury is expressed at the vegetal pole marking the margin of the blastopore [27][28][29] . Later, brachyury expression is confined to the most vegetal part of the hindgut surrounding the anus (proctodeum) and the stomodeum 30,31 . Sea urchin Brachyury has been functionally implicated in the regulation of gastrulation movements 32 . In protostomes, the role of Brachyury shows remarkable similarities to deuterostomes during gastrulation, for instance in the annelid 33,34 and mollusc embryos 35,36 .
Brachyury has also been studied in non-bilaterian animals, including cnidarians. For instance, in the freshwater polyp Hydra, the brachyury homologue Hybra1 is expressed at the hypostome and acts as an early marker of head formation during budding, regeneration and embryonic development [37][38][39] . 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 40 . Functional studies in N. vectensis and the coral Acropora digitifera suggest that brachyury may be involved in the formation of the pharynx and the endoderm-ectoderm boundary 41,42 . Strikingly, when Hybra1 mRNA is injected into blastomeres of the animal pole of Xenopus embryos, it is capable of inducing mesoderm formation with high efficiency, showing that in the context of the vertebrate, the cnidarian transcription factor behaves like its cognate vertebrate homologue 43,44 . 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, for this transcription factor to acquire such diverse developmental roles? Therefore, to shed light on the all Brachyury binding sites are located within 1 kb upstream or downstream of the TSS, compared to all three deuterostome species (5-12%) ( Fig. 2d-f). The smaller genome size also makes the gene density in Nematostella higher (about one gene per 10 kb) than that of other species concerned 47,48 . Different transcription factors tend to have strikingly different positional specificity within the enhancer regions 49 . From the binding profile of Brachyury ( Fig. 2d-g), Brachyury seems to be strongly concentrated around the TSS 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 50,51 , which was supported by ChIP-seq in Xenopus embryos 46 . We detected 700 and 238 palindromic Bra binding motifs in the sea anemone and sea urchin Brachyury ChIPseq peaks, respectively. This indicates that the capacity for dimerization is an ancestral feature of Brachyury (Fig. 2b,c). We detected significantly higher Brachyury peaks with palindromes compared to peaks with half-palindromes (Extended Data Fig. 4a), suggesting a higher binding affinity of the dimer compared to the half-palindromes. However, we found no significant correlation between the distribution of half-palindromes or palindromes with the expression level or gene ontology (GO) category of the target genes. Thus, the biological meaning of the distribution of half-palindromes 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 52 . This suggests that the DNA-binding domain and the corresponding binding motif have been conserved from protists to humans, indicating 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 43,44,53 .
To detect potential cofactors 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 ~50 base pairs (bp) from the peak centre, suggesting possible cobinding (Extended Data Fig. 4b). When scanning the peak sequences with the software tool Find Individual Motif Occurrences (FIMO) for known motifs matching to the enriched motifs identified in any species, we found a similar distribution of putative binding sites for homoeodomain, bHLH, Pax, HMG (Sox) and Fox proteins in all species (Extended Data Fig. 4b). Especially in Nematostella, a substantial number of putative Sox, Fox and Homeobox motif sites are found together with Brachyury (Extended Data Fig. 4b). This is of particular interest as the expression of foxA and foxB, as well as soxB1 and soxB2, is partially overlapping with brachyury in Nematostella (see below).

Brachyury can act as activator or repressor of target genes
Next, we wished to study the role of Nematostella and Strongylocentrotus Brachyury on gene expression by knockdown of brachyury function (Extended Data Fig. 3). We generated three to six replicate transcriptomes of gastrula-stage embryos each injected with control, splicing or translation-blocking morpholinos, respectively (Supplementary Information). In both species, in response to brachyury knockdown, we detected numerous upregulated or downregulated genes, a small fraction of which overlapped with the ChIP-seq target gene set (Supplementary Information; Extended Data Fig. 3) 21,31,[54][55][56][57] . GO analyses of the Brachyury ChIP targets found in Nematostella and Strongylocentrotus shows an enrichment in categories of transcriptional regulation, Wnt signalling, multicellular organism development, signal transduction and biological regulation (Extended Data Fig. 4c; Supplementary Information), which is similar to the GO categories of Brachyury ChIP targets reported in vertebrates 21,46,58 . This suggests that, in these species, Brachyury is a developmental regulator that mainly acts positively or negatively to control other developmental genes.

Brachyury and Wnt form an ancestral regulatory loop
When analysing the target genes in Nematostella and in Strongylocentrotus, we found that many members of Wnt, FGF, Notch and Bone Morphogenetic Protein (BMP) signalling pathways as well as several transcription factors were among the direct target genes (Figs. 3 and 5; Supplementary Information). Interestingly, at the tailbud stage of vertebrates Brachyury forms a feedback loop with canonical Wnt signalling as well as with FGF signalling in the NMPs 54,59,60 . In Nematostella and Strongylocentrotus, Wnt signalling is active in the blastopore region where brachyury is expressed 32,61-63 . Moreover, functional manipulation of the Wnt signalling pathway has demonstrated that brachyury is activated by Wnt/beta-catenin signalling [64][65][66] . To functionally investigate the regulation of key target blastoporal transcription factors and signalling pathway members by Brachyury, we carried out whole mount in situ hybridization (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 expression is upregulated upon morpholino-mediated knockdown, suggesting a negative feedback loop of Brachyury on its own gene expression (Fig. 4a,c).
In Nematostella, other genes with an overlapping expression with brachyury are either abolished in the overlapping domains (for example, wnt1, wnt3, foxA, foxB, fgf8a, ephrinB2 and myc2) or strongly downregulated (wnt4 and wntA) upon brachyury knockdown ( Fig. 4a and Extended Data Fig. 5). By contrast, wnt2, which shows a complementary expression pattern to brachyury, expands into the blastopore lip domain in brachyury morphants (Fig. 4a), indicating that Brachyury represses wnt2 from invading the oral territory in wild-type embryos. To correctly identify the target gene of the binding site, we make use of the orthologue information. Of the two genes closest to the binding site, we select the gene that is also a Brachyury target in other species (Extended Data Fig. 2d). b, Summary table showing the number of ChIP-seq peaks, their target genes and the presence of centrally enriched (Centrimo/ MEME-ChIP) palindrome or half-palindrome motifs within the peak region. (N/A not analysed) c, Brachyury motifs (palindromic and/or half-palindromic) in ChIP-seq peaks. In each comparison the upper motif is a vertebrate motif from the JASPAR database, while the lower motif is the one found in our respective datasets. For mouse, the enriched motifs are the same as the one found in the JASPAR database. d-g, Brachyury binding profile in mouse (d), Xenopus (e), Strongylocentrotus (f) and Nematostella (g) around the TSS with 95% confidence interval (light grey shaded area) calculated with 1,000 bootstrap. x axis, distance to TSS; y axis, density of peaks overlapping each basepair (as in distance to/from TSS) on the x axis. The plot was generated using the plotAvgProf function from the ChipSeeker R package. h, Genomic distribution of Brachyury peaks with respect to the TSS over the whole genome. All genomic regions around TSSs were binned into six bins depending on their distance to the TSS and split into upstream (left side of the plot) and downstream (right side of the plot) of the TSS. The percentage of peaks overlapping each bin is plotted on the x axis. This plot was generated using the plotAnnoBar function of the ChipSeeker R package. i, Distribution of the Brachyury binding sites with respect to genomic features such as promoters (defined here as 1,000 bp upstream and 200 bp downstream of the TSS), introns, exons, untranslated regions (UTRs) and downstream (defined here as 300 bp downstream of a gene) and intergenic regions. The percentage of Brachyury peaks (x axis) overlapping each region as shown in the diagram above is shown for each species (y axis).  In Strongylocentrotus, zygotic brachyury expression initiates at the blastopore just before the onset of gastrulation (18 hours post fertilization (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 co-expressed with brachyury in these two domains (for example, hox11/13b is co-expressed with brachyury in the blastoporal region, while gsc is co-expressed with brachyury in the oral ectoderm ( Fig. 4c and Extended Data Fig. 5b)). We validated several ChIP target genes that showed a change of expression in the differential gene expression analyses in brachyury morphants. Among the brachyury co-expressed genes of the presumptive endoderm, foxA 67 and wnt16 68 are downregulated, while hox11/13b 69 and otx 70 are upregulated in Brachyury knockdowns. Moreover, the mesodermal markers ets1 and ese 71 are ectopically expressed in the presumptive endoderm, suggesting a repressive function of Brachyury on these genes in this domain of wild-type embryos ( Fig. 4c and Extended Data Fig. 6b). These results show that, as 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 show overlapping expression with brachyury in the blastoporal region (for example, soxC, wnt8 and eve) 70,72,73 remain unaffected by the brachyury knockdown (Extended Data 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-β-catenin and the Wnt-Planar Cell Polarity (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 signalling 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 signalling at the blastopore and its derivative tissues.

Nematostella Brachyury represses neuronal genes in the oral domain
One of the surprising recent findings 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 21 . 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 74 . Notably, many of these aborally expressed target genes are involved in neurogenesis, for example the achaete-scute homologue ash-A, islet-1, tbx2/3, masterblind, rfx4-like, noc and lhx-1, which are expressed in single scattered cells, typical for early neurogenic factors [74][75][76][77] . 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 and Extended Data Fig. 5). Of note, the putative neuronal markers that show no significant change of expression also do not show a single-cell pattern but rather a global gastrodermal expression pattern, suggesting that the function of Brachyury on the expression of neuronal genes is restricted to the ectoderm (compare Fig. 4b and Extended Data Fig. 5).

Sea urchin Brachyury is involved in ectodermal patterning
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 78 , which shows an overlapping expression with brachyury in this region in wild-type embryos, is severely reduced in Brachyury morphants (Extended Data Fig. 6b,d). Similarly, the ANE marker Nkx2.1 is downregulated in Brachyury morphants as revealed both by differential RNA sequencing and immunohistochemistry (Extended Data Fig. 6c,d). Notably, several other genes that are also expressed in the ANE (for example, fgf9/16/20, fzd5/8, six3/6 and soxB2) [79][80][81] are downregulated in Brachyury morphants ( Fig. 4d and Extended Data Fig. 6b-d), suggesting that Brachyury might have an indirect role in patterning the anterior neural ectoderm. Moreover, genes involved in the specification of the oral-aboral (ventral-dorsal) axis such as nodal and bmp2/4 (ref. 82 ) that partially colocalize with brachyury seem to be severely affected, with the expression of nodal reduced and the expression of bmp2/4 abolished (Fig. 4d). Fgf9/16/20 that is expressed in the oral side 83 is also abolished (Fig. 4d). This suggests that the proper formation of the oral-aboral axis is compromised in brachyury morphants. Notably, the expression of the aboral ectoderm marker spec1 (ref. 84 ) is also reduced in brachyury morphants, confirming the general aberrant patterning of the ectoderm. Such a role of sea urchin Brachyury at the oral-aboral axis organizer is in line with what was recently suggested by single-cell transcriptomics data 85 .
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 (for example, hbn, nkx2.1 and nkx3.2; Extended Data Fig. 5c,d) appear to be downregulated upon Bra knockdown, however, mostly indirectly due to the fact that in the absence of Brachyury the ANE is generally mispatterned. 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 Brachyury targets
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 first generated a reliable set of orthologues in the species under consideration using OMA (Orthologous Matrix; Methods). We then compared cnidarian and sea urchin Brachyury targets identified by ChIP-seq with the published target gene sets from the urochordate C. intestinalis 86,87 and two vertebrates, the mouse Mus musculus 21 and the frog Xenopus tropicalis 46 . We added to this comparison the putative Brachyury targets from the protist Capsaspora owczarzaki determined by the search of Brachyury binding motifs in ATAC-seq peaks 52 . The OMA allowed us to assign pairwise homologues between these species and thereby reconstruct the ancestral target genes at each of the   Table 1). Since the mouse has the best-annotated genome, we used it as a reference for the functional 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 first focused our analysis on Article https://doi.org/10.1038/s41559-022-01905-w 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 (Extended Data Fig. 7). There is ample evidence that transcription factor binding sites and hence the corresponding target genes can diverge rapidly 55 . 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, to reconstruct how the Brachyury target gene set was conserved or was 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 homologue is shared between the earliest branching organism (outgroup) and at least one of the other ingroup species (Methods). To unravel shared ancestral and lineage-specific developmental regulators, we focused the analyses on target genes coding for transcription factors (TFs) (Fig. 5), while the whole set of target genes is found in Supplementary Table 1. The pairwise analysis suggested that six genes were TFs coding targets in the last common ancestor of Capsaspora and Metazoa (Fig. 5). Node II is of particular interest, as it reflects the common ancestor of the diploblastic sea anemone Nematostella with the deuterostomes. Notably, we found numerous homologous target TFs shared between Nematostella and at least one of the deuterostome species (Strongylocentrotus, Ciona, Xenopus and Mus) (Fig. 5).
To understand the functional consequences of shared or species-specific Brachyury target genes, we annotated the target genes encoding TFs by their differential expression in different germ layers and tissues during mouse gastrulation. As a reference, we used a recently published single-cell transcriptome dataset from mouse embryos corresponding to the gastrula stage (stage E8.0 and 8.5; Methods) 56 . When we calculated the log 2 normalized endodermal versus mesodermal as well as neuronal versus mesodermal gene expression levels, we found that ancestral node II genes (Nematostella plus deuterostomes) tend to have significantly higher expression levels in endodermal and neuronal cell types compared to mesodermal cell types (Fig. 5). A similar trend of more endodermal targets compared to mesodermal genes was detected at node III (deuterostomes) and node IV (chordates), albeit less pronounced than at node II. Only at node V (genes exclusively shared among the two vertebrates), we detected more mesodermal genes than neuronal or endodermal genes, in line with the conserved function of Brachyury in mesoderm specification in vertebrates (Fig. 5). However, in this single-cell expression dataset brachyury and noto-two notochord marker genes-are annotated as definitive endoderm, although notochord is commonly regarded as dorsal mesoderm (but see Discussion), which might influence the bias of endodermal over mesodermal targets at node II. We therefore removed the definitive endoderm in an alternative analysis (Extended Data Fig. 9). Yet, endodermal targets are still enriched over mesodermal targets. Moreover, another, independent recent single-cell dataset from the mouse (stage E8.5; Extended Data Fig. 9) 57 also showed an enrichment of neuronal and endodermal over mesodermal targets at node II, in contrast to node V (Extended Data Fig. 9) 88 . Thus, among the shared target genes coding for transcription factors, genes with a role in endodermal or neuronal cell types are more ancestral than the mesodermal target genes, indicating a shift of target genes of Brachyury in the chordate and vertebrate lineage.
Next, we were interested as to whether there are some crucial TFs that would make the difference in the evolutionary change of Brachyury function or whether it is rather a stepwise and gradual shift. We found 16 TFs shared between Strongylocentrotus and the chordates, of which only few have a defined mesodermal function in vertebrates. When Ciona was compared with the vertebrates (node IV), we detected 21 TFs as shared Brachyury targets. Given that Brachyury is expressed in the notochord in all chordates, 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, in line with previous studies 89 , although this might be again influenced by the annotation of notochordal cells as definitive endoderm 88 . Although Ciona shares with vertebrates the expression of Brachyury in the notochord, only five (six1/2, tbx18, tbx6, cdx1 and lmx1) have assigned 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 [90][91][92][93] . There is no tbx6 homologue in Cnidaria, nor in any other non-bilaterian phylum. By contrast, there are putative homologues in the deuterostome Strongylocentrotus and representatives of the protostomes (for example, the mollusc Lottia gigantea; Extended Data Fig. 8a), however, our data suggest that in the sea urchin tbx6 does not seem to be a target of Brachyury. We conclude that tbx6 probably evolved by a gene duplication event in the bilaterian ancestor, although our phylogenetic analysis (Extended Data Fig. 8a) 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, that is, complementary to brachyury 94 . This suggests that tbx6 is probably negatively regulated by Brachyury in C. intestinalis. Thus, tbx6 was recruited as a new target only in the chordates.
Node V target genes (32 of which are TFs) are shared between Xenopus and mouse. Since homologues of these genes are not found as targets in the other organisms, they supposedly have evolved newly or have been recruited in vertebrates. Unlike node II genes, we find node V genes to be more highly expressed in mesodermal cell types ( Fig. 5 and Extended Data Fig. 9). In addition, several of these genes are known as crucial regulators involved in mesoderm development. Among others, this includes mesp1/2, mesogenin, twist, mef2c and smad6. 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.
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, whereas a large fraction of key mesodermal target genes was only acquired at the level of the vertebrates. 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 the respective node. The expression of the target genes in neuronal versus mesodermal cell types was annotated using single-cell RNAseq data from E8.0 and E8.5 mouse embryos. a, Heatmap of the log 2 -fold change (logFC) of neuronal (turquoise) or endodermal (yellow) versus mesodermal (purple) gene expression for each gene using the corresponding mouse gene symbols as they appeared in the single-cell dataset in each node (Methods). b, Boxplot of the log 2 -fold change (y axis) of endodermal versus mesodermal gene expression for all genes per node (x axis). Note that node II is enriched in endodermally expressed genes while node V is enriched in mesodermally expressed genes (P = 0.008, Wilcoxon rank sum test). c, Boxplot of the log 2 -fold change (y axis) of neuronal versus mesodermal gene expression for all genes per node (x axis). Note that node II is enriched in neuronal expression while node V is enriched in mesodermal expression (P = 0.008, Wilcoxon rank sum test). The boxes range from the 25th to the 75th percentile and the horizontal lines represent the median. Outliers are shown as dots.

Discussion
In this study, we revealed the genome-wide targets of Brachyury in two invertebrates-the diploblast sea anemone N. vectensis and a basally branching deuterostome, the sea urchin S. purpuratus-and compared them with known sets of target genes from similar ChIP experiments in other organisms.

Limitations of this study
We took particular care to use comparable stages of Nematostella and sea urchin embryos for our ChIP-seq experiments at early gastrulation, however, as with any such comparative analyses, we cannot fully rule out that both missing or detected target genes are due to specific differences in the developmental stage, tissue origin and the experimental design. For example, all ChIP-seq data were generated using gastrulation stage embryos, except for the mouse dataset 21 , which was generated using in vitro differentiated NMP cells from embryonic stem cells and which might lack notochordal genes. ChIP-seq studies focusing on the mouse notochord might therefore add more ancestral genes. Also, each ChIP-seq experiment has a specific sensitivity and therefore leads to more or less detected target genes, which is probably the main reason why we detect fewer target genes in the sea urchin. The cross-species analyses of ancestral target genes (Fig. 5) also relies on the annotation of cell clusters of the mouse single-cell RNA-seq data. For instance, others have annotated early notochordal cells as definitive endoderm 88 , which might be considered a false annotation. However, we would like to note that in the protochordate Amphioxus the notochord derives by delamination from the dorsal archenteron midline, independent of the paraxial mesoderm, indicating a close link of notochord and endoderm. In summary, as a cautionary note, we cannot make statements and draw conclusions on individual target genes that might be missing in a given dataset. On the other hand, the mouse and frog datasets share a large number of target genes, many of which have also been validated experimentally at other stages. This suggests that a large fraction of the conserved target genes can be captured.

A conserved blastoporal Brachyury and Wnt feedback loop
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 (Fig. 6c) 53,95,96 . The function of brachyury in these organisms is not known but the T-domain of Capsaspora is sufficiently conserved to induce mesoderm when expressed in the frog embryo 53 . Our analyses of the binding sites confirm that the binding motif of Brachyury is deeply conserved throughout the animal kingdom and probably even extending to the protist C. owczarzaki, hence over >700 million years. As in vertebrates, we do find both palindromic sites indicating binding of dimers, as well as half-palindromes (Extended Data Fig. 3; Supplementary Information). Co-occurring motif analyses suggest either competitive or cooperative binding of other transcription factors with a bias for Sox, Fox and homoeodomain factors (Extended Data Fig. 3; Supplementary Information).
One of the striking features of Brachyury is the highly conserved expression domain around the blastopore in almost all animals studied, with only few exceptions (Supplementary Information). This raises the question of what might stabilize the expression domain over hundreds of millions of years in so many diverging phyla.
In vertebrates, several Wnt pathway members like dkk1, wnt5b, wnt8a, wnt9a and wnt11b (PCP) are targets of Brachyury 21,46 and brachyury is also regulated by Wnt3 and by eFGF 97 in a feedback loop. It was recently shown that Wnt signalling acts upstream of brachyury in Nematostella 64,65,98 . In this study, our ChIP-seq data in Nematostella have further revealed that numerous wnt ligands, all four frizzled receptors, dishevelled and β-catenin are direct targets of Brachyury, demonstrating the direct regulation of the canonical as well as the putative PCP pathway by Brachyury. Similarly, in Strongylocentrotus, brachyury is activated by the canonical Wnt pathway and, in turn, Brachyury targets several genes of the Wnt pathway, highlighting a deeply conserved feedback loop of Wnt and Brachyury (Fig. 6a) 27,68,99 . Since brachyury expression and Wnt signalling are colocalized at the blastopore lip in most bilaterian species, we postulate that brachyury forms an ancestral gene regulatory feedback loop at the blastopore. Furthermore, because Wnt signalling is known to play a crucial role in establishing the primary body axis in vertebrates, ambulacrarians (echinoderms and hemichordates) and cnidarians 100 , we propose that one of the ancestral roles of Brachyury in metazoans was in axial patterning (Fig. 6b). In support of this idea, knockdown of brachyury abolishes the axis induction capacity of the blastopore lip in Nematostella 64 . While we do not have evidence for a control of brachyury by FGF signalling at present, we also find 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. 6a,b). Animal icons are either self-generated (e.g. Nematostella) or from Phylopic (http://phylopic.org/): Echinodermata (Creative Commons Attribution 3.0 Unported) was created by H. N. Eyster; mammals (Creative Commons Attribution 3.0 Unported) by S. Miranda-Rottmann.

Evolution of the mesoderm inducing role of Bra in chordates
Since Brachyury has a conserved function in mesoderm development in vertebrates and is required for notochord development in urochordates (and possibly in cephalochordates), 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 (for example, tbx2/3, tbx20, mbnl1 and msx, gata1/2/3; Fig. 5, Supplementary Table 1 and Extended Data Fig. 6). 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. Likewise, a considerable number of genes that are endodermally expressed in the mouse are targets in Nematostella and hence were supposedly target genes in the common ancestor of cnidarians and bilaterians (node II). Surprisingly, in Strongylocentrotus, the few mesodermal Brachyury target genes (for example, ets1, gcm and foxn2/3; Supplementary Table 1 and Fig. 4c,d) are repressed by Brachyury in the presumptive endoderm, in stark contrast to what is observed in chordates. This difference highlights the plasticity of Brachyury in acquiring new functions and targeting new genes.
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 101 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 specification 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 90,92 . Tbx6 is also a Brachyury target gene in the ascidian C. intestinalis 87 , expressed in the future paraxial muscle, hence mutually exclusive to the brachyury expressing notochord and therefore likely to be negatively regulated by Brachyury. 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.
Article https://doi.org/10.1038/s41559-022-01905-w 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 and 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 specification of the mesoderm in vertebrates. While the myogenic regulatory factors (MRFs) (myf5, myf6, mrf4 and myoD) are not present in the cnidarian genomes, twist exists but is neither a direct target of Brachyury in Nematostella nor in Strongylocentrotus, nor does it seem to play a major role in axis formation or germ layer formation 102,103 . 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 new genes in the vertebrate lineage.

Brachyury and the evolution of neuromesodermal progenitors
The tailbud stage of vertebrates has been viewed as secondary axis formation 104 or even as a continuation of the gastrulation process itself 105 . The tailbud harbours a population of NMPs, which has the potential to differentiate into neural tube or into mesoderm, depending on the concentration of Wnt3a. NMPs express both brachyury and sox2, a marker for pluripotency and early neural fate 15,106,107 . Several studies support the view that sox2 expression drives the cells towards neural fate and a high level of brachyury expression promotes mesodermal fate 21,108,109 . ChIP-seq as well as functional studies have shown that Brachyury and Sox2 mutually inhibit each other, which eventually leads to a segregation of neuronal and mesodermal cell fates 21 . Brachyury also represses a number of other neuronal target genes, in line with its role in suppressing neuronal fate 46 .
Sox2 belongs to the subfamily of SoxB1, together with its vertebrate-specific paralogues Sox1 and Sox3. There is no one-to-one orthologue of sox2 in cnidarians and sea urchin but there is a cnidarian and sea urchin soxB1 homologue 110 . Interestingly, while there is currently no evidence that soxB1 is a target in Strongylocentrotus (although both soxB1 and soxB2 are differentially expressed in Bra knockdowns; Extended Data Fig. 5b,c), at least in Nematostella, soxB1 is a direct target of Brachyury (Fig. 3d). But unlike in vertebrate NMPs, there to date is no evidence for mutual inhibition of SoxB1 and Brachyury in Nematostella. Furthermore, although soxB1 is expressed in the domain of aboral neurogenesis, functional perturbation of soxB1 does not argue for a direct role of soxB1 in promoting neural development in Nematostella (Extended Data Fig. 4c). However, Brachyury contributes to the repression of neuronal differentiation during embryogenesis. This effect might be in concert with the role of Wnt signalling. Since a significant number of these neuronal genes are also negatively regulated by Brachyury in the NMPs of vertebrates, we conclude that repression of neuronal differentiation is an (unexpected) ancestral role of Brachyury, regardless of its axial position. This is not to say that NMPs have their origin in the common ancestor of cnidarians and bilaterians but that part of the NMP gene regulatory network and the role of Brachyury therein has a deeper origin.

Concluding remarks
Our comparative approach of genome-wide target genes of a key developmental regulator, Brachyury, has revealed both deeply conserved kernels as well as key acquisitions of new target genes, that were either recruited or evolved newly by gene duplications. Brachyury has a surprisingly conserved expression pattern around the blastopore in most animal phyla. Since this blastoporal expression is part of the ancestral kernel, we postulate that such feedback loops can evolutionarily stabilize a given expression pattern and its ancestral role. Nevertheless, the function of the gene can still drastically evolve, by the acquisition of only a few, but crucial, new target genes, which in the case of Brachyury, convey the new function in feed-forward loops. At present, the number of comparative studies of genome-wide target-gene screens of conserved transcription factors is very limited but future approaches over a wide range of distantly related organisms may offer insights into the general mechanisms of robustness and evolvability of gene regulatory networks in the evolution of animal body plans.

Animal culture
Nematostella polyps were kept as previously described 111,112 at 18 °C and largely in the dark. Spawning was induced by a combination of light and elevated temperature 111 . 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 filtered Mediterranean seawater diluted 9:1 with deionized water 81,101,113,114 .

Antibody generation
A 6×His-tagged full-length clone of Nematostella brachyury was expressed in Escherichia coli strain BL21. After transformation, bacteria were grown overnight in LB supermedium + ampicillin. Then, 10 ml of LB medium/ampicillin were inoculated with 500 µl of overnight culture, grown for 3 h (optical density 600 = 0.5-0.6), isopropyl-β-d-thiogalactoside was added to a final concentration of 1 mM and bacteria culture was incubated for 21 h at 4 °C. Brachyury protein was purified 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 two different rabbits were raised against the recombinantly produced C-terminal domain excluding the T-box domain. The antibodies were affinity purified by PRIMM Biotech before their use.
For both Nematostella and Strongylocentrotus the reactivity of polyclonal antibodies from each rabbit were assessed by both Western blot and immunohistochemistry, confirming their specificity (Extended Data Fig. 1). For the ChIP replicates, antibodies from the two different rabbits were used for each.

Chromatin immunoprecipitation and library preparation
ChIP-seq was carried out essentially as previously published 45 . In brief, embryos at gastrula stage were fixed with 2% formaldehyde for 12 min. Nuclei were collected by Dounce homogenizing the embryos. Chromatin was fragmented to an average size of 100-200 bp using a Covaris S1 ultrasonicator with the following settings: duty cycle, 20%; intensity, 5; cycles/burst, 200; duration, 240 s; mode, frequency sweeping. To carry out the immunoprecipitation, 300 µg of protein (as quantified by Bradford assay) of the fragmented chromatin was incubated with 5 µg of Nematostella or Strongylocentrotus anti-Brachyury. Chromatin from ~20 immunoprecipitations was pooled to obtain ~2 ng 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 (catalogue no. IP-102-1001). The quantity and quality (size distribution) of the libraries was confirmed using an Agilent Bioanalyzer. Deep sequencing was performed at the Vienna BioCenter Core Facilities (https://www.viennabiocenter.org/facilities/) with 50 bp of SE HiSeq 2500 reads. We obtained 45

Orthologue detection
For orthologue detection we used the OMA database [115][116][117] . Before we subjected our in silico translated protein sequences to OMA, Article https://doi.org/10.1038/s41559-022-01905-w we excluded peptide sequences smaller than 100 amino acids. To detect the orthologues 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 21 (GEO: GSE93524) and Xenopus 46 (GEO: GSE48560). Reads from M. musculus, X. tropicalis, S. purpuratus and N. vectensis were mapped to genome versions GRCm38 (mm10), v.9.1, v.5.0 and v.1.0, respectively, with BWA using the BWA-MEM algorithm 118 . ChIP-seq peaks were called using Peakzilla 119 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 filter 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 <2 for mouse and Xenopus, <0.5 for Strongylocentrotus and <5 for Nematostella. Of note, the number of detected target genes in Strongylocentrotus is noticeably lower than in the other species, which could be either due to biological reasons (fewer target genes) or technical reasons (less sensitive antibody) or a combination of both. This may lead to an underestimation of the number of conserved targets in our comparative analyses.

Target gene selection
To identify target genes, a commonly used procedure is to select the closest gene to an identified peak. However, in more compact genomes like Nematostella or Strongylocentrotus (five to ten genes per 100 kb) peaks in the intergenic space are often only marginally less distant to a gene in the other direction. Therefore, we identified the two closest genes to the binding site. Next, using orthologue information obtained earlier with OMA 115,116 , we determined whether either of the two closest genes has an orthologue 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 orthologue 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 (Extended Data Fig. 2d). In case an equal number of the orthologous target genes was found, we kept both putative targets and marked them with ++ (Supplementary Table 1).

Motif analysis
To find and compare the binding motifs of Brachyury we used the MEME-ChIP (4.11.2) suite 120 . For the number of motifs to be predicted in a sequence we selected 'anr ' (any number of repetitions) mode and 1 × 10 −10 for the e-value cutoff criterion. The discovered motifs were scanned against the non-redundant sequence profiles obtained from the JASPAR database 88 . The other parameters were set to default (a detailed list of parameters can be found on https://github.com/ technau/brachyury_grn).

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 finds all the motifs in a peak and their absolute distances from the centre of the peak. It also groups the motifs into motif groups on the basis of motif families in JASPAR 88 . In case of overlapping motifs, it selects the motif with better False Discovery Rate (FDR)-adjusted P value assigned by FIMO from the MEME-ChIP suite 120 .

Brachyury knockdowns and RNA-seq
To knockdown Brachyury function in Nematostella, we used two non-overlapping antisense morpholinos (MO) against brachyury: the translation-blocking morpholino BraATG-MO TCGTCCGAGTG-CATGTCCGACTATG and the splice-blocking morpholino BraSplic-eMO TCCCTGGTTGTCAACCATACCGTCC. A total 3-5 pl of both morpholinos were injected at the concentration of 500 µM. A total 50 µg ml −1 of Dextran-Alexa Fluor 488 MW 10000 (ThermoFisher) was co-injected as tracer to ensure proper delivery of MO 121 . Standard morpholino (Genetools) StdMO CCTCTTACCTCAGTTACAATTTATA was injected as a control at the same concentration. PolyA-enriched RNA samples were collected and processed for library production using the Lexogen kit and sequencing (50 bp single-end HiSeqV4, ~30 million reads per sample (for more details see https://github.com/ technau/brachyury_grn).
To knockdown Strongylocentrotus brachyury, we injected 2-4 pl of a 200 µM solution of morpholino previously characterized with the sequence: CGCTCATTGCAGGCATAGTGGCG 31 . We quantified the mortality at 24 hpf and discarded all experiments with <80% survival. Embryos were either fixed for WMISH or used for isolation of RNA extraction at 24 h. PolyA-enriched RNA samples were collected and processed for library production using the Illumina Truseq RNA sample preparation kit (unstranded) and sequencing (100 bp paired-end HiSeq, ~50 million reads per sample).

Differential gene expression analyses
For Nematostella, we generated six biological replicates of each translational and splice morpholino. Principal component analysis showed that the replicates cluster together (Extended Data Fig. 3a). For Strongylocentrotus, we generated three biological replicates (Supplementary Fig. 3b). This was compared to three replicates of morpholino knockdown in Xenopus 46 (GSE48663) (Extended Data Fig. 3c) and two replicates of mouse knockout mutants 21 (GSE93524) (Extended Data Fig. 3d).
RNA-seq reads for Mus (mouse), Xenopus, Strongylocentrotus and Nematostella were subjected to initial quality control with FastQC (v.0.11.5) 122 and a quality control report was generated with MultiQC (v.1.6) 123 . The index for Mus, Xenopus and Nematostella species was generated for STAR aligner (v.2.5.3a) 124 and the reads were aligned to the respective genomes. Reads aligned to genomic features were counted with featureCounts (Subread v.1.6.2) 125,126 . The count files generated were used as input to SARTools (v.1.6.0) 127 , which is a wrapper R script around two widely used differential expression analysis tools, namely DESeq2 (ref. 128 ) and edgeR 129 . Genes with an FDR-corrected P value (q-value) of 0.05 are considered as differentially expressed. The RNA-seq reads for Strongylocentrotus were quasi-mapped to the S. purpuratus genome v.5 associated transcriptome using Salmon (v.0.11.3) 130 . The resulting count files were used in DESeq2 (ref. 128 ) for differential expression analysis, differentially expressed genes with P-adjusted value of <0.05 were considered significant. A quantitative analysis of shared orthologous direct or indirect targets that are differentially expressed in the respective organisms is found in Extended Data Fig. 3f.

Whole mount in situ hybridization
WMISH in Nematostella was carried out as described before 131 with a few minor modifications. Briefly, embryos were fixed in 4% paraformaldehyde (PFA) in phosphate-buffered saline (PBS) for 1 h at room temperature (RT) and then washed 5× with 0.2% Tween in PBS (PTw) at RT. The embryos were then washed once with 50% methanol in PTw and lastly transferred into 100% methanol and stored at −20 °C until further processing. Hybridization was performed at 1 ng riboprobe µl −1 , overnight at 63 °C and detection with alkaline phosphatase conjugated anti-digoxigenin (1:2,000) (anti-digoxigenin-AP, Fab fragments, Roche, REF 11093274910). If no strong signal and mostly background Article https://doi.org/10.1038/s41559-022-01905-w began to develop in the NBT/BCIP solution, the embryos were washed 5× in PTw and then left to stand in PBS with 2% w/v Triton X100 until the background was at acceptable levels. This was repeated until a strong signal was observed.
In Strongylocentrotus, WMISH was done as previously described 132 . Briefly, embryos were fixed 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 prehybridized for 3 h in the hybridization buffer. Hybridization was carried out with 0.1 ng µl −1 of probes and incubated for 1 week. The signal for colorimetric in situ hybridization was developed using anti-digoxigenin AP (alkaline phosphatase) conjugated antibody (anti-digoxigenin-AP, Fab fragments, Roche, REF 11093274910) and the AP substrate, while fluorescent in situ hybridization was developed with fluorophore-conjugated tyramide (1:400 reagent diluents, TSA Plus Cyanine 3/Cyanine 5, PerkinElmer, NEL752001KT).

Immunohistochemistry
Immunohistochemistry in Nematostella was essentially performed as described in ref. 131 with the following changes. Briefly, embryos were fixed at 4 °C in 4% PFA in PTwTxD (1× PBS with 0.2% Tween, 0.2% Triton X100 and 0.2% DMSO). All further washing steps were performed with PTwTxD. They were then washed once with PTwTxD and ice-cold acetone and put on ice. Once the embryos had settled down, the acetone was removed and they were washed 10× with the cold PTwTxD. The embryos were then blocked for 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 pre-adsorbed in the blocking 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 10× with cold PTwTxD and blocked once again for 2 h at RT in the blocking solution. Blocking solution was then replaced with pre-absorbed goat anti-rabbit (Alexa Flour 568, 1:2,000, Invitrogen) mixed with DAPI (final concentration 5 µg ml −1 ) and Alexa Fluor 488 phalloidin (final concentration 5 U ml −1 , Invitrogen) and incubated overnight at 4 °C. Next day, the embryos were washed 10× in PTwTxD, infiltrated overnight with Vectashield (Vector laboratories) and imaged.
Immunohistochemical detection of Brachyury in Strongylocentrotus was performed as described in ref. 132 . Briefly, the embryos were fixed in 4% PFA in filtered seawater for 15 min at RT, dehydrated in 100% ice-cold methanol for 1 min and washed several times with 1× PBS. PBS was replaced with a blocking solution (4% sheep serum, 1% BSA in 1× PBS) and the samples incubated for 1 h 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 (ref. 133 ), SoxB2 1:500 (ref. 134 )) 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 Alexa Fluor 488 (Invitrogen) diluted in PBS was added for one hour in PBS and then washed again several times.

Phylogenetic analysis
Phylogenetic analysis of T-Box, Sox, Zic, TFAP, NR2, RNF and Rfx family members was performed using maximum likelihood methods. Orthologues and paralogues 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 135 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 136 as implemented in IQ-TREE with 1,000 (T-Box and Sox) and 10,000 (Zic, TFAP, NR2, RNF, Rfx) bootstrap samples, respectively.

Target gene comparison between species
We used OMA for its ability to report strict orthologues by verifying pairs of genes 115,116,137 . All the translated peptide sequences were filtered against proteins shorter than 100 amino acids. This filtering resulted in 8,381 (C. owczarzaki), 25,729 (N. vectensis), 28,658 (S. purpuratus), 25,425 (C. intestinalis), 39,662 (X. tropicalis) and 21,317 (M. musculus) proteins. Our orthology analysis resulted in 1,882 orthologous groups shared by all species (4,791 orthologous groups where at least five species have orthologues, 7,922 orthologous groups where at least four species have orthologues, 13,149 orthologous groups where at least three species have orthologues and 28,918 orthologous groups where at least two species have orthologues).
We then selected the target genes coding for transcription factors using InterProScan 138 annotation on the basis of the presence of the DNA-binding domain from the Pfam database, then sorted target genes into five nodes on the basis of evolutionary lineages representing Metazoa, Bilateria, Deuterostomia, Chordata and Vertebrata. Genes were annotated according to the mouse gene annotations. Homologous target genes were determined by OMA. Afterwards the OMA annotation was manually improved for several genes by more detailed phylogenetic analyses (see sections above). Homologous target genes that are shared between at least two organisms were placed at the node where the two lineages diverged, regardless of how many other members of the ingroup also shared this target gene. For instance, a target gene detected in Nematostella and at least one of the bilaterian species would be placed at node II (common ancestor of cnidarians and bilaterians). Vertebrates often have multiple paralogues that have only one homologue in the earlier diverging species, which can lead to double assignments in the analysis. For instance, foxA, which is a conserved target in Nematostella, sea urchin and vertebrates, has duplicated into several paralogues in vertebrates 139 . OMA finds foxA2 as a shared target of Nematostella and mouse and foxA3 as a shared target of sea urchin and mouse. In such cases, we only kept the homologue at the lower node.
To compare the expression of ancestral target genes in different nodes in mouse gastrulation stage cells of neuronal, endodermal and mesodermal origin, we extracted the mouse ENSEMBL gene IDs for all genes in each node. In the case where we extended the OMA annotation by phylogenetic analysis and one ancestral gene was present in node II (for example, foxA) we used all mouse paralogues (for example, foxa1, foxa2 and foxa3). We then used the Mouse Gastrulation Data R package (v.1.8.0) to extract single-cell RNA-seq data for mouse embryos at stages E8.0 and E8.5 for all node target genes described above. We used the gene counts to calculate log normalized expression values (logNormCounts function of R package scuttle) to categorize genes as neuronal, endodermal or mesodermal (and others). Then we generated a pseudo-bulk expression matrix per mesodermal, endodermal and neuronal cell category for all genes and normalized the resulting counts to counts per million reads (cpm). Node genes were selected from this pseudo-bulk matrix and the log-fold-change of neuronal versus mesodermal expression or endodermal versus mesodermal expression was calculated as log 2 of neuronal (or endodermal) cpm + 1/ mesodermal cpm + 1.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The raw files from the ChIP-seq and RNA-seq experiments have been deposited to the NCBI GEO database under the accession Article https://doi.org/10.1038/s41559-022-01905-w numbers GSE182573 and GSE198320 for N. vectensis and S. purpuratus, respectively.

Code availability
All scripts developed for this study are available at GitHub at https://github.com/technau/brachyury_grn, https://github.com/ xxxmichixxx/MouseEmbryoSingleCell and https://github.com/ fmi-basel/gbuehler-MiniChip. and Xenopus is from morpholino induced knockdown of Bra transcripts while mouse data is a result of Bra knockout. (F) Overlap of direct (ChIP-seq detected) and indirect targets (DEGs) across different species. Each species was used as a query species (query dataset) and the genes determined that are differentially expressed, and also ChIP targets ('column 'direct targets in query') or not ChIP targets (column 'indirect targets in query'). The numbers refer to the number of orthologs in each one of the other species. The numbers are low, since the overlap between ChIP-seq targets and DEGs from RNA-seq experiments is only 1-10% within a given species. Fig. 4 | Motif analyses of Brachyury ChIP-peaks and Gene Ontology analyses of target genes. Motif analyses of Brachyury ChIP-peaks and Gene Ontology analyses of target genes. (A) Brachyury ChIP metaplots around Brachyury peaks containing palindrome or half-palindrome (single) or no Brachyury binding motifs. The average read count (normalized to a million reads) was calculated for Brachyury ChIP-seq reads for regions around peak summits spanning 2kb in 20bp windows. The shaded area around the lines represents the 95% confidence interval across peaks in a category. Note that in all four species, peaks with palindromes show significantly higher ChIP-seq read counts, which may serve as a proxy for the strength of Brachyury binding. (B) Intersection of motifs with the peak region. The Brachyury peak region (as identified using the Peakzilla algorithm) was scanned for presence of other transcription factor binding sites using Fimo (MEME-ChIP suite) with default settings. The resulting binding motifs were grouped by the transcription factor families Paired box (Pax), basic helix-loop-helix (bHLH), Forkhead box (Fox), homeobox (homeo), High mobility group (hmg), T-box (tbox). Presence of these motifs together with the Brachyury motif is highlighted in orange in the upset plots (iii, vi, ix, xii). Distance of these motifs from the peak centre was also tracked and is shown in adjoining plots (I, ii, iv, v, vii, viii, x, xi). (C) ChIP/DE gene set overlap and GO analysis. Overlap (dark grey) between ChIP targets (black) and differentially expressed (light grey) genes in M. musculus (i), X. tropicalis (iv), S. purpuratus (vi), N. vectensis (x). In C. owczarzaki (xiii) only the number of ATAC-seq peaks with Bra motifs is shown. Gene ontology analysis of Brachyury ChIP targets for M. musculus (ii), X. tropicalis (v), S. purpuratus (vii), N. vectensis (xi) and C. owczarzaki (xiv) and of differentially expressed genes after Bra KO/KD (iii, vi, ix, xii). Only gene ontology terms for biological process and molecular function were reported. The colour of the dot represents the score (−log(p-value)) assigned by topGO while the size of the dot represents the number of genes associated with the term. (D) GO analyses of target genes of peaks with or without a Bra consensus motif. Note that no significant difference can be detected.

Extended Data
Article https://doi.org/10.1038/s41559-022-01905-w Extended Data Fig. 9 | Expression of apomorphic and synapomorphic target genes of Brachyury in mouse E8.5 neuronal, endodermal, and mesodermal cell types. Expression of apomorphic and synapomorphic target genes of Brachyury in mouse E8.5 neuronal, endodermal, and mesodermal cell types. This is the same analysis as shown in Figure 6, except that definitive endoderm cells were removed from the analysis (A-C) or the single-cell gene expression dataset from Grosswendt et al. 121 (E, F) was used. A-C: The expression of the target genes in neuronal versus mesodermal cell types was annotated using single-cell RNA-seq data from E8.0 and E8.5 mouse embryos, as in Figure 6. However, since definitive endoderm cells in this dataset express notochord marker genes, we decided to remove these cells from the analysis before continuing as in Figure  6. (A) Heatmap of the log2 fold change (logFC) of neuronal (turquoise) or endodermal (yellow) vs mesodermal (purple) gene expression for each gene using the corresponding mouse gene symbols as they appeared in the singlecell dataset in each node (see Methods for details). B) Boxplot of the log2 fold change (y-axis) of endodermal vs mesodermal gene expression for all genes per node (x-axis). Note that node II is enriched in endodermal expression while node V is enriched in mesodermal expression (p-value = 0.003, Wilcoxon Rank Sum Test). (C) Boxplot of the log2 fold change (y-axis) of neuronal vs mesodermal gene expression for all genes per node (x-axis). Note that node II is enriched in neuronal expression while node V is enriched in mesodermal expression (p-value = 0.008, Wilcoxon Rank Sum Test). D-F: The expression of the target genes in neuronal versus mesodermal cell types was annotated using single-cell RNA-seq data from E8.5 mouse embryos (read counts per gene/cell from GEO accession GSE122187). To annotate cells as neuronal, endodermal or mesodermal, we used the information from Supplementary Table 2 of Grosswendt et al. For endodermal, we used Lineage = Eendo, for mesodermal we used Lineage = Emeso, and for neuronal we used the following cell states: 1,11,24, and 39, which, according to Supplementary Fig. 1i of Grosswendt et al corresponds to: neural ectoderm anterior, neural ectoderm posterior, fore/midbrain, and future spinal cord. (D) Heatmap of the log2 fold change (logFC) of neuronal (turquoise) or endodermal (yellow) vs mesodermal (purple) gene expression for each gene using the corresponding mouse gene symbols as they appeared in the single-cell dataset in each node (see Methods for details). (E) Boxplot of the log2 fold change (y-axis) of neuronal vs mesodermal gene expression for all genes per node (x-axis). Note that node II is enriched in neuronal expression while node V is enriched in mesodermal expression (p-value = 0.037, Wilcoxon Rank Sum Test). (F) Boxplot of the log2 fold change (y-axis) of endodermal vs mesodermal gene expression for all genes per node (x-axis). Note that when, as is the case in this dataset, 'gut' is the only annotated endodermal cell type, node II is not more enriched in endodermally expressed genes compared to node V (p-value = 0.23, Wilcoxon Rank Sum Test). The boxes range from the 25 th to the 75 th percentile and the horizontal lines represent the median. Outliers are shown as dots.
Corresponding author(s): Ulrich Techanu, Maria Ina Arnone Last updated by author(s): Jul 25, 2022 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

March 2021
Data Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy The raw files from the ChIP-seq and RNA-seq experiments have been deposited to the NCBI GEO database under the accession number GSE182573 and GSE198320 for N. vectensis and S. purpuratus respectively.

Human research participants
Policy information about studies involving human research participants and Sex and Gender in Research.

Reporting on sex and gender
Use the terms sex (biological attribute) and gender (shaped by social and cultural circumstances) carefully in order to avoid confusing both terms. Indicate if findings apply to only one sex or gender; describe whether sex and gender were considered in study design whether sex and/or gender was determined based on self-reporting or assigned and methods used. Provide in the source data disaggregated sex and gender data where this information has been collected, and consent has been obtained for sharing of individual-level data; provide overall numbers in this Reporting Summary. Please state if this information has not been collected. Report sex-and gender-based analyses where performed, justify reasons for lack of sex-and gender-based analysis.

Recruitment
Describe how participants were recruited. Outline any potential self-selection bias or other biases that may be present and how these are likely to impact results.

Ethics oversight
Identify the organization(s) that approved the study protocol.
Note that full information on the approval of the study protocol must also be provided in the manuscript.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences
Behavioural & social sciences Ecological, evolutionary & environmental sciences