Brachyury and SMAD signalling collaboratively orchestrate distinct mesoderm and endoderm gene regulatory networks in differentiating human embryonic stem cells

1The three primary germ layers (ectoderm, mesoderm and endoderm) arise from the pluripotent epiblast during gastrulation in the amniote embryo; this can be modelled in vitro using pluripotent stem cells. The brachyury gene (T, BRA) encodes a T-box transcription factor that plays an essential role in mesoderm formation. Brachyury is expressed first in the primitive streak during gastrulation, and later in the notochord and tailbud. Mice lacking one copy of the gene have a short tail, while homozygous embryos die around embryonic day (E) 9-10. The latter develop only the first seven somites and lack a proper notochord; in addition, they display defects in left-right asymmetry and in cell migration. Importantly, the mutant phenotype also includes severe defects in cardiovascular and placental development. Brachyury is widely used as the earliest marker of mesodermal and also endodermal differentiation in embryonic stem cell (ESC) studies and during gastrulation, because both these cell lineages derive from the primitive streak. Indeed, definitive endoderm progenitors co-express Bra, Sox17, Foxa2, Gsc and other endoderm markers. Based on morphological analyses, mouse Bra homozygotes seem to form a normal foregut but an abnormal hindgut. However, on a molecular level, the role of BRA in endoderm formation remains poorly understood.

2Genomic targets of BRA orthologues have previously been identified in zebrafish and Xenopus embryos, and in mouse embryoid bodies using ChIP-chip or ChIP-seq approaches. Recently, BRA target genes were identified in differentiating hESCs, but the regulatory impact of BRA was not functionally characterized in these experiments. Thus, despite its important role in vertebrate development and its widespread use in stem cell biology, the precise regulatory role of BRA and its genome-wide functional role are poorly understood in humans.

3Human embryonic stem cells (hESCs) represent the best available system in which to study the molecular mechanisms underpinning human embryonic development. We have developed two in vitro protocols (activin or BMP4 based) that can induce hESCs to differentiate into distinct cell populations, both expressing BRA. Strikingly, these populations have the characteristics of the anterior and posterior regions of the early primitive streak from which, respectively, endoderm or mesoderm cells arise in vivo.

4Here, we perform chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) to identify genomic targets of BRA in these two cell populations. Interestingly, these genome-wide binding events differ in activin-treated or BMP4-treated cells, suggesting that BRA interacts with the genome depending on the signalling environment and cell identity. We provide functional validation of these genomic binding events by modulating levels of BRA within different signalling environments in differentiating hESCs and by analysing the expression pattern of BRA target genes in mouse embryos that lack brachyury.

5Importantly, our experiments also indicate that BRA physically interacts with downstream effectors of activin or BMP4 signalling: SMAD2/3 in endoderm progenitors and SMAD1 in mesoderm progenitors. We conclude that BMP4-SMAD1 signalling and BRA action are essential for proper mesodermal differentiation, while simultaneously repressing endodermal fates. By contrast, in partnership with eomesodermine (EOMES) and activin-SMAD2/3 signalling, BRA proves to be sufficient, but not necessary, to activate endodermal gene expression.

6We have previously optimised chemically defined conditions that cause hESCs to differentiate as progenitors of endoderm or of mesoderm. When hESCs are induced to differentiate (36 h of treatment with Fgf2 and the PI3 kinase inhibitor Ly294002) in an activin A-based medium (called FLyA), they upregulate endoderm markers such as CER1, FOXA2, GSC and SOX17. By contrast, hESCs similarly induced to differentiate in a BMP4-based medium (called FLyB) express mesodermal genes such as CDX2, TBX6, FOXF1 and BMP4. Thus, FLyA-treated hESCs resemble the anterior region of the early primitive streak, whereas FLyB-treated hESCs resemble the posterior region (Fig. 1A).

7Significantly, these two differentiated populations express different levels of BRA and EOMES, two T-BOX transcription factors that are key regulators of gastrulation. FLyA-treated cells are predominantly BRAlow/EOMEShigh and upregulate SOX17 (Fig. 1B,C; supplementary material Fig. S1A,B), whereas FLyB-treated cells are mainly BRAhigh/EOMESlow and upregulate CDX2 (Fig. 1B,C; supplementary material Fig. S1A,B). As expected, activin A or BMP4 treatment triggers the phosphorylation of their downstream targets, SMAD2/3 or SMAD1, respectively (Fig. 1B). This in vitro differentiation system allows us to capture and analyse the transient progenitor populations that give rise to two primary embryonic tissue lineages: endoderm and mesoderm.

8FLyA- and FLyB-treated hESCs express, respectively, genes that are characteristic of anterior and posterior regions of the primitive streak. They also express different levels of BRA, which reaches its peak at around 36 h of differentiation. We used this differentiation system to identify genomic targets of BRA by ChIP-seq and to ask whether these differ between endoderm and mesoderm progenitors (Fig. 2A).

9Two replicate ChIP-seq experiments, both with FLyA- and FLyB-treated cells, were carried out using two different anti-BRA antibodies, one from R&D Systems and the other from Santa Cruz Biotechnology. These replicates showed high correlation coefficients for BRA binding in both datasets (R=0.72 for FLyA conditions and R=0.73 for FLyB conditions) (supplementary material Fig. S2A). Overlapping ChIP-seq peaks (present using either antibody) are shown in supplementary material Table S1. However, the ChIP-seq experiments using the Santa Cruz antibody had lower signal-to-noise ratio and yielded lower ‘unique read’ counts. We therefore used the datasets generated with the R&D Systems antibody to perform all subsequent analyses.

10ChIP-seq analysis indicated that there are 25,836 BRA-binding events in FLyA-treated hESCs and 23,631 in FLyB-treated hESCs (Fig. 2B). Although there is considerable overlap between these two datasets (Fig. 2B; supplementary material Fig. S2B), a large proportion of peaks are unique to each cell population: 68% in the FLyA condition and 64% in the FLyB condition (Fig. 2B). For example, BRA binds in the vicinity of genes such as: (1) CDC6, FOXI1, HES5 and JAG1 when cultured in FLyA but not so strongly in FLyB (Fig. 2D; left blue rectangle, FLyA>FLyB); (2) AXIN2, MESP2, RIPPLY1 and TBX6 when cultured in FLyB but not so strongly in FLyA (Fig. 2D; right red rectangle, FLyA<FLyB); and (3) CER1, FOXA2, HOXB13 and MESP1 when hESCs are cultured in FLyA or FLyB (Fig. 2D; central purple rectangle, FLyA≈FLyB). Moreover, the ‘fold enrichment’ correlation coefficient for peaks detected in both FLyA and FLyB conditions was low (R=0.22) (Fig. 2C), further underscoring the differences in BRA binding between FLyA- and FLyB-treated hESCs. Importantly, these differences in BRA binding in FLyA or FLyB conditions were confirmed by ChIP-qPCR on a set of selected targets (supplementary material Fig. S2C,D). These observations reveal that BRA has distinct genome-wide binding landscapes in hESC-derived endoderm and mesoderm progenitors.

11Using a ‘nearby gene’ peak annotation approach (up to 50 kb on either side of a gene), BRA binding was detected in the vicinity of 10,074 genes in FLyA-treated hESCs and 8983 genes in FLyB-treated hESCs (Fig. 3A; supplementary material Table S1). As suggested by the overlap analysis (Fig. 2A), BRA binds to many of its target genes in a cell type-specific manner. Thus, in the FLyA dataset, 35% of targets were unique to this condition and 27% were unique to the FLyB dataset (Fig. 3A). These percentages are lower than those concerning the binding overlap (Fig. 2B) because several peaks are often located around a single locus (in both promoter and enhancer regions) and because we excluded the most distal intergenic peaks, the assignment of which to the nearest gene can lead to a higher false discovery rate (data not shown).

12Recently, BRA genomic binding was investigated by ChIP-seq in hESC-derived ‘mesendoderm’ progenitors (12 h treatment with WNT3A and activin A). This dataset shows a substantial overlap with ours both in terms of binding regions (supplementary material Fig. S3A-C) and target genes (Fig. 3A). However, many unique binding events are detected in each dataset, further suggesting that the regulatory role of BRA is context dependent.

13To investigate the developmental significance of our findings, we performed gene ontology analyses using GREAT (supplementary material Table S2), distinguishing between FLyA-only putative targets (Fig. 3B), FLyB-only putative targets (Fig. 3C), and targets common to both FLyA and FLyB (Fig. 3D). In all subsets, there was clear enrichment for genes that are co-expressed with Bra during mouse embryogenesis and involved in all aspects of its complex mutant phenotype (Fig. 3, Table 1; supplementary material Table S2).

14The gene ontology analysis of the recent human dataset included few highly significant terms (P-value <1×10−9) (supplementary material Table S2). However, the top term in the ‘mouse expression’ category was ‘TS9_primitive streak’ (P=4.92×10−9) and we note that several BRA targets that are common to the Tsankov et al. dataset and our dataset make biological sense (Table 1).

15Strikingly, whereas FLyA-only BRA targets were enriched for genes expressed in anterior primitive streak derivatives, such as the node and endoderm (Fig. 3B), FLyB-only BRA targets were enriched for genes expressed in mid/posterior primitive streak derivatives, such as tail mesoderm, lateral plate mesoderm, extraembryonic mesoderm, allantois and somitic/paraxial mesoderm (Fig. 3C). BRA targets that were common to both FLyA and FLyB conditions are expressed in all germ layers (Fig. 3D).

16Overall, our results show that in FLyA- and FLyB-treated cells, BRA binds to the genome in a manner that is correlated with the transcriptional and developmental identity of each cell population, endoderm and mesoderm, respectively.

17Having identified genomic targets of BRA in differentiating hESCs, we asked whether BRA-binding regions (200 bp sequences centred on ChIP-seq peaks) were enriched for specific DNA sequences. To this end, we performed de novo DNA motif analyses using the MEME suite. For FLyA-treated cells, these revealed enrichment of a motif resembling the T-BOX consensus-binding sequence in 63% of peaks (Fig. 4A). Motifs for other protein families were also enriched in BRA FLyA peaks (Fig. 4A; supplementary material Fig. S4A), including FOX, GATA/GSC, SMAD/ZIC, SOX, POU and STAT, suggesting that BRA might interact with members of these families. Likely candidates, based on their expression patterns in vivo and in FLyA-treated cells, include FOXA2, GATA4/6, GSC, SMAD2/3, SOX17 and POU5F1/OCT4 (supplementary material Fig. S4C,D). Significantly, we were able to show by co-immunoprecipitation that both BRA and EOMES interact with SMAD2/3 (Fig. 4B), suggesting cooperation with activin signalling. BRA-binding peaks show a remarkable overlap with those of EOMES and of SMAD2/3 in hESC-derived mesendodermal progenitors (Fig. 4C,D), and share with BRA many putative common target genes (Fig. 4E). Importantly, the genomic binding of EOMES and SMAD2 in FLyA-treated cells was confirmed by ChIP-qPCR on a set of selected BRA target regions (supplementary material Fig. S4B).

18Both EOMES and SMAD2/3 are essential for proper expression of endoderm markers in differentiating hESCs, but little is known about the role of BRA in this context. We therefore sought to discover which putative BRA targets require normal BRA levels for their correct expression in FLyA-treated hESCs. BRA knockdown (shRNA KD) hESCs were compared with their wild-type counterparts by transcriptional profiling at 36 (Fig. 4F; supplementary material Table S3) and 72 h of differentiation, in more mature endoderm cells (supplementary material Fig. S4F and Table S3). At the protein level, BRA was virtually absent in BRA knockdown cells (supplementary material Fig. S4E). Misregulated transcripts in BRA knockdown cells were compared with BRA-bound genes in the FLyA condition (Fig. 4G; supplementary material Fig. S4G). Up- and downregulated genes, both at 36 and at 72 h, showed enrichment for BRA targets (P<0.01, Pearson's Chi-squared).

19We proceeded to perform gene ontology analysis of misregulated BRA targets using GREAT (supplementary material Table S4). Although BRA was necessary for the upregulation of many target genes that are developmentally important and expressed under FLyA conditions (Fig. 4H; supplementary material Fig. S4H), this subset appears to be enriched for genes expressed in neurectoderm derivatives post-gastrulation (supplementary material Table S4), which are not gene ontology categories that are developmentally relevant to FLyA-treated cells. Interestingly, however, key genes involved in epithelial-to-mesenchymal transition (EMT) were downregulated at 36 h (SNAI2, TWIST2 and FOXD3), together with the upregulation of CDH1/E-cadherin, a classical hallmark of impaired EMT.

20Strikingly, the expression of many genes expressed in the anterior primitive streak and involved in endoderm formation was either unaltered or even upregulated in BRA knockdown cells (P=3.49×10−11 at 36 h, P=1.48×10−11 at 72 h; supplementary material Table S4). The latter included genes such as GSC, GATA6, HHEX, LEFTY1/2, CXCR4, OTX2 and LHX1 (Fig. 4H; supplementary material Fig. S4H). These observations were confirmed by qRT-PCR in another BRA knockdown clone (81% BRA knockdown efficiency in FLyA treatment at 36 h) when compared with a scrambled (mock shRNA) control line (supplementary material Fig. S4I). Under FLyA conditions, many of these genes were not only bound by BRA but also by EOMES and SMAD2/3 (Fig. 4I).

21Together, these data show that the genome-wide binding of BRA in endoderm progenitors broadly overlaps with that of EOMES and SMAD2/3, both of which are essential regulators of endoderm formation. However, BRA is not necessary for the expression of most key endoderm markers.

22We also performed de novo DNA motif analysis using the FLyB BRA peak dataset. As observed with the FLyA dataset, this analysis revealed enrichment (69% of peaks) of the T-BOX consensus binding sequence (Fig. 5A). Motifs characteristic of other protein families were also found in BRA FLyB peaks (Fig. 5A; supplementary material Fig. S5A), including PU-BOX, POU, SMAD/ZIC, FOX, KLF and GATA, again suggesting that BRA might interact with members of these families. Likely candidates, based on their expression patterns in vivo and in FLyB-treated cells, include SMAD1, FOXF1, GATA2/3 and POU5F1/OCT4 (supplementary material Fig. S5B,C). Previous work has shown that a BRA orthologue, Xbra, directly interacts with Smad1 in Xenopus embryos. Significantly, we were able to show by co-immunoprecipitation that both BRA and EOMES interact with SMAD1 (Fig. 5B), suggesting cooperation with BMP4 signalling. BRA binding peaks also showed close proximity to those of EOMES in hESC-derived mesendodermal progenitors (Fig. 5C), with a substantial number of common putative target genes (Fig. 5D). Importantly, the genomic binding of EOMES and SMAD1 in FLyB-treated cells was confirmed by ChIP-qPCR on a set of selected BRA target regions (Fig. 5E).

23We next asked which putative BRA targets require normal BRA levels for their correct expression in BMP4-treated cells by comparing the transcriptomes of BRA knockdown (shRNA KD) hESCs with their wild-type counterparts at 36 h of FLyB differentiation (Fig. 5F; supplementary material Table S3) and at 72 h of FLyB/FB differentiation (supplementary material Fig. S5D and Table S3), when cells resemble extraembryonic and lateral plate mesoderm. Misregulated transcripts in BRA knockdown cells were compared with BRA-bound genes in the FLyB condition (Fig. 5G; supplementary material Fig. S5E). Up- and downregulated genes, both at 36 and 72 h, showed enrichment for BRA targets (P<0.01, Pearson's Chi-squared).

24We then proceeded to perform gene ontology analysis of misregulated BRA targets using GREAT (supplementary material Table S4). Significantly, BRA was necessary for the normal upregulation of many genes expressed in BMP4-treated cells (at 36 and 72 h) that are essential for mesoderm development. These include ALDH1A2 (RALDH2), AXIN2, CDX2, FOXF1, KDR, LMO2, MSGN1, MEIS1, TBX6 and WNT3A (Fig. 5H; supplementary material Fig. S5F). Indeed, gene ontology analysis shows enrichment for phenotypes that are reminiscent of BRA mutant embryos, such as ‘abnormal somite size’ (P=9.26×10−24), ‘abnormal tail development’ (P=1.90×10−22), ‘abnormal gastrulation’ (P=1.56×10−17) and ‘abnormal vascular development’ (P=3.29×10−16) (supplementary material Table S4). Interestingly, several HOX genes were also downregulated upon BRA knockdown at 36 h (Fig. 5H) (P=1.30×10−12; supplementary material Table S4), an effect that was even more marked at 72 h of differentiation (supplementary material Fig. S5F) (P=3.26×10−6; supplementary material Table S4).

25Noticeably, the expression of many endoderm regulators was upregulated in BRA knockdown cells grown in FLyB conditions (P=3.87×10−19 at 36 h, P=8.36×10−16 at 72 h; supplementary material Table S4). These include CER1, CYP26A1, EOMES, FOXA2, GSC, GATA6, HHEX, LEFTY1/2, LHX1, MIXL1, OTX2 and SOX17 (Fig. 5H; supplementary material Fig. S5F). These observations were confirmed by qRT-PCR in another BRA knockdown clone (89% BRA knockdown efficiency in FLyB treatment at 36 h) when compared with a scrambled (mock shRNA) control line (supplementary material Fig. S5G). Strikingly, BRA was bound in the vicinity of some of these genes in regions only detected in the FLyB condition (e.g. LHX1, Fig. 5I).

26Together, these data show that the genome-wide binding of BRA in mesoderm progenitors overlaps with that of EOMES. However, unlike EOMES, BRA was necessary for the expression of many genes involved in mesoderm formation, while simultaneously repressing the expression of endoderm markers.

27Having established that BRA is required for the expression of many mesodermal but not endodermal genes, we asked whether BRA overexpression (BRA OE) (supplementary material Fig. S6A,B) in hESCs was sufficient to up- or downregulate its genomic targets. Bearing in mind that BRA has different targets in activin- or BMP4-treated cells, we analysed the phenotype of control and BRA OE hESC subclones in different signalling environments. All cells were grown in the presence of Fgf2 and Ly294002 as a basal differentiation medium (FLy) with the addition of activin (FLyA), BMP4 (FLyB), SB431542 (FLyS) to block activin-SMAD2/3 signalling or noggin (NOG) (FLyN) to block BMP4-SMAD1 signalling (Fig. 6A,B).

28Overexpression of BRA caused the upregulation of both endodermal (CER1, FOXA2, SOX17, MIXL1) and mesodermal (ALDH1A2, BMP4, CDX2, FOXF1, HOXD1, LMO2) target genes. Strikingly, however, strong upregulation of genes characteristic of endoderm/anterior primitive streak required activin signalling (Fig. 6C). Genes characteristic of mesoderm/posterior primitive streak either required or benefitted from BMP4 signalling and, interestingly, in some cases, activin antagonism (Fig. 6D). These experiments reveal that BRA expression alone is not sufficient to fully upregulate some of its target genes to wild-type levels and that it requires appropriate signalling cues to activate the expression of its genomic targets.

29BRA putative target genes have been identified by ChIP-seq in activin-treated mouse embryoid bodies. Our analysis reveals a considerable overlap between our datasets and that of Lolas et al., as 71% of the mouse targets are included in the human datasets (Fig. 7A). Indeed, many of these conserved targets make biological sense (Table 1). Importantly, our datasets also reveal the identity of BRA targets not identified in mouse (Table 1), including ALDH1A2, HES7, KDR, LMO2 and MEOX1.

30In order to investigate whether our observations made with BRA knockdown hESCs in vitro are also relevant in vivo, we analysed the expression of BRA targets in mouse Bra mutants compared with wild-type embryos at E7.5-8.0 (Fig. 7B; supplementary material Table S5) (fold change>2; RNA-seq data from ). E7.5 is a comparable developmental stage to our 36 h hESC differentiation protocol. These analyses show that orthologues of key BRA mesodermal targets such as Aldh1a2, Cdh2 (N-cadherin), Foxd3, Hes7, several Hox genes, Msgn1, Tbx4, Tbx6 and Wnt5b are downregulated in Bra−/− mouse embryos (supplementary material Table S5). In parallel, targets such as Cdh1, Cer1 and Eomes are upregulated in Bra−/− mouse embryos (supplementary material Table S5). Other anterior/endoderm markers, including Cyp26a1, Foxa2, Gata6, Otx2 and Sox17, are also upregulated in Bra−/− mouse embryos but to a lesser extent (fold change <2). More broadly, 823 putative BRA targets exclusively found in the human ChIP-seq datasets show misregulation in mouse mutant embryos (Fig. 7C), underscoring the value of our new data.

31However, RNA-seq data analysis does not provide information on spatial expression patterns in vivo. Thus, we analysed the expression pattern of two key BRA target genes at the protein level in both wild-type and Bra mutant mouse embryos: the mesoderm regulator Cdx2 and the endoderm regulator Foxa2. In Bra-null embryos, Cdx2 mRNA is downregulated 1.22 fold, while Foxa2 mRNA is upregulated 1.30 fold when compared with wild type.

32In agreement with our prediction, Cdx2 nuclear expression, which is visible in a small group of cells in the proximal primitive streak where extraembryonic mesoderm is being formed (Bra+/+, white arrowhead, Fig. 7D) is completely lacking in Bra mouse mutants (Bra−/−, white arrowheads, Fig. 7D). Foxa2 expression, which marks the distal (future anterior) region of the primitive streak (Bra+/+, white bracket, Fig. 7D) and is also expressed in extraembryonic endoderm is clearly present in Bra-null embryos (Bra−/−, white brackets, Fig. 7D; supplementary material Fig. S7A), and is seemingly upregulated in the mid/proximal primitive streak. Together, the results obtained with mouse embryos that lack Bra are consistent with our findings in differentiating hESCs: BRA upregulates key mesoderm genes and is dispensable for the expression of several endoderm markers.

33BRA ChIP-seq analysis of hESCs differentiating in FLyA-, FLyB- or WNT3A- and activin-containing media showed similar binding sites but also significantly different DNA occupancy levels of BRA. Chromatin accessibility is suggested to be a prerequisite for many sequence-specific transcription factors such as BRA to directly bind DNA. Thus, differences in BRA binding may arise from differential nucleosome occupancy, distinct histone modifications, or the presence or absence of specific protein partners. It may be that activin, BMP4 and WNT3A signalling cascades interact with chromatin remodelling complexes so as to affect DNA accessibility in a cell type-specific manner.

34It is unlikely that the distinct BRA binding in FLyA or FLyB occurs solely because BRA protein levels are lower in hESCs cultured in FLyA rather than FLyB, as there are many binding sites in FLyA-cultured cells with equal or even higher occupancy levels than in FLyB-culture cells. Furthermore, the strong correlation of BRA differential binding with the distinct transcriptional identities of FLyA- or FLyB-treated cells suggests that these events are biologically meaningful. For example, in the FLyB condition (BRAhigh), BRA binds close to many genes involved in posterior mesoderm development, consistent with fate-mapping studies showing that this tissue emerges from the posterior primitive streak, where Bra levels are higher. Gene ontology analysis reveals that BRA binds in the vicinity of genes involved in specific developmental contexts and that are expressed in tissues where BRA function is essential. These contexts include: epithelial-to-mesenchymal transition; node and notochord formation; the establishment of left-right asymmetry (including heart looping), axial extension and somitogenesis; and cardiovascular development (Fig. 3).

35Interestingly, the most highly enriched category in FLyA and FLyB common peaks was ‘neural tube development’. Indeed, BRA seems to repress neurectoderm genes in the tail bud ‘stem cell’ pool that gives rise to posterior structures such as somites, the notochord and the neural tube.

36We observed a substantial overlap between BRA targets identified in FLyA- or FLyB-treated hESCs and those identified in WNT3A/activin-treated hESCs or identified in activin-treated mouse embryoid bodies – the latter suggesting a considerable degree of evolutionary conservation. The fact that the BRA targets we uncovered differ from those of Tsankov et al. could be due to distinct signalling environments, developmental stages, ChIP protocols, antibodies used or other technical issues. We note that both in Tsankov et al. and in Lolas et al. the authors used activin treatment to induce ESC differentiation, thus probably enriching for endoderm progenitors. Indeed, their datasets show a greater overlap with FLyA treatment, rather than with FLyB. Possibly due to this, Tsankov et al. and Lolas et al. did not detect BRA binding in the vicinity of many important mesodermal loci (Table 1). These include key genes in posterior mesoderm and cardiovascular development such as ALDH1A2, HAND1, HES7, KDR, LMO2 and MEOX1, which thus seem to be human specific and unique to our datasets. In summary, our ChIP-seq datasets, together with those recently generated, have greatly expanded our knowledge of BRA targets in mammals and now provide a rich resource for future studies to discover new regulators of human gastrulation.

37In addition to the T-BOX consensus binding motif (T-site), DNA motif analyses revealed enrichment for binding sites for protein families with important roles during embryonic development such as FOX, POU, SOX and GATA. Of the many motifs enriched in BRA-binding sites, we focused on the SMAD/ZIC motif. SMAD proteins are frequently recruited by master transcription factors to regulate cell type-specific gene regulatory programs. Indeed, we found that BRA interacts with SMAD2/3 in FLyA-treated cells and with SMAD1 in FLyB-treated cells. Overexpression of BRA in hESCs upregulates several mesodermal and endodermal target genes, but this effect is particularly potent in the presence of either BMP4 or activin in the culture medium, thus reinforcing the idea that BRA collaborates with these SMAD signalling cascades to regulate its targets in a cell type-specific manner.

38BRA binding in FLyA-treated hESCs overlaps very significantly with EOMES binding. Like other T-BOX proteins in zebrafish (Bra, Ntl, Spt and Tbx6), in Xenopus (Xbra, Xbra3, VegT and Eomes) or in the mouse (Bra and Eomes), human EOMES and BRA might be functionally interconnected in distinct developmental contexts.

39Interestingly, the POU core motif was also enriched in BRA peaks. Although OCT4 is traditionally regarded as a pluripotency regulator, it also acts as lineage specifier. Indeed, OCT4 is expressed throughout the primitive streak during gastrulation, together with Bra. Although we (data not shown) and others have been unable to detect a direct physical interaction between BRA and OCT4, it is tempting to speculate that BRA might cooperate with OCT4 and/or other pluripotency factors during the early stages of gastrulation by co-regulating some common target genes, possibly through dynamic sequential binding to these loci. Recent ChIP-seq data analysis seems to corroborate this hypothesis.

40Comparison of our ChIP-seq data with results obtained from transcriptional profiling of BRA knockdown cells showed that, although significantly enriched for BRA targets, around half of up- and downregulated transcripts do not appear to be direct targets of BRA and, conversely, that the expression of most putative direct targets is not affected by loss of BRA function. This outcome may derive in part from the inefficiency of our shRNA approach to completely eliminate BRA (supplementary material Fig. S4E), but it is also likely to reflect the functional complexity of transcriptional networks. For example, BRA shares a vast number of targets with EOMES. This complicates the analysis of BRA-driven gene regulatory networks because a redundancy with EOMES might mask the number of functional target genes of BRA, as seen in Xenopus embryos.

41The number of functional targets of BRA might well be extended if BRA knockdown hESCs were tested in other differentiation protocols where BRA has proposed functions, including primordial germ cells, axial mesoderm, and paraxial and cardiac mesoderm. Indeed, our BRA ChIP-seq data revealed that BRA binds in the vicinity of genes involved in the formation of these lineages (Fig. 3).

42Importantly, the comparison of our human datasets with mouse Bra ChIP-seq data and RNA-seq from Bra-null gastrulae provides in vivo validation for key BRA targets and underscores the novelty of our work, because around 800 BRA putative targets exclusively identified in our ChIP-seq datasets are either down- or upregulated in Bra mutant embryos (Fig. 7C).

43Gene ontology analysis of transcriptional changes for BRA target genes observed in BMP4-treated BRA knockdown cells revealed a clear phenotypic scenario: while several BRA target genes involved in mesoderm formation were markedly downregulated, many genes important for endoderm development were upregulated. In agreement with our results, it has been shown that Bra, in collaboration with Mixl1, is able to repress endoderm or anterior primitive streak markers such as Gsc and Pdgfra. However, the mechanism by which BRA regulates endoderm differentiation is likely to be complex, because elevated levels of BRA in FLyA-treated cells cause an increase in expression of endodermal marker genes (Fig. 6). Similar results have been obtained by, who showed that BRA knockdown and BRA overexpression both cause upregulation of the endodermal markers GSC and SOX17 in hESCs.

44Interestingly, the forced expression of several BRA orthologues in Xenopus animal caps leads to mesoderm but not to endoderm formation, except when using the Drosophila and Ciona orthologues. In these two organisms, the N-terminal domain of Bra, which is essential for the interaction with SMAD1, is not conserved. Indeed, a truncated form of Xbra, which can bind DNA but is unable to interact with SMAD1, can induce the expression of Gsc, an anterior/endoderm marker. These observations suggest that Bra induces mesoderm formation when cooperating with SMAD1, and endoderm formation when the interaction with SMAD1 is absent. This may account for the aforementioned paradox that BRA can both repress and activate endoderm gene expression.

45By contrast, gene ontology analysis of downregulated transcripts in activin-treated BRA knockdown cells at either 36 or 72 h of differentiation did not reveal any significant biological insights into the role of BRA as an activator in FLyA-treated cells. However, one important aspect that emerged from this profiling was the observation that genes classically involved in EMT were misregulated in FLyA-treated BRA knockdown cells and also in vivo. This provides a new molecular insight into the known migratory defects of Bra knockout cells and agrees with the EMT-promoting role of BRA in human cancer samples. Indeed, the focus on BRA as an important player in cancer biology, particularly in chordoma, has been growing over the past decade.

46Our results provide the basis for the following model explaining the establishment of the gene expression patterns in the primitive streak that result in cells acquiring different fates (Fig. 8). Anterior early primitive streak cells, which give rise to definitive endoderm and axial mesoderm, experience high levels of NODAL signalling. This induces the expression of high levels of EOMES and low levels of BRA (Fig. 1). EOMES then cooperates with NODAL-SMAD2/3 signalling to induce the expression of anterior markers such as CER1, FOXA2 and SOX17, while repressing posteriorly expressed genes. Interestingly, by targeting the same genomic locations as EOMES, BRA is also able to activate the expression of endodermal genes, by cooperating with SMAD2/3 (Fig. 4). However, BRA is not required for anterior primitive streak gene expression (Fig. 4) where EOMES is the main player.

47At the opposite end of the early primitive streak, which gives rise to extraembryonic and posterior mesoderm, cells experience high levels of BMP4 signalling. This induces high levels of BRA and low levels of EOMES (Fig. 1). BRA, in cooperation with BMP4-SMAD1 signalling, then induces the characteristic expression of posterior markers such as HOX genes, CDX2, TBX6 and KDR, while repressing the expression of anterior genes (Fig. 5). Whether EOMES can activate mesodermal gene expression in collaboration with BMP4-SMAD1 signalling remains an unresolved issue.

48Together, our findings illuminate the function of BRA in the human species. We have shown that BRA is indeed both necessary and sufficient to regulate the transcription of many of its putative targets, which are key players during mesoderm or endoderm development. Moreover, our study shows that the regulatory role of BRA is context dependent, thus establishing an intimate collaboration of BRA with SMAD2/3 in an activin/NODAL-dominated context and with SMAD1 in a BMP4-dominated context.

49Our study thus reinforces the value of hESCs as tools to model human embryonic development and emphasizes the importance of analysing the functions of individual members of complex transcription factor networks in distinct cellular and signalling contexts.

50Human ESCs (H9 line, WiCell) were grown in a chemically defined medium (CDM) as previously described. For differentiation, cells were grown in CDM containing PVA instead of BSA and supplemented as described in the text and in the methods in the supplementary material. Transfection and selection of stable knockdown or overexpression clones were carried out as described in the methods in the supplementary material.

51Total RNA was extracted using the RNeasy Mini kit (Qiagen). Half a microgram of RNA was reverse transcribed using the Maxima First Strand cDNA Synthesis Kit for RT-qPCR (Thermo Scientific). Quantitative reverse transcription polymerase chain reaction (qRT-PCR) mixtures were prepared using Fast SYBR Green Master Mix (Applied Biosystems). PCR reactions were performed in a 7500 Fast Real-Time PCR System (Applied Biosystems). All procedures followed manufacturer's instructions. See also methods and Table S6 in the supplementary material.

52Nuclear extract preparation, co-immunoprecipitation and western blotting are described in the methods in the supplementary material.

53Cells were fixed for 10 min at room temperature in 4% paraformaldehyde and immunostained as described in the methods in the supplementary material. Fluorescent images were taken using an Olympus IX71 microscope.

54Cells were fixed and immunostained using the Cytofix-Cytoperm kit (BD Biosciences) as described in the methods in the supplementary material. Cells were analysed using a Beckman Coulter CyAnADP flow cytometer and FlowJo software (Becton Dickinson).

55Chromatin immunoprecipitation was performed as previously described, with some modifications. Sequencing libraries were prepared using the ChIP-seq DNA sample kit (Illumina, IP-102-1001) with some modifications and sequenced with a Genome Analyzer II (Illumina). Data were further processed using the MACS peak finder. DNA motifs were analysed using the MEME suite. ChIP-seq data were visualised using the UCSC Genome Browser. Gene ontology analysis was performed using GREAT with default parameters. Data are available in the GEO database (www.ncbi.nlm.nih.gov/geo) under accession number GSE60606. See methods in the supplementary material for further details.

56Sample preparation was performed according to manufacturer's instructions (Illumina). Labelled extracts were hybridised to whole-genome bead array (HumanWG-6 v3.0 Expression BeadChip) on an Illumina BeadArray reader. Gene expression heat-maps were generated by importing subsets of processed microarray data as described in the methods in the supplementary material. Data are available in the ArrayExpress database (www.ebi.ac.uk/arrayexpress) under accession numbers E-MTAB-2912 and E-MTAB-464.

57T/+×T/+ mice were mated for embryo collections. Late gastrulae were dissected (E6.75-7.0) for further analysis. All mouse studies were performed under a UK Home Office project license and complied fully with the UK Animals (Scientific Procedures) Act 1986 as implemented by the University of Cambridge and the Medical Research Council. Embryos were fixed for 20-30 min at room temperature in 4% paraformaldehyde (PFA) supplemented with 0.1% Tween-20 (Sigma) and 0.01% Triton-X100 (Sigma). Immunostaining was performed following standard procedures as described in the methods in the supplementary material. Fluorescent images were captured using a Zeiss LSM 710 microscope.