A gene regulatory network for neural induction

In the chick embryo, a graft of Hensen’s node to the inner area opaca, an extraembryonic region of competent non-neural ectoderm (Gallera and Ivanov, 1964; Dias and Schoenwolf, 1990; Storey et al., 1992; Streit et al., 1998), induces the formation of a mature, patterned neural tube after about 15 hr of culture from the host ectoderm (Figure 1A). To characterize the events over this period, the induction of several neural markers was assessed. OTX2 is first expressed in the pre-primitive-streak stage epiblast at EGKXII-XIII (Albazerchi and Stern, 2007; Pinho et al., 2011) and is induced by 7/7 grafts after 3 hr (Figure 1B–F). SOX2 is first expressed at HH4+/5 (Rex et al., 1997; Streit et al., 1997; Streit et al., 1998; Sheng et al., 2003) and requires 9 hr for induction (8/8) (Figure 1G–K). SOX1 starts to be expressed weakly in the neural plate around HH7-8 (Stavridis et al., 2010; Uchikawa et al., 2011) and is induced by 50% of grafts after 12 hr (Figure 1L–P). This time course confirms that a sequence of events occurs in response to a grafted node over 0–12 hr (Figure 1—figure supplement 1A), culminating in the acquisition of neural plate/tube identity (Streit et al., 1997; Streit et al., 1998; Streit and Stern, 1999; Streit et al., 2000; Sheng et al., 2003; Stern, 2005; Albazerchi and Stern, 2007; Papanayotou et al., 2008; Gibson et al., 2011; Pinho et al., 2011; Papanayotou et al., 2013; Trevers et al., 2018).

Figure 1 with 1 supplement see all

Download asset Open asset

Figure 1—source data 1

RNAseq galGal3 analysis key and differentially expressed transcripts.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig1-data1-v2.xlsx

Download elife-73189-fig1-data1-v2.xlsx

Figure 1—source data 2

RNAseq galGal4 analysis key and differentially expressed transcripts.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig1-data2-v2.xlsx

Download elife-73189-fig1-data2-v2.xlsx

Figure 1—source data 3

Nanostring probe codeset.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig1-data3-v2.xlsx

Download elife-73189-fig1-data3-v2.xlsx

To uncover genes whose expression changes over this period, we performed RNAseq of the responding tissue at three time points following a graft of Hensen’s node: 5 hr (to identify early ‘pre-neural’ responses), 9 hr (when SOX2 expression defines neural specification), and 12 hr, as the host tissue starts to express SOX1. Hensen’s node was grafted from HH4- chick donors to a region of competent epiblast (inner area opaca) of chick hosts at the same stage. Embryos were cultured for the desired period of time, after which the graft was removed and the underlying ‘induced’ tissue was collected. Control (‘uninduced’) ectoderm from the corresponding region on the contralateral side of the same embryos (Figure 1Q) was also collected, as well as competent area opaca at HH4-, representing a ‘0 hr’ starting point for the time course. When ‘induced’ and ‘uninduced’ reads at each time point are compared by DESeq analyses, 8673 differentially expressed transcripts were identified (see volcano plots, Figure 1R–T). Of these, 4130 were upregulated (enriched in ‘induced’ tissue) and 4543 were downregulated (depleted in ‘induced’ tissue) relative to the ‘uninduced’ counterpart.

To construct a GRN we focused on 482 transcription factors and chromatin modifiers that are differentially expressed within these data. Of these ‘transcriptional regulators’, 255 are upregulated, 203 are downregulated, and 24 have more complex expression dynamics (Figure 1U and Figure 1—figure supplement 1B). Grouping transcripts based on the timing of their response (Figure 1V–W) reveals that 120 are differentially expressed throughout (53 upregulated and 67 downregulated), while others are associated with particular time points. Therefore, transcriptional responses to organizer grafts involve changes in the expression of many genes over a relatively short period.

Due to this complexity, we chose to increase the time resolution of the analysis. NanoString nCounter was used to quantify the gene expression changes of transcriptional regulators at six time points (1, 3, 5, 7, 9, and 12 hr after the node graft) in ‘induced’ tissue compared to ‘uninduced’ ectoderm. By consolidating the data from RNAseq and NanoString, a set of refined expression profiles was established for 213 transcriptional regulators – 156 that are enriched and 57 that are depleted in induced tissues (Figure 3—source data 1). These represent the core components of our GRN.

We next sought to identify the regions of non-coding chromatin that govern these transcriptional responses to signals from the organizer. Histone modifications can regulate gene expression by altering chromatin structure. For example, H3K27ac is associated with actively transcribed genes and their enhancers, whereas transcriptionally inactive regions are often marked by H3K27me3 (Tiwari et al., 2008; Heintzman et al., 2009; Creyghton et al., 2010; Kharchenko et al., 2011; Rada-Iglesias et al., 2011; Tolhuis et al., 2011; Zentner et al., 2011; Bonn et al., 2012). Therefore, ChIPseq and ATACseq (Buenrostro et al., 2013; Buenrostro et al., 2015; Corces et al., 2017) were performed on induced and uninduced ectoderm following 5, 9, or 12 hr of a node graft, to detect histone and chromatin conformation changes during neural induction.

H3K27ac and H3K27me3 enriched chromatin (ChIPseq) were identified genome-wide by comparison to matched genomic input samples. Chromatin sites were then categorized according to their histone signatures across induced and corresponding uninduced tissues at each time point (Figure 2—figure supplement 1A). Sites that become acetylated and/or demethylated in induced tissues, compared to uninduced, were considered to undergo ‘activation’ (Indices 1–3) while those that become deacetylated and/or methylated undergo ‘repression’ (Indices 4–6). Chromatin marked by both H3K27ac (activation) and H3K27me3 (repression) marks in either tissue were described as ‘poised’ (Indices 7–12). The remaining indices (13–16) define sites that do not change marks between the two conditions.

Chromatin sites undergoing activation are enriched for H3K27ac marks in induced tissues at each time point (Figure 2B–D). Sites of repression are more varied: those belonging to this category at 5 hr lose acetylation, whereas several sites become methylated in induced tissues at 9 and 12 hr. Very few sites are ‘poised’ at any time. Overall, as neural induction progresses, the number of sites undergoing activation increases, along with a reduction in those that do not change state. In contrast, ATACseq suggests that there is little difference between induced and uninduced samples in terms of chromatin accessibility (Figure 2—figure supplement 2) at 5 hr, but the differences become more marked at later time points. At 12 hr, the ATACseq profile resembles that of the endogenous neural plate of the embryo (Figure 2—figure supplement 2).

Figure 2 with 3 supplements see all

Download asset Open asset

Figure 2—source data 1

ChIPseq tissue dissociation details.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig2-data1-v2.xlsx

Download elife-73189-fig2-data1-v2.xlsx

Figure 2—source data 2

ATACseq indices.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig2-data2-v2.xlsx

Download elife-73189-fig2-data2-v2.xlsx

Next, we focused on epigenetic changes associated with the 213 genes encoding transcriptional regulators for which refined expression profiles were established by NanoString analysis (see above). ‘Constitutive’ CTCF-bound sites (putative insulators) flanking these genes were obtained from chicken CTCF-ChIPseq data (Khan et al., 2013; Kadota et al., 2017). They were used to predict the boundaries within which we considered H3K27 marks as being associated with the gene of interest; these regions (‘loci’) are up to 500 kb in length and may contain several genes (Figure 2—figure supplement 1B). Within each locus, the H3K27ac or H3K27me3 enriched peaks were categorized as before (Figure 2—figure supplement 1A). The 213 selected loci contain a total of 6971 sites that change in response to a node graft (Indices 1–10; Figure 2—figure supplement 1A). We noticed that the flanking regions (proximal promoters) of transcriptional regulators whose expression increases are enriched for H3K27ac marks in induced tissue (Figure 2—figure supplement 3). In contrast, there is no difference between H3K27ac and H3K27me3 marks around the flanking regions of repressed transcriptional regulators at 5 and 9 hr – they only become enriched for methylation at 12 hr. Our analysis therefore identifies putative regulatory elements that are controlled epigenetically during the process.

Having identified transcriptional responses to signals from the organizer and their accompanying chromatin changes, we combined them to construct a GRN to describe the time course of these events and to illustrate predicted interactions between transcriptional regulators. The open source platform BioTapestry (Longabaugh et al., 2005; Longabaugh et al., 2009; Paquette et al., 2016) (http://www.biotapestry.org) is used extensively to provide a visually intuitive representation of developmental networks (Davidson, 1990; Arnone and Davidson, 1997; Davidson et al., 1998; Betancur et al., 2010; Simões-Costa and Bronner, 2015; Thiery et al., 2020). Importantly, this software allows dynamic changes in the interactions between regulators and their target genes to be represented.

A BioTapestry model for regulatory gene interactions

To integrate these complex time course data into a GRN that models the interactions during neural induction, we developed a custom computational pipeline (Figure 3). The 213 transcriptional regulators previously identified were selected as candidate components of the network because of their responses to a node graft in both RNAseq and NanoString data, and their associated chromatin changes in ChIPseq. Interactions were predicted by screening putative regulatory sites that belong to Indices 1–3 and 7 (activation) (Figure 2—figure supplement 1C–E) for the presence of putative binding sites for transcription factors that are differentially expressed in the dataset, and expressed at the appropriate time point to fit with a putative regulatory role (i.e. simultaneously or just before) for the target locus.

Figure 3

Download asset Open asset

Figure 3—source data 1

Integrated GRN gene expression profile.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig3-data1-v2.xlsx

Download elife-73189-fig3-data1-v2.xlsx

We then used the regulatory rules shown in Figure 4A–B to predict positive or negative regulatory events between regulators and their putative targets, overlaying the predicted transcription factor binding profiles with the time course of expression from RNAseq and NanoString. At each time point, changes in expression of a putative regulator that could lead to the same change in a candidate target are represented as positive interactions. Negative, or inhibitory, interactions are predicted when changes in expression of the regulator could cause the opposite change in a downstream target. Genes lacking both regulatory inputs and outputs with other genes in the network were excluded.

Figure 4

Download asset Open asset

Figure 4—source data 1

List of GRN transcription factor binding sites within regulatory sites at 5h, 9h and 12h following a node graft.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig4-data1-v2.xlsx

Download elife-73189-fig4-data1-v2.xlsx

The resulting GRN (Figure 4C), depicting interactions between 175 genes, was represented using BioTapestry (Supplementary file 1 and Supplementary file 2). Most targets are predicted to receive multiple positive and/or negative inputs, which are initiated at, and can act across, multiple time points. For example, the non-neural gene GATA2 is predicted to be repressed by TFAP2C after just 1 hr of a node graft, followed by increasing repression via HIF1A, MEIS2, and numerous other factors over the remaining time course. On the other hand, BLIMP1 is predicted to be induced by TFAP2C after 1 hr, alongside SNAI1, ETV1, ETV4, MGA, and SOX4 (Supplementary file 2). In total, 5614 interactions are predicted to occur between these components, highlighting the intricate and highly dynamic sequence of regulatory events that are triggered by signals from the organizer.

Incorporating individual regulatory sites and their dynamics into the network

BioTapestry networks usually represent each gene once, with multiple regulatory inputs, and each regulator as a single input into the target gene. However, it is now known that gene expression is controlled by multiple elements, each with characteristic spatial and temporal activity. In the past, such elements were identified within non-coding regions that are conserved across species by testing their ability to direct appropriate expression in vivo. This approach was used, for example, to identify many enhancers controlling SOX2 expression (Uchikawa et al., 2003; Okamoto et al., 2015) – the spatiotemporal pattern of expression of most developmentally regulated genes is likely to be controlled by multiple elements. Our epigenetic analysis predicts that in almost all cases there are indeed multiple putative regulatory sites associated with genes that change expression, which can be in different chromatin states (Figure 5A). To illustrate the contributions of different chromatin elements, we generated a subnetwork for BRD8, which is located in a relatively gene sparse locus with only a few putative regulatory elements associated with it. The changes that occur as sites undergo epigenetic ‘activation’ or ‘repression’ (using Indices 1–7; Figure 4A) were modelled across the 12 hr time course. Since BioTapestry lacks a notation for discrete regulatory elements of individual genes, the subnetwork shows seven candidate cis-regulatory regions of BRD8 (BRD8_site1-BRD8_site7 in Figure 5B) at 5, 9, and 12 hr represented as separate targets (Figure 5B–D). This depicts which specific elements contain binding sites for each regulating transcription factor. One element (BRD8_site1 in Figure 5B) is initially acetylated; however, it does not seeem to receive inputs from other GRN components. This site is then methylated at 9 and 12 hr as a number of other elements become active, presumably to stabilize and maintain BRD8 expression. In situ hybridization (Figure 5E) confirms that BRD8 is not expressed in the pre-primitive-streak stage epiblast (EGKXII) or in the prospective neural plate at HH3-4. BRD8 expression becomes enriched later in the neural plate and neural tube, as predicted by the network regulatory dynamics. This subnetwork illustrates the complexity of regulatory dynamics that underlie gene expression even at a relatively simple locus. Although it is not practical to represent the entire network including all the elements in their different states in this way using BioTapestry, all identified elements and their predicted activity are presented in the associated UCSC browser tracks (https://genome.ucsc.edu/s/stern_lab/Neural_Induction_2021), together with the predicted binding sites relevant to the network contained within each site; a full list is also given in Figure 4—source data 1.

Figure 5

Download asset Open asset

Prediction of core transcriptional regulators of the network

Genes with high ‘outdegrees’ and ‘betweenness centralities’ are often considered to be ‘core genes’ in a network. The ‘outdegree’ of a network component A is the number of outgoing interactions from component A to other components in the network (including component A itself), whereas the ‘betweenness centrality’ of a network component A is a measurement that captures the importance of component A based on the frequency at which interactions between pairs of other genes do so only by passing through component A (White and Borgatti, 1994). To predict key hubs in the network that may play more important roles during the process of neural induction, we calculated the ‘outdegrees’ and ‘betweenness centralities’ of network components across seven time points (0, 1, 3, 5, 7, 9, and 12 hr) (Figure 6). In the 0 hr network, ESRRG, GRHL2, and TFAP2A appear to be core genes in ectoderm cells prior to their exposure to a node graft (Figure 6B). Soon after grafting a node, these genes are immediately downregulated and repressed, with predicted repressive inputs from the core genes in the GRN. At the 1 hr time point following a node graft, TFAP2C, SOX4, and BLIMP1 start to be expressed and are predicted to act as core genes, acting as repressors of the next layer of core genes in ectoderm cells (Figure 6C–E). Other predicted neural induction core genes, such as SOX3, TCF12, and MYCN, are highlighted as key transcriptional regulators acting throughout the process of neural induction at multiple time points (Figure 6D–H), and are predicted to bind to the enhancers of many genes in the GRN.

Figure 6

Download asset Open asset

To what extent is the induction of an ectopic neural plate by an organizer graft comparable to the events of neural plate development in the embryo? To address this, we used three complementary approaches. First, we explored the spatial and temporal expression patterns of network components during normal embryo development. Then, we conducted scRNAseq analysis of stages of neural plate development concentrating on network components and their temporal hierarchical organization. Finally, we tested the activity of several of the putative regulatory elements in normal embryos to assess whether their activity in the neural plate in vivo resembles the hierarchy revealed by node grafts.

Spatiotemporal expression of GRN components during normal neural plate development

A previous study identified genes whose expression differs between EGKXII-XIII (pre-primitive-streak) epiblast, neural plate at HH6-7, and non-neural ectoderm (Trevers et al., 2018). When these embryonic tissues (Figure 7A) are compared to neural induction, we find that the pre-streak epiblast is most similar to induced ectoderm after 5 hr of a node graft, while 9 and 12 hr time points are more closely related to mature neural plate (Figure 7B). This encouraged us to generate a detailed spatiotemporal map of transcripts during embryonic neural plate development. Genes that are differentially expressed in induced tissues compared to uninduced in RNAseq were selected. In situ hybridization was performed for 174 transcriptional regulators (including 123 that are represented in the GRN) at four stages: pre-streak (EGKXII-XIII), primitive streak (HH3-4), head process/early neural plate (HH5-6), and neural fold/tube (HH7-9).

Figure 7 with 8 supplements see all

Download asset Open asset

Figure 7—source data 1

Sumary of changes for upregulated transcriptional regulators based on in situ expression patterns.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig7-data1-v2.xlsx

Download elife-73189-fig7-data1-v2.xlsx

Figure 7—source data 2

Summary of changes in downregulated transcriptional regulators based on in situ expression patterns.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig7-data2-v2.xlsx

Download elife-73189-fig7-data2-v2.xlsx

Figure 7—source data 3

Details of probes generated by digestion of cDNAs that were used for in situ hybridisation.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig7-data3-v2.xlsx

Download elife-73189-fig7-data3-v2.xlsx

Figure 7—source data 4

Details of probes generated by PCR that were used for in situ hybridisation.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig7-data4-v2.xlsx

Download elife-73189-fig7-data4-v2.xlsx

Figure 7—source data 5

PCR primer sequences.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig7-data5-v2.xlsx

Download elife-73189-fig7-data5-v2.xlsx

There is a striking correspondence between neural induction by a node graft and normal neural development in terms of the spatiotemporal patterns we observed. In addition to confirming the expression of known neural markers (ERNI, SOX3, SOX2, and SOX1; Figure 7C–F), the predictions from node grafts uncover many novel responses such as CITED4 (ERNI-like), MAFA (SOX3-like), GLI2, and ZNF423 (SOX2-like) (Figure 7G–J). In total, 84/89 transcriptional regulators that are induced by a node graft were detected in neural tissues at some stage (Figure 7 and Figure 7—figure supplements 1–4). Only 5 of the 89 genes (HEY1, RFX3, STOX1, TAF1A, and ZIC2) could not be detected by in situ hybridization – these are expressed at very low FPKM levels according to RNAseq. Most of the transcriptional regulators that are induced by a node graft (79/84) have not previously been associated with neural induction (Figure 7—source data 1). Likewise, a high proportion (65/85; 76%) of the transcripts that are downregulated after a node graft are depleted in neural tissues relative to non-neural territories. Among these are known non-neural markers such as GATA2 and MSX2 and the extraembryonic marker HIF2A (Figure 7K, L and O). Neural/non-neural border markers (such as DLX5, MSX1, and TFAP2C; Figure 7M and Figure 7—figure supplements 5–8) were also downregulated, supporting suggestions that neural induction involves the repression of an early border identity (Trevers et al., 2018). We also uncover several other genes (such as ISX, GRHL3, and VGLL1) that are excluded from neural tissues (Figure 7P–R).

This approach reveals even subtle temporal changes in expression. For example, CITED4 is sharply upregulated after 5 hr of a node graft but is later downregulated (Figure 7G). In the embryo, it is expressed in pre-streak epiblast and prospective neural plate at HH3-4, but CITED4 transcripts are then absent from the neural plate at HH5-6. CDX4 (Figure 7N) and HOXA2 (Figure 7—figure supplement 6J) are initially repressed by a node before their expression levels increase at 12 hr – just as they become enriched in the posterior neural tube from HH7-9. These observations confirm that transcriptional responses to a node graft closely follow the events that occur during development of the embryonic neural plate of the normal embryo.

The GRN describes a temporal hierarchy occurring in single cells during neural development

To follow the expression of multiple genes during normal neural plate development and to compare this to the temporal hierarchy of our GRN, we assessed the expression of network components in individual cells at different stages of normal neural development. scRNAseq was performed on ectoderm dissected from chick embryos at HH4, HH6, HH8, and HH9+. These broad regions included prospective or mature neural plate and neural tube as well as non-neural ectoderm and neural plate border (Figure 8A).

Figure 8 with 1 supplement see all

Download asset Open asset

The non-linear dimensionality reduction method UMAP was then used for unbiased cell clustering in low dimensionality space (Figure 8B). Cells collected from the two earliest developmental stages (HH4 and HH6) each form a distinct group, whereas cells from later stages (HH8 and HH9+) cluster primarily according to identifiable cell populations (Figure 8C). Cell identities were assigned to the 15 clusters using well-established markers (Figure 8—figure supplement 1A). This revealed groups corresponding to placodal (cluster 8), neural crest (clusters 10, 11), and neural tube identities (clusters 4, 6, 7, 9, 13, 14) within the HH8 and HH9+ cell populations. At HH6, neural plate cells expressing LIN28B, SETD2, and SOX2 (cluster 0) are distinct from non-neural plate progenitors (cluster 3), which express GATA2, DLX5, and TFAP2A. At HH4, prospective neural cells express MAFA, ING5, and LIN28B (clusters 1, 2, and 5) whereas cluster 12 expresses the node marker ADMP.

To identify groups of genes that are co-expressed across the entire dataset, gene module analysis was performed using Antler (Delile et al., 2019). Cell state-specific modules were then selected by applying the criterion that at least 50% of the genes in a module must be differentially expressed (log fold change (FC) >0.25, FDR <0.001) in at least one cell cluster relative to the rest of the dataset. This defines 17 distinct modules (Figure 8—figure supplement 1B), which support our cell type classifications. Module-3 shows co-expression of PAX2, WNT4, and NKX6-2 and is upregulated in prospective midbrain clusters 9 (at HH8) and 6 (at HH9+). Module-10 is defined by co-expression of PAX7, SOX8, SOX10, and TFAP2B and is associated with neural crest (clusters 10, 11). Module-25 has low levels of expression of GATA2 and ID3 (clusters 1, 2, 5, 12) at HH4, which are then further downregulated at later stages, except in non-neural cells at HH6 (cluster 3) and placodal cells at HH8 (cluster 8).

To explore the temporal hierarchy of GRN components, we assessed their expression in prospective neural plate and neural tube cell populations (clusters 0, 1, 2, 4, 5, 6, 7, 9, 13, 14). Within the ‘neural’ population, principal component (PC) analysis (PCA) revealed that PC1 inversely correlates with developmental stage (Figure 8—figure supplement 1C), therefore cells were ranked according to their position along PC1 to explore temporal expression dynamics. To visualize the expression of GRN components, genes were grouped according to the time point at which they are first induced by a node graft; their average (normalized and scaled) expression was then plotted (Figure 8D). This analysis reveals that components upregulated within 1–5 hr of a node graft are initially highly expressed in prospective neural cells and that their expression decreases progressively across PC1. GRN components that are induced after 7–12 hr are expressed at higher levels later, in cells of the maturing neural plate and neural tube. Therefore, GRN components identified from node graft experiments are expressed in a temporal hierarchy that closely replicates normal neural plate development.

Responses to neural induction reveal enhancers that drive neural gene expression during normal development

Do the regulatory elements predicted by the GRN generated from Hensen’s node grafts correspond to those that regulate normal development of the neural plate in the embryo? To define the putative enhancers more precisely, we combined our histone profiling (Figure 2 and Figure 2—figure supplement 1) and ATACseq (Figure 2—figure supplement 2) to identify regions of open/active chromatin from induced and uninduced tissues after 5, 9, and 12 hr after a node graft. These regions were then compared to ChIPseq (for H3K27ac/me3) and ATACseq performed on pre-primitive-streak epiblast and mature neural plate at HH6-7.

We chose six putative enhancers for more detailed exploration. These were identified within the previously described CTCF loci using the overlapping combination of H3K27ac marks and ATACseq reads to predict enhancers that are likely to be active (Figure 9—figure supplements 1–3) based on the indices and expression profiles as described earlier. Each region was cloned into the pTK-EGFP vector (Uchikawa et al., 2003) to drive GFP expression via a minimal promoter; the reporter was then electroporated into the pre-primitive-streak stage (EGKX-XI) epiblast or a broad region including the prospective neural plate at HH3/3+ (Figure 9A) and embryos cultured for 5 hr or to the desired stages.

Figure 9 with 3 supplements see all

Download asset Open asset

Figure 9—source data 1

Sequences of primers used for cloning enhancers.

https://cdn.elifesciences.org/articles/73189/elife-73189-fig9-data1-v2.xlsx

Download elife-73189-fig9-data1-v2.xlsx

Each enhancer drives neural-specific GFP expression in a spatiotemporal pattern strikingly similar to the gene that lies in its proximity (Figure 9B–O and Figure 9—figure supplements 1–3). For example, SOX11_enh1 and SETD2_enh1 are enriched for H3K27ac in induced tissues after 5 hr (Figure 9—figure supplement 1A–B). Both are active in the epiblast at EGKXII-XIII and prospective neural plate at HH3-4, when SOX11 and SETD2 are expressed (Figure 9—figure supplement 1 C-AY). SOX2_enh3 is ‘poised’ after 5 hr of neural induction, but acquires an active signature at 9 hr when it gains H3K27ac and ATAC reads (Figure 9—figure supplement 2A). As expected, this enhancer is not active in the prospective neural plate at HH4 (Figure 9—figure supplement 2C–H); GFP positive cells are detected later when SOX2 is expressed in the neural plate proper at HH5-6 (Figure 9—figure supplement 2I–N). Likewise, GLI2_enh1 and SIX3_enh2 also become active after 9 hr; GLI2_enh1 is detected throughout the neural plate, whereas SIX3_enh2 is restricted to the forming forebrain at HH6 (Figure 9—figure supplement 3A–Z). These five enhancers remain active in the neural tube in correspondence with H3K27ac marks and ATAC reads that persist after 12 hr of a node graft.

CDX1_enh1 is also active in the caudal neural tube at HH7-9, when it also weakly labels some presomitic mesoderm cells. Its activity is detected at HH5-6, but only in a few cells of the caudal lateral epiblast adjacent to Hensen’s node and not in the neural plate (Figure 9—figure supplement 2U -AL). This sparse activity overlaps closely with the discrete region from which neuro-mesodermal precursors (NMPs), which are said to express CDX1 (Gouti et al., 2017), are thought to contribute to both neural and somite tissues (Brown and Storey, 2000; Tzouanacou et al., 2009). Since the CDX1 gene is expressed earlier (HH3-4), but only in mesoderm and primitive streak, the activity of this enhancer is likely to be associated with NMPs.

Therefore, the domain of reporter activity for all six predicted enhancers corresponds to that of the normal expression of the associated gene. Furthermore, the stage at which the reporter first becomes active closely matches the stage at which the gene starts to be expressed during normal development.

Finally, to validate predictions from the network, we tested the interaction between MYCN (one of the core network genes; Figure 6A), which is upregulated from 3 hr, and its predicted target ZNF423, which is upregulated at 5 hr onwards (Figure 9P). Electroporation of a morpholino against MYCN into the prospective neural plate at HH3/3+ (Figure 9Q) knocks down the expression of ZNF423 (Figure 9R–S; 5/10 embryos), which is usually expressed later in the neural plate and neural tube (Figure 7J), compared to a control morpholino (Figure 9T–U; 1/10 embryos).

Taken together, these comparisons between node grafts and normal embryos strongly suggest that histone profiles that change in response to a node graft (upon which our network is based) can predict regulatory elements relevant to neural development. This establishes our GRN as a robust tool to describe the genetic hierarchy associated with the acquisition of neural identity.

In addition to the BioTapestry GRN and other associated data described above, the results of this project, including data from node grafts in time course and from normal neural plate development, are made available on the UCSC browser at https://genome.ucsc.edu/s/stern_lab/Neural_Induction_2021. This incorporates RNAseq, ChIPseq, and ATACseq for all expressed and repressed genes (not limited to GRN components) and the annotated putative regulatory regions with predicted binding sites for GRN components. Finally, to help to predict potential conservation of the regulatory sites across species, we performed DREiVe analysis (Sosinsky et al., 2007; Khan et al., 2013; Streit et al., 2013) for domains within a 500 kb window around each of the genes represented in the GRN, using human (hg38), mouse (mm10), rat (rn6), golden eagle (aquChr2), zebrafish (danRer10), and chicken (galGal5). Blocks of sequence with conserved motifs are included as a track in the browser.

Although a graft of Hensen’s node is able to induce a neural plate, neural tube and their subsequent patterning into forebrain, midbrain, hindbrain, and spinal cord, as well as neural crest and placodes from the host epiblast of the area opaca, this is probably not how the process of central nervous system (CNS) development proceeds in the normal embryo. For a start, a large proportion of the cells that will form part of the neural plate (corresponding to the most cranial regions, forebrain, midbrain, and most of the hindbrain), never appear to be in close proximity to the node or anterior primitive streak. This implies that signals from tissues other than the node or streak are likely to play a role, especially at the earliest stages of neural development. We have shown that neural plate development is likely to start very early, before primitive streak formation, and probably under the initial influence of the hypoblast, involving FGF signals (Streit et al., 2000; Trevers et al., 2018). Therefore, although the node can induce an entire CNS including the earliest steps (as it also expresses FGF8), this is not what happens in normal development.

A second important issue arises by considering the formal definition of induction: ‘an interaction between one (inducing) tissue and another (responding) tissue, as a result of which the responding tissue undergoes a change in its direction of differentiation’ (Gurdon, 1987). The ‘change in the direction of differentiation’, in more contemporary terms, might be described as a ‘change of fate’ or ‘change in tissue identity’. By implication, to conclude that an inductive interaction takes place, we first need to define an inducing tissue (source of signals), the outcome of the induction (the final identity/fate resulting from the interaction) and importantly, a responding tissue whose normal fate, in the absence of induction, should not be the same as the outcome – otherwise there would be no ‘change’. This means that it is not possible to determine whether the prospective neural plate of the normal embryo develops as a result of one or more inductive interactions, especially since we cannot determine the earliest stage at which the first inducing signals are presented to that responding tissue.

For these reasons, the time course initiated by a transplanted tissue that is able to induce a complete nervous system in a responding tissue whose fate is not to form a neural plate is an invaluable tool to study the acquisition of neural plate identity. Here, we have taken advantage of this approach by first constructing a GRN with a fine time course in a competent responding tissue after such a transplant, defining the hierarchy and dynamics of expression of more than 150 transcription factors and other markers, along with changes in chromatin marks that accompany this. Then, we addressed the question of whether this hierarchy of events is in any way comparable to what occurs during normal neural development, using three complementary approaches:

First, using in situ hybridization for the majority of genes contained in the GRN, noting the time and location of their initial expression, we assessed whether markers upregulated by a node graft are expressed in the normal neural plate, and whether those repressed by a graft become excluded from the neural plate, and from what stage. This was then compared with the order of events after a node graft. We found a very strong correspondence: genes whose expression is induced by a node graft within 3 hr tend to be expressed in the pre-primitive-streak stage embryo (around stage EGK XIII), those induced at 5–7 hr in the central area pellucida anterior to and around the node at stages 4–6, and those induced at 9–12 hr are expressed in the elevating neural plate from stages 7–9, confirming previous observations made with just a few genes (Figure 1 and Figure 1—figure supplement 1).

Second, we generated reporters for a subset of the predicted enhancers from the GRN, containing a minimal promoter (TK), the putative enhancer, and the coding sequence of GFP. This was co-electroporated into a wide region of the early embryo along with a vector encoding a fluorescent protein under the control of a strong ubiquitous promoter (DS-Red-Express) to mark the entire territory that had been electroporated. We then assessed the area within the electroporated territory (red from DS-Red) where the enhancer was active (green from GFP). In all six cases examined, the domain of expression corresponded to that of the normal expression of the associated gene, and the stage at which the reporter first became active closely matched the stage at which the gene starts to be expressed.

Finally, we compared gene expression of the neural plate/tube and the neural plate border in the normal embryo at single-cell level at four stages of development (primitive streak stage HH4, head fold stage HH6, four-somite stage HH8, and nine-somite stage HH9+; the time interval between each pair of these stages of normal development is approximately 2–3 hr). Dimensionality reduction and unsupervised clustering revealed clusters of cells representing neural plate and its different anterior-posterior parts, non-neural ectoderm, neural crest, and placodal cell states. Modelling transcriptional dynamics of components of the GRN across inferred pseudo-time revealed a very close correspondence between normal neural plate development and the sequence of transcriptional changes which take place following a node graft.

Taken together, these three approaches all strongly indicate that the sequence of events that characterize the response of area opaca epiblast cells to a graft of the organizer closely follows the dynamics of gene expression in the normal neural plate.

Our analysis reveals thousands of sites across the genome containing marks that change with the progression of neural induction. Each locus contains genes whose expression changes during this cascade, and their associated putative regulatory elements. As previously described for other systems such as embryonic stem (ES) cells in culture, regulatory elements may be decorated simultaneously by marks associated with activation (e.g. H3K27ac) and with repression (e.g. H3K27me3); these have been called ‘poised’ or ‘bivalent’ enhancers (Cui et al., 2009; Heintzman et al., 2009; Xu et al., 2009; Creyghton et al., 2010) and have been interpreted to indicate that the site is in transition from one state to the other, or perhaps able to go either way. Our fine time course reinforces this concept; in many cases we observe that cells pass from one state (e.g. H3K27me3 without acetylation) to the opposite, sometimes passing through a ‘poised’ state at the intermediate time point. This could suggest that these changes can occur sequentially in a cell. However, other interpretations are also possible. One possibility is that the apparent contradictory marks are due to different modifications of the two alleles of the locus in question, which could also be a transitional state. Another possibility is that this progression may reflect cells within a tissue that do not undergo the change synchronously, but some cells start earlier than their neighbours, the tissue gradually becoming more homogeneous. This sort of progression is sometimes seen by in situ hybridization or time-lapse movies of reporter activity (e.g. to visualize the activity of enhancer N2 associated with the SOX2 gene [Uchikawa et al., 2003; Papanayotou et al., 2008]). Likewise, in the present paper, SOX11 expression begins as a mosaic (salt-and-pepper) in the early pre-streak epiblast before becoming more homogeneous, and the activity of the SOX11_enh1 (Figure 9—figure supplement 1C–H) shows a similar pattern at the same stage. Most likely these apparently ‘poised’ or ‘bivalent’ sites reflect a combination of these reasons.

Among the many regulatory elements predicted by our analysis are some that have been described previously (Uchikawa et al., 2003; Uchikawa et al., 2011; Okamoto et al., 2015; Iida et al., 2020). Their analysis was based on conserved stretches of non-coding sequence, followed by testing the activity of fragments in a vector containing a heterologous minimal promoter. In some cases, there are differences in the length of the sequence considered to act as an enhancer by both studies. For example, our SOX2_enh3 partially overlaps with, and its behaviour upon electroporation mimics the activity of, the previously identified SOX2 D1 enhancer (Okamoto et al., 2015; Iida et al., 2020). Therefore, our ChIPseq analysis can verify known enhancers as well as identifying many new candidates.

Calculations of ‘outdegrees’ and ‘betweenness centralities’ of the components in the network across seven time points (0, 1, 3, 5, 7, 9, and 12 hr) reveals a relatively small subset of transcriptional regulators with high values in both parameters (Figure 6). These are defined as ‘core genes’ (or ‘network hubs’) in the network and may play more important roles during the process of neural induction. We tested whether one of these, MYCN, is required to regulate its predicted target, ZNF423 at the appropriate time point (Figure 9R–S) using morpholino electroporation. The knockdown reduced expression of the target gene in a significant proportion of embryos, confirming MYCN as a regulator of ZNF423 expression. However, it is worth considering that in the case of other enhancers such as the N2 element of SOX2, it appears that no single TF binding to the element is individually essential – rather, several TFs need to bind but the precise combination may not be so critical (Uchikawa et al., 2003; Uchikawa et al., 2011; Okamoto et al., 2015; Iida et al., 2020). This may be a general feature of such developmentally regulated elements, therefore testing whether a factor is ‘essential’ can give a misleading view of gene regulation in the biological context being studied.

The GRN has great complexity over a relatively short period of time, with many transcription factors changing their expression and thousands of predicted interactions between them undergoing concomitant changes. We initially expected that a detailed view of the responses to signals from the organizer might somehow highlight specific changes, or states, that could be associated with some of the key concepts of developmental biology, such as ‘induction’ and ‘commitment’ (Slack, 1991; Stern, 2004) to a neural plate identity. These concepts have classically been defined based on specific embryological manipulations but it has hitherto not been possible to identify specific molecular signatures associated with these transitions. Could we use the present GRN to find objective molecular criteria that correspond to such changes in developmental properties?

Based on timed organizer grafts followed by their removal, a critical time was established at 12–13 hr after a graft (Gallera and Ivanov, 1964; Gallera, 1965; Gallera, 1971b; Gallera, 1971a), using similar conditions to those of our node grafts here. At this time point, removal of the node graft does not impair the subsequent development of the neural plate. This represents the endpoint of our GRN; at this time SOX2 is robustly expressed, and SOX1 is just starting to be expressed, in the responding tissue (Pevny et al., 1998). Although this functional assay does not test for ‘irreversibility’ of the cell fate transition, the finding that only after this point it can proceed autonomously without further signals from the organizer suggests that this time point might correspond to at least some degree of ‘commitment’ (or ‘determination’) to a neural plate identity.

Another earlier study has defined an initial state (named the ‘common state’) reflecting the earliest responses (1–3 hr) to a node graft into the area opaca (Trevers et al., 2018). A similar state is elicited by a graft of lateral head mesoderm (which induces placodes) and by a graft of the node (which induces the CNS), and the genes comprising this ‘common state’ have some similarity to those expressed by ES cells. In normal embryos, these genes are expressed in a broad region of the embryonic epiblast; moreover, a graft of the hypoblast into the area opaca also induces a similar set of genes. These observations suggested that the responses to inducing tissues start by a ‘rewinding’ (or reprogramming) cells to a state similar to that of the very early embryonic epiblast, from where cells can then be redirected to their new fate: neural plate, placode, neural crest, epidermis, and perhaps others. A similar state can also be elicited by exposure of area opaca to a source of FGF8 (a factor normally secreted by both the hypoblast and the node, but lost by the node at the time when its inducing ability starts to decline; Storey et al., 1992; Streit et al., 2000).

Functional studies have established that competent ectoderm cannot respond to BMP inhibitors like Chordin unless it has previously been exposed to a graft of the organizer for at least 5 hr (Streit et al., 1998). This led to the identification of a few genes whose expression is associated with this transition (Streit et al., 2000; Sheng et al., 2003; Gibson et al., 2011; Pinho et al., 2011; Papanayotou et al., 2013), but the present experiments greatly expand and enrich this by placing genes in the context of a comprehensive network of gene interactions. Future experiments can take advantage of this information to determine the mechanisms that cause epiblast cells to become sensitive to BMP inhibition at this time point, and the upstream signals responsible.

To define possible ‘key’ steps, if they exist, we have attempted to identify ‘core’ factors that appear to represent important hubs in the GRN and their downstream targets and interconnections (Figure 6). This revealed a small number of genes that appear to occupy nodal, highly connected, positions in the network, and include, at successive time points, ESRRG, GRHL2, and TFAP2A in epiblast not exposed to a graft (0 hr) but then almost immediately downregulated, then TFAP2C, SOX4, and BLIMP1 just 1 hr after a node graft. Indeed, BLMP1 (PRDX1) has recently been shown to be required transiently at an early stage for neural plate, crest, and placode development, and its subsequent downregulation is also essential for further differentiation (Prajapati et al., 2019). Other later putative core genes include SOX3, TCF12, and MYCN. Future experiments could address whether these do indeed represent key steps accompanying neural induction. For the time being, we feel that functional definitions of the developmental concepts will remain key to defining the properties of developing tissues, especially given the complexity of the changes in gene expression accompanying neural induction.

The gene signature that defines the ‘common state’ is induced by a graft of Hensen’s node after approximately 3 hr. The same genes are expressed in the epiblast of the normal embryo before primitive streak formation (stage EGK XIII), long before the node has formed; they can be induced by a graft of the hypoblast, as well as by the node. This suggests that in the normal embryo, it is the hypoblast rather than the node that initiates this process (Foley et al., 2000; Streit et al., 2000; Trevers et al., 2018). However, a node graft can mimic both these initial steps and those that follow. Therefore, although a node graft is able to generate a complete neural plate when grafted into a region of competent epiblast that normally does not give rise to this tissue, neural plate development in the normal embryo is likely to be directed by signals emanating from more than one tissue. A large part of the prospective CNS (forebrain, midbrain, and anterior hindbrain) is never close to the node but extends far anteriorly from the primitive streak. Indeed it has been suggested that the most anterior nervous system may be induced by signals other than the node, as well as earlier in development, and that the node and its derivatives may be involved mainly in the induction of the more caudal parts of the CNS (posterior hindbrain and spinal cord) (Beddington, 1994; Thomas and Beddington, 1996; Foley et al., 2000; Streit et al., 2000; Metzis et al., 2018). Armed with the GRN presented here, it will be particularly interesting in future to assess which sub-modules (‘kernels’) of the GRN are regulated by which secreted factors expressed in the node and whether and when similar factors are also normally expressed in other tissues that display partial or full neural induction activity.

A growing literature makes use of protocols to direct cultured mouse ES cells or human induced pluripotent stem cells into neurons, and calls this process ‘neural induction’. After maintaining the stem cell population in a medium containing FGF to sustain proliferation and pluripotency, a common current protocol involves withdrawal of the FGF, followed by activation of the hedgehog pathway (Gouti et al., 2014; Metzis et al., 2018). Other protocols in the literature include BMP inhibition, manipulation of the Wnt pathway, and/or retinoid signalling, in various combinations. Almost all of these protocols need at least 5 days, and sometimes as long as 14 days, for the culture to generate neuronal cells. This is a much longer period than we see in vivo, which raises the question of to what extent they represent equivalent processes. Our network should help to compare them and to assess the degree of similarity as well as identify any differences. We have already shown that the initial (‘common’) state in response to induction (the first 3 hr or so in our assays) generates a state comparable to ES cells, suggesting that the remaining part of the differentiation protocols correspond to the period after this.

In conclusion, our study provides the first GRN to represent the changes that accompany the process of neural induction after grafting an organizer to an ectopic site, and how it relates to normal neural plate development, in fine time course. We provide a comprehensive resource (deposited in the UCSC genome browser) that includes not only the network components (transcriptional regulators with predicted cross-regulatory interactions), but also a genome-wide view of all transcriptional changes, H3K27 acetylation/methylation, and chromatin conformation (ATACseq) during this process, from the time cells are initially exposed to neural inducing signals to the time at which definitive neural plate markers start to be expressed. We envisage that the GRN and the accompanying resource will be invaluable tools to study other important aspects of early neural development, such as the signalling inputs responsible for these changes, the molecular basis of the competence of cells to respond to neural inducing signals at different locations and stages, and aspects of neural patterning leading to regional diversification of the CNS. We also envisage that this resource will be useful to understand the process of neural induction not only in the chick but also in other vertebrate species, as well as in protocols for neural differentiation in ES cells in vitro. Indeed, we have noticed that out of 84 upregulated genes included in this network, the majority (74/84) are expressed in the prospective mouse neural plate (Figure 7—source data 1). These 84 genes include 69 not previously associated with neural induction in any species, of which the majority (56/69) are expressed in the future mouse neural plate. These observations suggest a high degree of evolutionary conservation. To aid in cross-species comparisons and possible studies on evolutionary changes that have occurred in neural induction, the browser resource also includes the results of computational prediction of conservation of the key GRN loci among several vertebrate species, using DREiVe, a computational tool designed to identify conserved sequence motifs and patterns rather than precise sequence conservation (Sosinsky et al., 2007; Khan et al., 2013; Streit et al., 2013).

Fertile hens’ eggs (Brown Bovan Gold; Henry Stewart & Co., UK), were incubated at 38°C in a humidified chamber to the desired stages. Embryos were staged according to Hamburger and Hamilton, 1951, in Arabic numerals or Eyal-Giladi and Kochav, 1976, in Roman numerals for pre-primitive streak (pre-streak) stages. Embryos were harvested and dissected in Pannett-Compton saline (Pannett and Compton, 1924) or 1× PBS, before being fixed in 4% paraformaldehyde (PFA) overnight at 4°C or prepared for ex ovo culture. All experiments were conducted on chicken embryos younger than 12 days of development and therefore were not regulated by the Animals (Scientific Procedures) Act 1986.

Chicken embryos were cultured ex ovo using a modified New culture method (New, 1955; Stern and Ireland, 1981) as previously described (Voiculescu et al., 2008). In brief, eggs were opened and the thin albumin was collected. Intact yolks were transferred to a dish containing Pannett-Compton saline. The vitelline membrane was cut at the equator while keeping the embryo central. The membrane was then slowly peeled away from the underlying yolk to keep the embryo attached. Membranes were placed on a watch glass with the embryo orientated ventral side up. A glass ring was positioned over the embryo and the edges of the membrane wrapped over the ring. This assembly was then lifted out of the dish and adjusted under a dissecting microscope while keeping the embryo submerged in a small pool of saline. The vitelline membrane was gently pulled taut around the glass ring before the excess was trimmed. Excess yolk was cleared from the membrane and embryo by gently pipetting a stream of saline and then replaced with fresh liquid, leaving the embryo on an optically clear membrane. At this point, embryos can be cultured as they are or with the addition of a node graft which were transferred to the ring and positioned while the embryo is submerged. To complete the culture, saline was removed from around and within the ring without disturbing the embryo and grafts. The dry ring was then transferred to a 35 mm Petri dish containing a shallow pool of thin albumin. The edges of the ring were pressed down to prevent it from floating, leaving the embryo supported by a shallow bubble of albumin beneath the membrane. Completed cultures were incubated in a humidified chamber at 38°C for the desired length of time. After culture, embryos were fixed on the membrane with 4% PFA or submerged with ice-cold saline to allow tissue dissection.

Neural induction assays were performed as previously described (Stern, 2008; Streit and Stern, 2008). For Hensen’s node grafts, chick donors and hosts at Hamburger-Hamilton (HH) 3+/4- were used and New cultures incubated for 1, 3, 5, 7, 9, or 12 hr at 38°C in a humidified chamber. One or two nodes were grafted per embryo, placed contralaterally within the inner third of the area opaca, at or above the level of the host node. This region is competent to respond to neural inducing signals but only contributes to the extra-embryonic membranes and not the embryo proper (Streit et al., 1997). Nodes were grafted with their endodermal surface in contact with the host epiblast.

For RNAseq, HH4- chick nodes were grafted to the area opaca of HH4- chick hosts. A single node was grafted to the left or right side per host, and embryos were cultured for 5, 9, or 12 hr. Node grafts were removed by submerging the embryo in saline while still attached to the membrane. Fine syringe needles (27 G or 30 G, BD Microlance) were used gently to lift the grafted tissue away from the underlying epiblast. Where grafts were firmly attached, 0.12% trypsin dissolved in saline was gently pipetted over the graft site to assist in removal. After dissection, trypsin activity was neutralized by treating embryos briefly with heat-inactivated goat serum (Stern, 1993), which was removed by washing with fresh saline. The induced epiblast directly beneath the graft, which appeared greyish and thickened, was dissected using syringe needles and mounted insect pins. Uninduced epiblast tissue from same position on the contralateral side was also dissected. Tissue samples per condition were pooled, frozen on dry ice, and stored at –80°C. At each time point (5, 9, or 12 hr) a total of 50 induced and contralateral uninduced pieces of tissue were collected. A further 50 pieces of uninduced tissue were collected from HH4- embryos that had not been cultured, representing a 0 hr control. Tissue samples for each condition were then pooled and lysed in 1 mL TRIzol (Invitrogen) for RNA extraction.

RNA sequencing was conducted by Edinburgh Genomics (formerly ARK Genomics; Roslin Institute, University of Edinburgh, UK). Total RNA was extracted and the quality was assessed using the Agilent 2100 Bioanalyzer – all samples registered an RIN value between 9.0 and 10.0. Libraries were constructed using the Illumina TruSeq mRNA library preparation kit. In brief, mRNAs were purified using oligo-dT conjugated magnetic beads before being chemically fragmented to on average 180–200 bp. Fragments were then transcribed using short random primers and reverse transcriptase to produce single-stranded cDNA, from which double-stranded cDNA was generated by DNA polymerase I and RNase-H. After synthesis, cDNA was blunt-ended and a single A-base added to the 3’ end. Sequencing adapters were ligated via a T-base overhang at their 3’ and these libraries were then purified to remove unincorporated adapters before being enriched by 10 cycles of PCR. Library quality was checked by electrophoresis and quantified by qPCR. The seven RNA libraries were sequenced over two lanes via 100-cycle, paired-end sequencing using the Illumina HiSeq 2000 system.

Antisense riboprobes were generated from cDNA plasmid templates or Chick Expressed Sequence Tag (ChEST) clones (Source Bioscience). Templates were linearized by restriction digest (Figure 7—source data 3) or by PCR (Figure 7—source data 4) using M13 forward (5’-GTAAAACGACGGCCAGT-3’) and M13 reverse (5’-GCGGATAACAATTTCACACAGG-3’) primers (Figure 7—source data 5). The products were transcribed using SP6, T7, or T3 (Promega) and DIG-labelled nucleotides. DNA templates were digested using RNase-free DNase and the RNA probe was precipitated and resuspended in water. Working probes were diluted to 1 μg/mL in Hybridization (HYB) buffer and stored at –20°C.

Whole-mount in situ hybridization was performed as described previously (Stern, 1998; Streit and Stern, 2001), but omitting pre-adsorption of the anti-DIG-AP antibody (Roche 11633716001, used at 1:5000). Stained embryos were imaged from the dorsal perspective using an Olympus SZH10 Stereomicroscope with an Olympus DF PlanApo 1× objective and an Olympus NFK 3.3× LD 125 photo eyepiece. Images were captured using a QImaging Retiga 2000R Fast 1394 camera and QCapture Pro software as 24-bit, colour TIFF files at 300 dpi with dimensions of 1600×1200 pixels.

A spatiotemporal expression map was generated for genes that are differentially expressed after either 5, 9, or 12 hr of a node graft. These were selected from the galGal3 and galGal4 RNAseq analyses according to the criteria as described in ‘Differential expression analysis’.

Uninduced tissue or induced epiblast that had been exposed to node grafts was dissected as previously described. Four to eight pieces of tissue were collected in triplicate per condition. Tissues were dissected in ice-cold 1× PBS and collected on ice. All excess solution was removed using a fine needle and promptly processed by adding a total of 1 µL of lysis buffer from the RNAqueous-Micro Total RNA Isolation kit (Thermo Fisher). Tubes were immediately snap-frozen on dry ice and stored at –80°C.

NanoString experiments were run on the nCounter Analysis System using a custom codeset and following NanoString guidelines. The codeset consisted of 386 probes (Figure 1—source data 3) belonging to various categories. Transcriptional regulators that are upregulated (coloured red) or downregulated (coloured green) at 5, 9, or 12 hr were selected from RNAseq screening (details in ‘Differential expression analysis‘). Other genes that respond to 5 hr exposure to a node graft (Streit et al., 2000; Sheng et al., 2003; Gibson et al., 2011; Pinho et al., 2011; Papanayotou et al., 2013) were also included. Probes were also designed against housekeeping genes (ACTB, GAPDH), markers of apoptosis and proliferation, transcriptional readouts of FGF, BMP, WNT, retinoic acid, Notch and Hedgehog signalling; epithelial, mesodermal, endodermal, neural plate border, pre-placodal, neural crest and Hensen’s node markers, and standard positive and negative NanoString controls.

Tissue samples were processed using the NanoString master kit and following the manufacturer’s instructions. Probes were hybridized to lysates overnight for 17 hr at 65°C before the reactions were transferred to the NanoString prep-station robot and processed using the ‘high-sensitivity’ program. Probe-target complexes were digitally counted from 600 fields of view on the NanoString Analyzer.

Enhancers were identified based on changing active, inactive, or poised chromatin signatures over 5, 9, and 12 hr of neural induction. These were referenced against the 5, 9, and 12 hr ATACseq data to select putative enhancers with accessible chromatin. The primers used for cloning are detailed in Figure 9—source data 1.

Enhancers were cloned from chick gDNA as previously described (Chen and Streit, 2015). The correct sized fragments were amplified by PCR using PCRBio Ultramix (PCRBiosystems) and ligated into the pTK_JC_EGFP plasmid (Chen and Streit, 2015) using T4 ligase (NEB). Ligation products were transformed into DH5α Escherichia coli and grown on agar plates with Ampicillin selection. Clones were screened by colony PCR using PCRBio Ultramix and pTK_Fwd (5’-GTGCCAGAACATTTCTCTATCG-3’) and pTK_Rev (5’-GTCCAGCTCGACCAGGATG-3’) primers. Positive clones were cultured overnight in 5 mL LB broth plus Ampicillin and the constructs were purified by mini-prep (Qiagen). Enhancer inserts were sequenced using Citrine_Fwd (5’-TGTCCCCAGTGCAAGTGC-3’) and Citrine_Rev (5’-TAGAACTAAAGACATGCAAATATATTT-3’) primers. Plasmids with sequence verified inserts were maxi-cultured in 200 mL, of which 1 mL used to make a 50% glycerol stock and stored at –80°C. Bulk plasmid DNA was extracted using the Endofree Plasmid Maxi-prep kit (Qiagen) following the manufacturer’s instructions and resuspended in 100 µL of sterile water. Further purification was conducted to reduce residual salts which interfere with electroporation efficiency and embryo damage. Samples were spun for 15 min/15,000 rpm/4°C to pellet any particulates. The supernatant was removed to a fresh tube, to which 10 µL of 3 M NaOAc pH 5.2 and 250 µL of 100% EtOH was added to precipitate the plasmid DNA. This was pelleted by centrifugation at 15 min/15,000 rpm/4°C and the pellet was washed with 70% EtOH and spun again for 10 min/15,000 rpm/4°C. The supernatant was removed and the pellet allowed to air-dry completely, before being resuspended in 50 µL of sterile water and the concentration checked by nanodrop and adjusted to 5 mg/mL.

Enhancer plasmids were co-electroporated with pCMV-DsRed-Express-N1 (Clontech) which labels all electroporated cells. Tissues were electroporated using an anode electroporation chamber and wire cathode according to Voiculescu et al., 2008, with minor modifications. Electroporation mix was made fresh by combining the following in water: enhancer_TK_EGFP plasmid (2 mg/mL), pCMV-DsRed-Express-N1 (1 mg/mL), 6% sucrose, 0.02% Brilliant Blue FCF. This was backloaded into a capillary needle with a fine tapered tip. The anode chamber was filled with 1× PBS and the electrodes positioned ≈2 cm apart.

Harvested embryos were positioned directly between the electrodes with the dorsal side (epiblast) upwards, closest to the cathode. Electroporation mix was applied using a mouth aspirator attached to the capillary needle, and by bringing the needle tip in close contact with the embryo. For pre-streak EGKX-XI embryos, electroporation mix was applied broadly across the central epiblast region and immediately electroporated using the following current settings: 3.5–4.0 V, 5 pulses, 50 ms duration and a 100 ms gap. At HH4- stages, electroporation mix was applied broadly across the prospective neural plate and primitive streak and immediately electroporated using 4.5–5 V, 5 pulses, 50 ms duration and 100 ms gap. After electroporation embryos were transferred to vitelline membranes, prepared for New culture and cultured to the desired stages (Voiculescu et al., 2008).

Enhancer activity was imaged using a Leica SPEinv microscope with a 5× objective and LAS-X software. Images were captured at 2048×2048 pixels, speed 400 ms, and 1× zoom with an averaging of 4. EGFP was excited using a 488 nm laser and detected with a 492–543 nm wavelength filter. DsRed was excited using a 532 nm laser and detected with a 547–653 nm wavelength filter. Files were saved as.lif format and composite images were stitched together in Fiji (Schindelin et al., 2012) using the Pairwise Stitching plugin (Preibisch et al., 2009).

Fluorescein-labelled morpholinos at 500 µM (MYCN: TCTTGCTGATCATTCCCGGCATGGC, or Standard control: CCTCTTACCTCAGTTACAATTTATA) were co-electroporated with pCMV-DsRed-Express-N1 (1 mg/mL) as a carrier. Constructs were targeted to the lateral prospective neural plate of embryos at HH3/3+ and electroporated using the following conditions: 4.0 V, 5 pulses, 50 ms duration and 100 ms gap. Embryos were then cultured overnight for 12 hr to reach stages HH5-8.

In situ hybridization was performed (as described above) to reveal ZNF423 expression using an anti-DIG-AP antibody and NBT/BCIP. Embryos were then post-fixed overnight in 4% PFA. They were then washed extensively in PBST and incubated at 70°C for 1 hr. Further PBST washes were performed to remove residual traces of NBT-BCIP. Immunohistochemistry was then performed to detect fluorescein-labelled morpholinos using anti-FITC-AP (Roche 11426338910, used at 1:5000) and BCIP alone.

In brief, tissues were dissected in ice cold 1× PBS and transferred to low binding PCR tubes containing 50–100 µL of 1× Protease Inhibitor Cocktail in 1× PBS (cOmplete mini EDTA-free; Roche) on ice. Where necessary, 0.12% trypsin or 1× PBS + 200 mM EDTA was used to aid dissection. Treated tissues were rinsed with PBS to remove residual trypsin. Tissues were collected in batches every 30 min before tubes were spun at 100 × g/3 min/4°C to gently pellet the tissue. Tubes were snap-frozen in liquid nitrogen before storage at –80°C for up to 6 months. For pre-streak EGK XII-XIII embryos, the hypoblast was removed and central epiblast was collected from six embryos in triplicate. After removal of the underlying mesoderm and endoderm, the neural plate at HH6-7 was harvested from 20 embryos in triplicate. For 0, 5, 9, and 12 hr node induced and contralateral uninduced area opaca, 14–30 pieces of tissue were collected in triplicate. A further 60 pieces of 9 and 12 hr uninduced tissue were collected and processed separately after it was discovered that 30 pieces of uninduced tissue at these time points did not generate enough material for ChIP.

ChIP was performed by micrococcal nuclease (MNase) digestion to fragment the chromatin following a previous protocol (Brind’Amour et al., 2015), followed by immunoprecipitation according to the Low Cell ChIPseq kit (Active Motif, 53084) with some modifications. Corresponding induced and uninduced samples were processed side-by-side, each as three reactions: H3K27ac, H3K27me3, and input. In brief, 200 µL of protein G agarose beads were added to 2×1.5 mL low binding tubes and spun for 3 min at 1250 × g and 4°C. The supernatant was removed and the beads were washed with 690 µL of TE buffer at pH 8.0. These tubes were spun again for 3 min at 1250 × g and 4°C and the supernatant removed before bead blocking reactions were set up. To one tube, 200 µL of TE buffer pH 8.0, 20 µL of blocking reagent AM1 and 20 µL of BSA were added. This tube of ‘pre-clearing’ beads was incubated for 3–6 hr on a 4°C rotator. To the second tube, 180 µL of TE buffer pH 8.0, 20 µL of blocking reagent AM1, 20 µL of BSA, and 20 µL of Blocker were added. This tube of ‘IP’ reaction beads was incubated overnight on a 4°C rotator.

While beads were blocking, tissues were prepared for nuclei isolation and MNase digestion. A pair of samples (induced and contralateral uninduced) were thawed on ice for 5–10 min and the tissues pooled into one tube per condition using low binding tips. Pooled samples were spun at 100 × g/3 min/4°C and the excess supernatant was removed without disturbing the tissues, to leave ≈5 µL remaining. To each tube, 20 µL of nuclei lysis buffer was added, containing 16 µL of nuclei isolation buffer (Sigma N3408), 2 µL 1× protease inhibitor (Roche cOmplete mini, EDTA-free) and 2 µL of 1% IGEPAL CA-630. Tubes were gently tapped to mix and incubated on ice for 5 min. Nuclei were dissociated gently by pipetting 15 or 30 times without making bubbles, then incubating samples on ice for 15 min before repeating these steps four to five times as listed in Figure 2—source data 1.

Chromatin was then fragmented using MNase. This was carried out in two steps to fragment the easily digestible chromatin and then re-digest the more compact chromatin in a second longer digestion. A 3× master mix containing 169.5 µL of sterile H₂O, 30 µL of 10× MNase buffer (NEB), 30 µL of 50% PEG 6000, and 6 µL of 100 mM DTT (made fresh) was prepared. In a separate tube, 3 µL of 10× MNase buffer (NEB M0247) was diluted in 25 µL of H₂O, before adding 2 µL of MNase, to give a 1:15 dilution of MNase enzyme. Of this, 4.5 µL of MNase 1:15 was added to the 3× master mix, to give a total volume of 240 µL. Then, 50 µL of 3× MNase master mix was added to each tube containing 25 µL of isolated nuclei. The mixture was pipetted 15 times to mix without making bubbles, while tubes were kept on ice. This mixture was transferred to a fresh 0.2 mL low binding PCR tube at room temperature and the digestion incubated at room temperature for 2.5 min. This tube was returned to ice and 6 µL of 200 mM EDTA added. Tubes were spun at 500 × g and 4°C for 5 min to pellet larger undigested chromatin and the supernatant (containing smaller, readily digestible chromatin) was transferred to a fresh 1.5 mL low binding microfuge tube. A further 30 µL of the MNase master mix was added to the remaining chromatin pellets on ice and resuspended by pipetting 60 times, before tubes were incubated for 6 min at room temperature. To break up remaining nuclear debris, 4 µL of 200 mM EDTA and 4 µL of 1% TX-100 +1 % sodium deoxycholate were added and the mixture pipetted 30 times. All 38 µL of this second digestion was collated with the first supernatant and kept on ice (total volume ≈110 µL).

Once the ‘pre-clearing’ beads were ready (after 3–6 hr) they were spun for 3 min/1250 × g/4°C. The supernatant was removed and the beads washed with 665 µL of ChIP Buffer (from Active Motif kit) by inverting the tube several times. Beads were spun again for 3 min/1250 × g/4°C and supernatant removed again, before 232 µL of ChIP buffer was added. An even suspension of 100 µL of beads was added to each tube of MNase digested chromatin (≈110 µL). Then, 10 µL of proteinase inhibitor cocktail, 10 µL PMSF (from the Active Motif kit), and a further 290 µL of ChIP buffer were added, giving a total volume of 520 µL. Tubes were then placed on a rotator at 4°C for 3 hr to pull down chromatin fragments that bind non-specifically to the protein G agarose beads.

After ‘pre-clearing’, samples were spun at 3500 rpm and 4°C for 3 min and the supernatant was divided into three separate 0.5 mL low binding PCR tubes: 200 µL each into H3K27ac and H3K27me3 tubes and the rest into the input tube (≈60 µL). To the experimental tubes, 2 µL of H3K27ac (Active Motif 39133, lot# 31814008) or H3K27me3 (Active Motif 39155, lot# 31814017) antibody was added, and these were incubated overnight at 4°C. To the input, 80 µL of low EDTA TE buffer was added and stored overnight at 4°C.

The following day, ChIP reaction tubes were spun briefly at 4°C to collect the solution. ‘IP’ reaction beads were removed from the rotator and mixed to an even suspension by pipetting with a wide bore tip. Keeping the beads well mixed, 50 µL of bead slurry was added to each immunoprecipitation and the tubes returned to the 4°C rotator for 3–4 hr to pull down ChIPed chromatin. Next, samples were washed according to Section G of the Active Motif Low Cell ChIPseq manual and eluted with 100 µL of Elution Buffer AM4.

Chromatin was de-crosslinked by incubating ChIP and input samples at 65°C for 2 hr. Then, 125 µL of phenol:chloroform:isoamyl alcohol 25:24:1 and 64 µL of chloroform:isoamyl alcohol mixture was added and samples shaken vigorously for 15 s before centrifuging at 16,000 × g/15 min/4°C. The top aqueous layer was transferred to a fresh low binding tube. An additional 100 µL of nuclease free water was added to the original tube, which was shaken and centrifuged again. Then the aqueous top layer was removed and added to the first. To this, 2 µL of blue glycogen and 20 µL of 3 M NaOAc pH 5.2 were added and samples flicked to mix. Finally, 660 µL 100% ethanol was added, and tubes flicked to mix well before precipitating overnight at –20°C. The next day, the chromatin was centrifuged at 16,000 × g/30 min/4°C. Supernatant was removed and the pellet washed with 900 µL of cold 70% EtOH. The pellet was spun at 16,000 × g/30 min/4°C and washed with 70% EtOH again, a total of four times. Then, chromatin pellets were air-dried completely for 10–15 min and resuspended in 20 µL of TE low EDTA buffer overnight at 4°C, before being stored at –20°C until library preparation.

Libraries were prepared using the Next Gen DNA Library kit (Active Motif 53216) and Next Gen Indexing kit (Active Motif 53264) but using ProNex size-selective beads (Promega) to enrich for mono- and di-nucleosomes (≈200–500 bp). Each sample of ChIPed chromatin was thawed on ice, to which 20 µL of low EDTA TE Buffer was added. For each library, ‘Repair 1’ reaction mix was prepared according to the Active Motif manual by combining 13 µL of low EDTA TE, 6 µL of Buffer W1, and 1 µL of Enzyme W2. This was added to the 40 µL of chromatin and mixed by pipetting. Reactions were incubated for 10 min at 37°C in a thermocycler with the lid open. To clean up the reaction, 180 µL (3× volume) of evenly suspended ProNex beads were added and incubated for 10 min at room temperature before placing on a magnetic rack for 5 min. The supernatant was removed and beads were washed three times with 100 µL of wash buffer, for 30–60 s each. The final wash buffer was removed and the beads were air-dried for 5 min. Next, Repair II reactions were prepared by combining 30 µL of low EDTA TE, 5 µL of Buffer G1, 13 µL of Reagent G2, 1 µL of Enzyme G3, and 1 µL of Enzyme G4. This was added to dry beads and pipetted to mix. Tubes were incubated on a thermocycler for 20 min at 20°C with lid open. Next, 90 µL (1.8× volume) of PEG NaCl was added to each reaction and mixed by pipetting 15 times. Samples were incubated at room temperature for 10 min and then on the magnetic rack for 5 min and beads washed and air-dried as before. Next, Ligation I mix was prepared by combining 20 µL of low EDTA TE, 3 µL of Buffer Y1, and 2 µL of Enzyme Y3, which was added to each tube of dry beads and mixed by pipetting. To each sample, 5 µL of Y2 index was added to uniquely barcoded libraries. Reactions were mixed by pipetting 60× and incubated in a thermocycler for 15 min at 25°C with the lid open. To clean up, 39 µL (1.3×) of PEG NaCl was added to each tube and mixed by pipetting 15×. Samples were incubated for 10 min at room temperature and 5 min on a magnetic rack, before beads were washed and air-dried as before. Ligation II mix was prepared by combining 30 µL of low EDTA TE, 5 µL of Buffer B1 2 µL of Reagent B2-MID, 9 µL of Reagent B3, 1 µL of Enzyme B4, 2 µL of Enzyme B5, and 1 µL of Enzyme B6. This was added to the dried beads and mixed by pipetting 30×, before incubating in a thermocycler for 10 min at 40°C with the lid open. Then, 55 µL (1.1×) of PEG NaCl was added and mixed by pipetting 30×. Samples were incubated for 10 min at room temperature and 5 min on a magnetic rack before beads were washed and air-dried as before. DNA was eluted by adding 50 µL of elution buffer and pipetting 15×. Tubes were incubated tubes at room temperature for 5 min, and then on a magnetic rack for 5 min. The eluate was transferred to a clean tube and the beads washed with a further 10 µL of elution buffer. This was collected and combined with the first 50 µL of eluate.

DNAs were size-selected using 1.2×/0.4× ratios of ProNex beads to purify ≈200–500 bp fragments. First, 72 µL of ProNex beads were added to the eluate and mixed by pipetting 30×. Samples were incubated for 10 min at room temperature and for 5 min on a magnetic rack. The supernatant (containing fragments <500 bp) was transferred to a clean tube. Then, a further 24 µL of ProNex beads were added to this supernatant and mixed 30× by pipetting. This was incubated for 10 min at room temperature followed by 5 min on a magnetic rack. The supernatant was removed and beads (with DNA fragments >200 bp bound) were washed three times with 100 µL of wash buffer, each time soaking beads for 30–60 s. The final wash buffer was removed and beads air-dried for 5 min. DNA was eluted by adding 10 µL of low EDTA TE, pipetting 10× to mix, and incubating for 10 min at room temperature. After 5 min on a magnetic rack, the eluate was transferred to a fresh tube and the beads incubated for a further 10 min with 10 µL of TE. This was collected and combined with the first collection to give a total of 20 µL.

Libraries were amplified by adding amplification mix containing 10 µL of low EDTA TE, 5 µL of Reagent 1, 4 µL of Reagent 2, 10 µL of Buffer R3, and 1 µL of Enzyme R4 to each tube. This was mixed by pipetting and incubated on a thermocycler using the following conditions: Denature1 – 98°C for 30 s, Denature2 – 98°C for 10 s, Anneal – 60°C for 30 s, Extension – 68°C for 60 s, Repeat to Denature2 for total 13 cycles, Final extension – 68°C for 60 s, Soak – 4°C forever.

Libraries were purified by adding 75 µL (1.5× volume) of ProNex beads and pipetting 60×. Samples were incubated at room temperature for 10 min, and 5 min on a magnetic rack. The supernatant was discarded and beads washed three times with 100 µL of wash buffer, before being air-dried for 5 min. Libraries were finally eluted in 20 µL of TE and their quality checked by Tapestation. They were sequenced by 75 bp single end reads to a depth of ≈20 million reads per library using the NextSeq 500 system.

OMNI-ATACseq was conducted according to Corces et al., 2017, with some adjustments, which was in turn based on Buenrostro et al., 2013; Buenrostro et al., 2015. In brief, 0, 5, 9, and 12 hr node induced and contralateral uninduced tissues were dissected in ice-cold 1× PBS and transferred to low binding PCR tubes containing 1× Protease Inhibitor Cocktail in 1× PBS (cOmplete mini EDTA-free; Roche) on ice. Where necessary, 0.12% trypsin was used to aid dissection, and treated tissues rinsed with 1× PBS to remove residual trypsin activity. Tissues were collected in batches every 30 min before tubes were spun at 100 × g/3 min/4°C to gently pellet tissues. Tubes were snap-frozen in liquid nitrogen before storage at –80°C for up to 4 months. At each time point, 14–22 pieces of tissue were collected in duplicate.

Once sufficient tissues were collected, they were thawed on ice and the nuclei isolated in 2 mL of 1× Homogenization Buffer Unstable (HBU) with a Dounce homogenizer on ice. HBU 1× contains 5 mM CaCl₂, 3 mM Mg(Ac)₂, 10 mM Tris pH 7.8, 320 mM sucrose, 0.1 mM EDTA, 0.1% NP-40, 0.02 mM PMSF, and 0.3 mM β-mercaptoethanol, dissolved in H₂O. Tissues were ground using 15 strokes with the loose pestle and 20 strokes with the tight pestle. Homogenate was then passed through a 30 µm nylon mesh filter into a fresh 2 mL low binding microfuge tube. Next, nuclei were gently pelleted by centrifugation for 10 min/500 × g/4°C and the supernatant removed. The pellet was resuspended in 200 µL of 1× HBU by gentle pipetting. Nuclei were counted by taking 6×10 µL homogenate samples, staining them with DAPI and viewing on a haemocytometer. After counting, nuclei were washed with 1 mL of ATAC-RSB + 0.1% Tween-20 and spun for 10 min/500 × g/4°C to pellet. The supernatant was carefully removed and the nuclei pellet was resuspended in 50 µL of transposition mix by pipetting. The volume of Tn5 transposase (Illumina, 20034197) added to each reaction was scaled to 0.1 µL per 1000 nuclei, so that ATAC samples with different numbers of nuclei were treated equivalently. Therefore, transposition mix contained: 25 µL of 2× TD Buffer, 0.5 µL of 1% Digitonin, 0.5 µL of 10% Tween-20, 16.5 µL of 1× PBS, X µL of Tn5 transposase (0.1 µL of Tn5 per 1000 nuclei), topped up with sterile H₂O up to 50 µL. Transposition reactions were incubated for 30 min at 37°C in a thermomixer with 1000 rpm shaking. Immediately after samples were purified using the Zymo DNA Clean and Concentrator-5 Kit, eluted using 21 µL of elution buffer and stored at –20°C until library preparation.

Libraries were prepared exactly as described (Corces et al., 2017), using the adapter sequences listed in Figure 2—source data 2 (Buenrostro et al., 2013). After qPCR amplification, qPCR profiles were manually assessed to determine the additional number of cycles to amplify according to Buenrostro et al., 2015. Libraries were purified using the Zymo DNA Clean and Concentrator-5 Kit and eluted in 11 µL of H₂O. Concentrations and banding profiles were checked by Qubit and Tapestation. Sequencing (50 bp paired-end) was conducted by Oxford Genomics to ≈15–20 million reads –sufficient for open chromatin profiling of the chick genome.

ATACseq was also performed on pre-primitive-streak stage epiblast dissected from EGKXII-XIII embryos and on mature neural plate dissected from HH6-7 embryos (15–20 pieces of tissue per sample). Samples were collected in 1× Protease Inhibitor Cocktail in 1× PBS (cOmplete mini EDTA-free; Roche) and dissociated using a Dounce homogenizer. ATACseq was performed as described (Buenrostro et al., 2013; Buenrostro et al., 2015).

Fertilized chicken eggs (Henry Stewart & Co. Ltd, Norfolk UK) were incubated at 38°C. The embryos were removed from the egg, staged, and pinned ventral side up on a resin plate in Tyrode’s saline. Embryos were dissected with a G-31 syringe needle. The endoderm and mesoderm were removed with a small volume of dispase (10 mg/mL) before dissecting an ectodermal region spanning the neural plate and neural plate border (Figure 8A). Tissue samples were pooled from multiple embryos of the same stage prior to dissociation (HH4=~55 embryos; HH6=~45 embryos; HH8=10 embryos; HH9+=7 embryos). For cell dissociation, samples were incubated in FACSmax cell dissociation solution (Amsbio, T200100) with 30 U/mL Papain (Sigma, P3125) at 37°C for 20 min. The tissue was pipetted 10 times every 5 min in order to facilitate mechanical dissociation. To stop enzymatic dissociation, an equal volume of resuspension solution (HBSS with 0.1% non-acetylated BSA (Invitrogen, 10743447)), 1% HEPES, 1× non-essential amino acids (Thermo Fisher Scientific, 11140050) and 10 μM rock inhibitor (StemCell Technologies, Y-27632) was added before passing cells through a 20 μm filter (Miltenyi Biotech, 130-101-812). Samples were then centrifuged at 100 rcf at 4°C for 5 min and resuspended in resuspension solution twice. Cells were then FAC sorted with 1 µL 0.1 mg/mL DAPI to remove dead cells and any remaining cell doublets. Immediately after FACS, samples were centrifuged and placed in 90% MeOH with 10% resuspension solution for storage.

Given the limited number of cells collected from individual dissection rounds and to obtain the required cell concentration for 10×, multiple rounds of dissection were carried out over successive days. Immediately prior to loading on the 10× Chromium controller, samples were pooled by stage and underwent three rounds of centrifugation at 1500 rcf at 4°C for 10 min followed by rehydration with DPBS with 0.5% non-acetylated BSA and 0.5 U/μL RNAse inhibitor (Roche 3335399001).

Cells were loaded on the 10× Chromium controller with the aim of capturing 3000–5000 cells per stage. cDNA synthesis and library preparation were carried out using the Chromium Single-Cell 3′ Reagent v3 Kit (10× Genomics, Pleasanton, CA, USA) according to the manufacturer’s protocol. This was carried out by the advanced sequencing facility at the Francis Crick Institute, London. Barcoded libraries were multiplexed and sequenced on an Illumina HiSeq 4000.

The initial RNAseq survey was conducted when galGal4 was the latest available chicken genome assembly. galGal4 and galGal3 compensated each other with the gene annotations (some genes were only annotated in one of the assemblies), therefore the RNAseq data was analysed with galGal3 and galGal4 as described below. Differentially expressed genes were then selected to represent on the NanoString probe codeset for further gene expression screening in triplicate with fine time course.

Raw sequencing data in FASTQ format underwent quality control analysis according to the pipeline published by Blankenberg et al., 2010. Files were converted to Sanger format using FASTQ Groomer. The quality scores of reads were calculated using FASTQ Summary Statistics. Quality scores were measured as phred = −10 log₁₀(p), where ‘p’ is the estimated probability of a base being called incorrectly. Reads with a phred score of equal to or greater than 20 (i.e. 99% probability of a base being correctly identified) were kept, while poor quality reads with average phred scores of less than 20 were filtered out using FASTQ Quality Filter. The remaining paired-end reads were trimmed at the 3’ ends where the phred score drops below 20. All reads that passed quality control analysis were initially aligned to the chicken genome (assembly version Gallus_gallus-2.1, GenBank Assembly ID GCA_000002315.1) using TopHat (Trapnell et al., 2009). The raw data was re-analysed when the chicken genome was updated to assembly version gallus_Gallus-4.0, GenBank Assembly ID GCA_000002315.2 in 2013.

To build the GRN, RNAseq data was re-analysed using galGal5 to allow incorporation of the ChIPseq, ATACseq, and scRNAseq data on the same genome assembly. The node graft time course data together with previously published embryonic neural tissue RNAseq data (Trevers et al., 2018) were re-analysed using the galGal5 genome assembly. Sequencing adapters and poor quality or N bases at 5’ and 3’ ends were trimmed from each sequence read using Trimmomatic-0.36 (Bolger et al., 2014). Unpaired reads and those that were less than 36 bases long after trimming were also removed. The remaining reads were aligned to the galGal5 genome using TopHat2 (Kim et al., 2013), alignment rates were 67.4% ± 12.9%. A custom GTF file was generated by adding additional gene ERNI (NM_001080874) annotation, taken from RefSeq, to Gallus_gallus-5.0 Ensembl release 94 GTF file. Transcripts were counted and normalized based on the custom-edited GTF file using Cufflinks (Trapnell et al., 2012) program cuffquant and cuffnorm, respectively. The output table contains transcript FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values across 11 samples (node graft four time points and embryonic neural tissues with their controls). Sequence read tracks were computed using deepTools (Ramírez et al., 2016) bamCoverage (parameters -bs 1 --scaleFactor 10⁶/Library size --extendReads --samFlagInclude 66 –ignoreDuplicates --effectiveGenomeSize 1230258557).

The initial RNAseq analysis of pooled samples was mainly used for selection of differentially expressed genes to design a probeset for NanoString expression analysis, for more detailed study of the timing of expression. Differentially expressed transcripts were identified by comparing uninduced and induced tissues at each time point (5, 9, and 12 hr). The ‘easyRNASeq’ (Delhomme et al., 2012) together with the Ensembl galGal3 Gene Transfer File (GTF) were used for transcript counts from aligned reads. Differential expression analysis was then performed using two different methods: Cufflinks cuffdiff (Trapnell et al., 2012) and DESeq (Anders and Huber, 2010). Genes that are upregulated by a log₂(FC of ≥1.2 or downregulated by ≤–1.2 in either Cufflinks or DESeq analyses, and were statistically significant with a p-value of ≤0.05 (i.e. coloured red or orange in Figure 1—source data 1) were selected). This identified 7745 transcripts that were differentially expressed across three time points (Figure 1—source data 1). Gene annotations were assigned to 4508 of these, corresponding to 2333 unique genes. Due to the incomplete nature of the galGal3 chicken genome, 3237 transcripts were left unannotated.

This process was repeated using the Ensembl and UCSC galGal4 GTF when they were released. For Ensembl transcripts, gene annotations and coordinates were obtained from Ensembl Biomart. For UCSC transcripts, gene annotations and coordinates were obtained from UCSC galGal4 GTF file and annotation data from the UCSC table browser. Fully annotated transcripts from either Ensembl or UCSC lists were combined to provide the most comprehensive set of annotations. Genes that were upregulated by a log₂(FC) of ≥1.2 or downregulated by ≤–1.2 in DESeq and were statistically significant using a p-value of ≤0.05 (i.e. coloured pale blue in Figure 1—source data 2) were selected. In situ hybridization was performed for transcriptional regulators that satisfied these criteria (Figure 1—source data 2). This reanalysis identified 8673 differentially expressed transcripts across the three time points. Gene annotations were assigned to 7184 of these, corresponding to 4145 unique genes, but 989 transcripts still remained unannotated. Volcano plots in Figure 1R–T were generated using the data generated from galGal4 assembly.

More stringent criteria were applied to select transcriptional regulators to represent on the NanoString for a more detailed study of the timing of expression. Transcriptional regulators that are upregulated at 5, 9, or 12 hr with log₂(FC)≥1.2, p<0.05, and an induced base mean ≥45 were selected from DESeq analysis were selected. Transcriptional regulators that are downregulated at 5, 9, or 12 hr with log₂(FC) ≤–1.2, p<0.05 and uninduced base mean ≥200 were also selected. A total of 213 transcriptional regulators were selected to be included in the NanoString probeset.

For the GRN pipeline, the expression data (Figure 3—source data 1) generated based on galGal5 genome assembly was analysed in R-3.5.1 environment. FC was calculated from the levels of expression in the induced relative to the corresponding uninduced tissue, or neural relative to non-neural tissue samples. Genes with FPKM >10 in an induced or neural tissue and FC >1.5 are defined as upregulated, whereas genes with FPKM >10 in an uninduced/non-neural tissue and FC <0.5 are defined as downregulated. For the time point 0, a gene is defined as upregulated when 0 hr uninduced FPKM >10, FC >1.5 against 5 hr induced tissue, and 5 hr uninduced FPKM value is over 10 and larger than the value in 5 hr induced tissue. The FPKM expression levels of the candidate GRN members are provided in Figure 3—source data 1. This gene expression information was then incorporated with the results from the NanoString analysis, which were performed in triplicate. For details of how the RNAseq data were used to build a GRN, please see the section ‘Pipeline for GRN’.

Raw NanoString data were analysed in Microsoft Excel according to NanoString guidelines. Data were checked to ensure 600 fields of view were counted and that binding density values fell within the 0.05–2.25 range. Next, the geometric mean (geomean) of six positive control probes was calculated for each assay and averaged across the entire dataset. From these, a positive lane normalization factor (PLNF) was calculated by dividing the average geomean by the geomean for each assay. Then, negative and endogenous probe counts were multiplied by their respective PLNF to normalize for differences in numbers of reporter and capture probes between assays. Next, the geomean of housekeeping (HK) genes ACTB and GAPDH was calculated for each assay and averaged across all assays. The average HK geomean was then divided by the HK geomean for each lane, to generate a lane-specific normalization factor (LSNF). Then, negative and endogenous probe counts were multiplied by the LSNF to normalize for cell number differences between assays. For each assay, the mean plus two times standard deviation was calculated from the eight normalized negative control probes. This value was subtracted from all normalized endogenous probe counts to remove background noise levels. Any transcript counts that became negative as a result were reset to zero. Next, the normalized expression counts were set to 1 when the counts are equal to zero. Average counts for each probe and standard deviation were calculated across triplicate assays. Finally, FC for each probe was calculated by dividing the induced average by uninduced average. Differential expression was defined as FC of ≥1.2 (upregulation) or ≤0.8 (downregulation). The statistical significance was calculated using a two-tailed Type 2 t-test with a p-value of 0.05. The result is provided in Figure 3—source data 1.

Single-end sequence reads were quality checked and trimmed using Trimmomatic-0.36 with the same criteria as RNAseq analysis and aligned to the galGal5 genome using bowtie2 (Langmead and Salzberg, 2012) with default parameters. The alignment rates were 92.3% ± 6.5%. PCR duplicates were removed using Piccard MarkDuplicates (parameter REMOVE_DUPLICATES = TRUE). Peak calling was computed using MACS2 (Zhang et al., 2008) callpeak (parameters -f SAM -B --nomodel --broad --SPMR -g 1230258557 -q 0.01) with genomic input as control. The MACS2 outputs were quantified and filtered for the peaks enriched with either H3K27ac or H3K27me3 by comparing to the corresponding inputs with cut-off p-value <1E-5, FC >1.2, and q-value >3. Signal tracks were computed using deepTools bamCompare (parameters --scaleFactorsMethod SES --operation log2 -bs 1 --ignoreDuplicates) and displayed using R package trackViewer (Ou and Zhu, 2019).

CTCF ChIPseq data (Kadota et al., 2017) were processed with the same pipeline. The resulting file (bed format) from peak calling was used as an input file to the pipeline for constructing our GRN.

Paired-end sequence reads were quality checked and trimmed using Trimmomatic-0.36 with the same criteria as RNAseq and aligned to the galGal5 genome using bowtie2 (parameters -X 2000 --sensitive-local). The alignment rates were 92.36% ± 5.33%. PCR duplicates were removed using Piccard MarkDuplicates (parameter REMOVE_DUPLICATES = TRUE). Coverage tracks were generated using deepTools bamCoverage (parameters -bs 1 --scaleFactor 10⁶/library size --extendReads --samFlagInclude 66 --ignoreDuplicates --effectiveGenomeSize 1230258557).

scRNAseq alignment and downstream analysis was run using Nextflow (version 20.07.1) (Di Tommaso et al., 2017) and Docker for reproducibility. All of the required packages and respective versions are found within the Docker container alexthiery/10×_neural_tube:v1.1. The full analysis pipeline, including documentation, can be found at https://github.com/alexthiery/10x_neural_tube, (Trevers et al., 2023a copy archived at swh:1:rev:30b789553714ccc4ed436f37fdf0efb4e332c888).

In order to obtain the recommended 50k reads/cell, libraries were sequenced twice, with reads from both flow cells pooled during sequence alignment. scRNAseq reads were demultiplexed, aligned, and filtered using Cell Ranger (version 3.0.2, 10× Genomics). Reads were aligned to galGal5, using a custom-edited GTF annotation file as described in ‘RNAseq analysis’ section. Prior to alignment, MT, W, and X chromosome genes in the galGal5 GTF file were prefixed accordingly for simple identification downstream. Coordinates for W chromosome genes in galGal6 were also transferred to galGal5 and prefixed in the GTF.

Quality control, filtering, dimensionality reduction, and clustering were carried out in Seurat (version 3.1.5) (Stuart et al., 2019). We first applied a modest filtering threshold, removing cells with greater than 15% UMI counts from mitochondrial genes and cells with fewer than 1k or greater than 6k unique genes. After filtering 8652 cells remain in the dataset (HH4=2745; HH6=1823; HH8=1750; HH9+=2334). Throughout subsequent clustering steps, we measure total UMI counts, total gene counts, and percentage mitochondrial content for each cluster and remove any clear outlier clusters. Following initial clustering, cell clustering was dominated by the expression of a few W chromosome genes and we therefore removed W chromosome genes. We then identified and removed contaminating cell clusters, including putative mesoderm and primordial germ cells. Percentage mitochondrial content, cell cycle, and sex were all regressed out during scaling of the final dataset.

PCA was used as an initial dimensionality reduction step, with the top 15 PCs used for clustering. We then calculated k-nearest neighbours in order to embed a shared nearest neighbour (SNN) graph. The Louvain algorithm was used to partition the SNN graph according to the Seurat default parameters. At each clustering step, we visualized multiple resolutions using the R CRAN clustree package (Zappia and Oshlack, 2018) and determined the optimal resolution manually based on cluster stability.

An unbiased approach was used to identify modules of co-correlated genes using the Antler gene modules algorithm (https://github.com/juliendelile/Antler) (Delile et al., 2019). For this we kept genes which have a Spearman correlation greater than 0.3 with at least three other genes. Ward’s hierarchical clustering is then carried out on the Spearman gene-gene dissimilarity matrix of the remaining genes. Iterative clustering by the Antler algorithm selects modules based on their consistency. In order to identify gene modules that explain the greatest amount of variation between our cell clusters, we filtered modules for which fewer than 50% of genes are differentially expressed (logFC >0.25, adjusted p-value <0.001) between at least one cell cluster and to the remaining dataset. Differential expression tests were carried out using Seurat FindAllMarkers function (Wilcoxon test).

Neural cell clusters were subset from the parent dataset using a candidate marker approach. PC1 inversely correlates with developmental stage, therefore cells were subsequently ranked according to their position along PC1. GRN components were grouped based on the timing of their initial expression following neural induction. Housekeeping genes (GAPDH; SDHA; HPRT1; HBS1L; TUBB1; and YWHAZ) were used as a control. A general additive model was used to estimate the changes in scaled expression across the ranking of PC1, with a smoothing spline fitted by REML.

To identify the core transcription regulators in the GRN, the betweenness centrality ($C\text{B}$) of each gene is computed as follows White and Borgatti, 1994:

where ∂ is the total number of directed shortest paths from one component to another component in the network, and $\partial \left(i\right)$ is the number of these directed shortest paths that pass through $i$ for each time point. Cytoscape (Shannon et al., 2003) was then used for the network representation.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

This study was funded by grants from NIH (R01 MH 60156), MRC (G0400559), Wellcome Trust (063988) and BBSRC (BB/R003432/1 and BB/K007742/1) to CDS and BBSRC (BB/K006207/1) to AS. The work of NML is supported by the Francis Crick Institute which receives its core funding from Cancer Research UK (FC010110), the UK Medical Research Council (FC010110), and the Wellcome Trust (FC010110). The scRNAseq data analyses were performed using infrastructure from the Crick Scientific Computing science technology platform. NML is a Winton Group Leader in recognition of the Winton Charitable Foundation’s support towards the establishment of the Francis Crick Institute.

© 2023, Trevers, Lu et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.