Self-organizing in vitro mouse neural tube organoids mimic embryonic development

Given the developmental importance of the neural tube and the biological richness of its morphogenesis and patterning, it is somewhat surprising that its in-depth, real-time study is hampered by the availability of only relatively few models. Indeed, studies of neural tube development to date have been possible mainly through the use of animal models, which require labor-intensive, technically challenging and low-throughput interventions. Additionally, the tracking of key cellular events in real time, which is important for understanding various cellular processes, is complicated in animals. However, 3D organoids have emerged as faithful in vitro models that recapitulate key developmental processes to better study the emergence and function of organ development with greater experimental accessibility and reproducibility in real time (Huch et al., 2017; Rossi et al., 2018; Takebe and Wells, 2019). The development of neural tube organoids and their accessibility should help address the cellular organization and connectivity of the neural network, which previously seemed impossible.

There are many specific and tightly controlled events that must take place in a model to replicate in vivo neural tube development. In the developing mammalian embryo, the neural tube emerges by the formation of a tube through an invagination of the neuroectoderm, which is derived from the embryonic epiblast after gastrulation (Sadler, 2005; Greene and Copp, 2009). The structure contains a single continuous lumen extending the entire length of the embryo, called the neural canal (Thouvenin et al., 2020), that forms the early ventricular system that is crucial for maintaining homeostasis and protecting central nervous system tissues (Fame and Lehtinen, 2020; Fame et al., 2020). The neural tube develops to specify mature, region-specific cell types along its anterior-posterior (AP) and dorsal-ventral (DV) axis in a remarkably sophisticated spatio-temporal patterning process controlled by morphogens and interactions with surrounding tissue types (Yamaguchi, 2001; Deschamps and Van Nes, 2005; Schilling, 2008; Briscoe and Small, 2015; Sagner and Briscoe, 2017, 2019). In particular, multipotent neural crest cells emerge, delaminate and migrate from the dorsal side of the neural tube and differentiate into a wide variety of cell types that contribute to the peripheral nervous system, the endocrine system and the facial skeleton, among many others (Le Douarin et al., 2004; Dupin, 2006). The abnormalities in their development can cause congenital neurocristopathies such as Treacher Collins (Vega-Lopez et al., 2018; Sato et al., 2019; Pilon, 2021). Even though several in vitro neural tube models have been reported recently (Meinhardt et al., 2014; Ranga et al., 2016; Zheng et al., 2019; Ogura et al., 2018; Duval et al., 2019; Veenvliet et al., 2020; Libby et al., 2021; Demers et al., 2016; Rifes et al., 2020), to date, there does not appear to be a suitable model that recapitulates the in vivo morphology, coordinated AP patterning and emergence of the neural crest.

We present here a self-organizing neural tube organoid that mimics the mouse embryonic neural tube in terms of its morphology, cell type composition and aspects of patterning. When exposed to a sequence of epiblast culture conditions and neural differentiation conditions in 3D Matrigel, single mouse embryonic stem cells (mESCs) develop into spontaneously elongating neuroepithelial tissues organized with key hallmarks of the neural tube. Single cell transcriptomics revealed the presence of regionalized cell types spanning the region from the midbrain to the spinal cord along the AP axis of the embryo. In long-term culture experiments for organoid maturation, we observed the emergence of multipotent neural crest cells and mature neurons. Tracking this development by single-cell transcriptomic analysis at serial time points revealed the transcriptomic details of neural crest cell formation and neural differentiation. Collectively, the accessibility and scalability of these organoids make them a potent model of the developing neural tube, allowing the study of key mechanisms involved in self-organization, cell fate specification and neural crest cell development in vivo.

To recapitulate neural tube development in vivo, we postulated the generation of a self-organizing neural tube organoid derived from the epiblast. To generate such structures, we first embedded a suspension of single mESCs in Matrigel and promoted their differentiation towards epiblast fates and neural lineages, sequentially. Specifically, to physiologically mimic the neuroepithelium development in vivo, the cells were cultured for the first 3 days in serum-free N2B27 medium containing bFGF and activin A (Brons et al., 2007; Hayashi et al., 2011) to guide them toward an epiblast state. The cells were then cultured in N2B27 to induce neural differentiation (Fig. 1A). We observed that mESCs form cysts and spontaneously elongate into structures harboring a continuous lumen and displaying characteristic features of the embryonic neural tube (Fig. 1B, Fig. S1, Movie 1). Time-lapse imaging between day (D) 4 and D6 using a live cell actin probe to capture the lumen of apicobasally polarized organoids demonstrated that the neuroepithelial cysts indeed elongate and such elongation was concomitant with that of their lumen (Fig. S1B, Movie 1). Time-course images from D0 to D6 confirmed that self-elongating organoids emerged from a single cell (Movie 1). Elongation efficiency was 52.0±4.5% (n=9 Matrigel drops, N=3 replicates) (Fig. S1C) and the length and elongation index increased about threefold between D4 and D6 (Fig. S1D). Neural tube organoid formation was reproducible in three different mESC lines that were tested, each line requiring an optimal cell seeding density (Fig. S2, Table S1).

To assess the similarity of the organoids to the embryonic neural tube, we examined the expression of key markers of neural tube development in vivo. This revealed progressive differentiation from a pluripotent state toward neuroprogenitor identities. Specifically, quantitative PCR analysis of organoids between D3 and D6 shows a gradual transition from epiblast to neuroprogenitor fate, with epiblast markers Pou5f1 (also known as Oct4), Otx2 and Fgf5 gradually decreasing over time, and neuroprogenitor markers Sox1, Sox2 and Pax6 gradually increasing as the organoids develop (Fig. 1C).

In addition to neuroprogenitor markers, our organoids also exhibit key hallmarks of the tissue architecture of the embryonic neural tube. For example, they also possess apicobasal polarity, as demonstrated by the apical localization of the tight junction marker ZO-1 and the presence of apical domain markers prominin 1 and PKCζ (Fig. 1D). Immunofluorescence staining confirmed the presence of neuroprogenitor cells expressing the markers SOX1, SOX2, PAX6 and nestin (Fig. 1E). In addition, confocal imaging revealed the presence of primary cilia (marked by ARL13B) at the apical membrane of the organoids (Fig. 1F, Movie 2), which are essential for the sonic hedgehog (Shh) pathway and the planar cell polarity (PCP) pathway during neural tube development in vivo (Bay and Caspary, 2012). Moreover, light-sheet microscopy showed a highly organized 3D structure of the neuroepithelium bounding the elongated lumen, which confirms that this is indeed continuous throughout the organoid (Fig. 1G, Movie 3). This elongated tubular lumen resembles the neural canal of a developing neural tube, which forms the early ventricular system that is essential for tissue homeostasis of the nervous system. Furthermore, the basal localization of EdU-labeled DNA synthesis phase (S phase) cells and the apical localization of pH3⁺ mitotic cells (M phase) are indicative of characteristic division patterns of interkinetic nuclear migration, a defining feature of cell behavior in the embryonic pseudostratified neuroepithelium (Baye and Link, 2007; Miyata et al., 2015) (Fig. 1G, Fig. S3, Movie 4). Quantification of the distance of M-phase cells (pH3⁺) and S-phase cells (EdU⁺) from the apical domain (prominin 1⁺) showed the apical localization of pH3⁺ mitotic cells (Fig. S3A). Time-lapse imaging captured and confirmed the interkinetic nuclear migration movement of the cells in the neural tube organoids (Fig. S3B, Movie 4). Altogether, these results suggest that the neural tube organoids both physiologically and morphologically recapitulate key features of developing neural tube in vivo.

To better characterize the cellular diversity of organoids, we performed single cell RNA sequencing (scRNAseq) of cells isolated from D6 organoids (n=4478 cells, N=2 replicates). A clustering analysis using a shared nearest neighbor modularity algorithm followed by an analysis of the markers characterizing each cluster identified a total of six cell types (Fig. 2A,B). The organoid is predominantly composed of neuroprogenitor cells from different AP axis regions of the neural tube, which include the midbrain (En1⁺ and En2⁺), hindbrain (Crabp1⁺ and Egr2⁺) and spinal cord (Pax6⁺ and Hoxb9⁺). The remainder of the cells consists of the caudal epiblast (Pou5f1⁺, Cdx2⁺ and Wnt8a⁺), surface ectoderm (Krt8⁺ and Krt18⁺) and mesenchyme (Twist1⁺ and Snai1⁺) (Table S7).

To verify the spatial organization of AP axial identity cell types in a single organoid, we performed single-molecule fluorescence in situ hybridization (smFISH). We examined two pairs of genes: Otx2 (midbrain) and Gbx2 (hindbrain); and Hoxb1 (hindbrain) and Hoxb9 (spinal cord) (Fig. 2C). Co-expression levels of Otx2 and Gbx2 projected on UMAP showed that these two genes are expressed in a mutual exclusive manner (Fig. 2D). smFISH revealed their spatial organization along the AP axis of the organoid (Fig. 2E; Otx2 in green and Gbx2 in magenta). On the other hand, a higher level of Hoxb1 and Hoxb9 co-expression was observed, which was confirmed by their spatial arrangement visualized using smFISH (Fig. 2G). Hoxb1 signal spots were detected throughout the organoid, whereas Hoxb9 signal spots were localized at one end of the organoid, indicating the posterior part of the organoid. The expression of a broad spectrum of Hox genes along the AP axis in the organoid analyzed by scRNAseq and qPCR supports these results (Fig. 2F, Fig. S4). Indeed, as shown in Fig. S4A, anterior Hox genes were expressed in the majority of cells that made up the organoids, whereas posterior Hox genes were expressed in more limited cell populations. Crucially, the range of Hox gene expression in the organoids extended from Hox1 to Hox9, corresponding to regions of the embryonic neural tube ranging from the cervical to the beginning of lumbar levels. The same range of Hox gene expression was confirmed by qPCR analysis (Fig. S4B). From these results, we conclude that the genes for AP axis identity are spatially organized in a single organoid. Importantly, we were able to obtain such patterning without the need for ‘artificial’ posterization, e.g. by using CHIR99021 or retinoic acid (RA), as in previous studies (Meinhardt et al., 2014; Ranga et al., 2016; Turner et al., 2017; Beccari et al., 2018; Veenvliet et al., 2020). Immunofluorescence imaging over time against the neuromesodermal progenitor (NMP) markers SOX2 and T/brachyury from D3 to D6, and the results of scRNAseq analysis also confirmed that there was no NMP in the organoid (Fig. S5). These results indicate that the axial elongation of the neural tube organoid is not driven by NMPs, as is the case with other axially elongated in vitro tissue models, such as gastruloids (Turner et al., 2017; Beccari et al., 2018; Veenvliet et al., 2020).

To understand the difference in cell type composition between the organoids that are elongated (D6 organoids) and not elongated (D6 NE organoids), we performed scRNAseq of these two organoid types (cultured under the same conditions) at day 6 (Fig. S6). One cell type that was preferentially present in D6 organoids compared with the D6 NE organoids was the spinal cord (24.4% of the D6 and 8% of the D6 NE organoids, respectively) (Fig. S6A). Visualization of Hox gene expression in these two groups reveals that the D6 organoids have relatively stronger expression of posterior Hox genes, whereas the D6 NE organoids express anterior genes more strongly (Fig. S6B). These data suggest that there may be a relationship between axial elongation and Hox gene expression that could explain the elongation of D6 organoids rather than D6 NE organoids. However, further studies on the mechanisms involved are needed to clarify the relationship between elongation and the development of spinal cord cell identities.

As the neural tube develops and expands, there are two major developmental events. First, neural crest cells emerge from the dorsal side of the neural tube, which then migrate out and contribute to the development of various tissues from the peripheral nervous system to the facial skeleton. Second, neuroprogenitors of the neural tube differentiate into more mature neurons, such as motor neurons and interneurons. As our neural tube organoids exhibit some of the hallmarks of the native neural tube, we hypothesized that they would also retain the capacity to generate neural crest cells and mature neurons in vitro. To test this, we harvested and re-embedded in a Matrigel drop a single organoid and tested whether it could further mature from D6 onwards (Fig. 3A). Time-lapse imaging of re-embedded organoids from D6 to D8 revealed two main types of maturation events: neural crest cells and neurons (Fig. 3B, Movie 5). By quantifying the proportion of emergence of these two cell types from D7 to D10, we observed that neural crest cell emergence happened earlier than neurite outgrowth (Fig. 3C). From D9 onward, the majority of organoids had both neural crest cells and neurons. These two modes of organoid maturation were further characterized.

First, the process of the outward emergence and migration of neural crest cells in 3D from the organoids was captured by the time-lapse imaging (Fig. 3B, left box, Movie 5). Its adequate accessibility allowed us to observe the delamination, migration, cell division and phagocytosis of neural crest cells (Fig. S7, Movie 6). Especially, it was intriguing to observe the phagocytosis event of neural crest cells, which has not been reported in mammals, only in zebrafish, due to its limited accessibility (Zhu et al., 2019) (Fig. S7C, Movie 6). This clearly suggests that our model could provide a new approach for answering questions in the field. Those migratory cells were later confirmed to be neural crest cells expressing neural crest cell markers such as AP2α, nestin and SOX10 by immunofluorescence on D8 (Fig. 3D, Movie 7). Furthermore, culture of the neural crest cells for 1 month revealed the multipotency of the neural crest cells (Fig. S8). The derivation of mature neurons (βIII tubulin⁺), peripheral nervous system neurons (peripherin⁺), glia (GFAP⁺ and S100β⁺), adipocytes (bodipy⁺ indicating lipid droplet) and smooth muscle cells (SM22α⁺) from neural crest cells suggests that the neural crest cells derived from organoids are multipotent, with the capacity to differentiate into not only ectodermal but also mesenchymal cell types.

The other type of maturation involved neural differentiation (Fig. 3B, right box, Movie 5). Immunofluorescence on D10 showed that these neurons were TUJ1⁺ (neuron marker) and BRN3a⁺ (dorsal neuron marker), indicating that the neural progenitors of the neural tube were differentiated into dorsal neurons. As shown in previous studies, neural tissue differentiated from pluripotent stem cells without ventral morphogens has a default dorsal fate (Meinhardt et al., 2014; Ranga et al., 2016), and our organoids cultured without ventral morphogens lack the ventral cell types of the neural tube (Fig. S9), supporting BRN3a⁺ neurons. Of note, ventralized neural tube organoids can be obtained by exposing organoids from D3 to D5 to 500 nM smoothened agonist (SAG, ventral signal) (Fig. S9A). qPCR analysis of D6 organoids for marker genes along the DV axis showed that control organoids have default dorsal identity with low expression of ventral marker genes, but ventralized organoids show upregulated ventral marker genes (Fig. S9B). Immunostaining against the ventral markers FOXA2 and OLIG2 on D6 organoids also showed that the organoids can adopt a ventral identity in response to ventral morphogens (Fig. S9C). Altogether, these data demonstrate that our organoid can be a new tool for studying the emergence, delamination and migration of neural crest cells, and the neural differentiation of the neural tube in an easily accessible and tractable fashion.

To better understand neural tube organoid maturation at the single cell level, we performed scRNAseq on D6, D7, D8 and D10 organoids, identifying the key cell types that emerge during maturation (Fig. 4A, Fig. S10, Table S8). The transition in the ratio of cell types from D6 to D10 shows a similar trend to the earlier observation that neural crest cells appear before neurons (Figs 4B, 3C). The major cell types of D6 organoid are mid/hindbrain and spinal cord cells, as shown before. Noticeably, we also identified pre-EMT (epithelial-mesenchymal transition) neural crest cell clusters (Bmp6⁺ and Msx1⁺) (Fig. 4A,C,D). Key to determine the regional AP identity of the pre-EMT neural crest is the expression of Hox genes, and both cranial and vagal/trunk neural crest clusters were identified (Fig. 4C, Fig. S11A). On D7, migratory neural crest cells (Sox10⁺ and Foxd3⁺) appear as a trajectory of pre-EMT cranial neural crest cells, and because they show low expression of Hox genes, they were identified as migratory cranial neural crest cells (Fig. 4A,C,D, Fig. S11A). Interestingly, single-cell transcriptomic analysis captured the sequence of transcriptional events during neural crest cell development. Neural crest cells delaminate by EMT between emergence and migration, and key markers of delamination (Dlx5⁺ and Pdgfra⁺) could be detected in cells in the pathway between the pre-EMT neural crest cells and the migrating neural crest cell cluster (Fig. S11B). Moreover, the cells of the pre-EMT neural crest cell cluster express roof plate markers, such as Lmx1a, Msx1, Msx2 and Wnt1 (Delile et al., 2019) (Fig. S12), suggesting that neural crest cells arise from roof plate cells as in vivo. To the best of our knowledge, this is the first time such in vivo-like transcriptomic events of neural crest development have been detected in an in vitro model.

In addition, single-cell transcriptome analysis revealed neural differentiation in the mature neural tube organoid in detail. Consistent with the standard dorsal fate of the neural tube organoid, two identified neural progenitor clusters in mature neural tube organoids were the most dorsal progenitor domains of the neural tube: dorsal progenitor 1 (dp1) and dorsal progenitor 2/3 (dp2/3) (Fig. 4A,C). All of these progenitors collectively express the neural progenitor marker Sox2 and are distinguished by their different markers: Atoh1 and Olig3 for dp1, and Pax7 and Ascl1 for dp2/3 (Fig. 4D, Fig. S13). These domains differentiate into dorsal interneurons: dp1 into dorsal interneuron 1 (dI1), and dp2/3 into dorsal interneuron 2 (dI2) and dorsal interneuron 3 (dI3). All dorsal interneuron clusters expressed both of the neuron markers Tubb3 and Elavl3; dI1 neurons could be identified by the markers Lhx2 and Lhx9, dI2 neurons by Lhx1 and Lhx5, and dI3 neurons by Tlx3 and Prrxl1 (Fig. 4C,D, Fig. S13). Another marker for dorsal neurons, Pou4f1 (also known as Brn3a), is also expressed in the dorsal interneuron clusters, and its expression was shown in the D10 organoid by immunofluorescence in Fig. 3D.

Overall, scRNAseq revealed the cell types and transcriptomic features of neural crest and neural development during maturation of the neural tube organoid. The migrating neural crest cells observed in the mature organoid were identified as cranial neural crest cells, and both cranial and vagal/trunk pre-EMT neural crest cells were detected. In addition, neurons in the mature organoid were identified as dorsal interneurons derived from dorsal neural progenitors. We believe that these results demonstrate that the neural tube organoid faithfully recapitulates the neural crest and neural development in the embryo.

Here, we present a novel organoid model that faithfully reproduces key features of the developing mouse embryonic neural tube, such as its elongated tubular architecture with a neural canal, AP patterning and the emergence of neural crest cells and neurons. Our system allows real-time observation of cell dynamics during morphogenesis and elongation, as well as patterning along the elongation axis with the emergence of region-specific cell identities.

A key parameter to promote self-elongation and -patterning in our model is the relatively low initial cell seeding density, which is ∼30% of that reported in previous neuronal cyst models (Meinhardt et al., 2014; Ranga et al., 2016). When stem cells are cultured at the cell density specified in previous neural cyst protocols, the tissue remains cystic, i.e. spontaneous expansion does not occur. The effect of cell density on organoid development should be further investigated.

Another intriguing aspect of our neural tube organoid is that AP-patterning occurs spontaneously, i.e. without ‘artificial’ posterization culture conditions induced by e.g. CHIR99021 or RA (Meinhardt et al., 2014; Ranga et al., 2016; Turner et al., 2017; Beccari et al., 2018; Veenvliet et al., 2020). For example, the generation of gastruloid or trunk-like structure (TLS) requires a CHIR99021 pulse that promotes emergence of SOX2⁺T/brachyury⁺ NMP-like cells that drive axial elongation of the system (Veenvliet et al., 2020), as in the tail bud in vivo. Reporter gene expression analysis, immunostaining and scRNAseq experiments do not indicate the presence of such NMP-like cells in our model. Considering that NMPs at the tail bud of the embryo contribute to secondary neurulation, forming the posterior spinal cord and expressing posterior Hox genes, and that our neural tube organoid lacks a NMP population and posterior Hox gene expression, we believe that elongation and AP patterning in our system occur via a different mechanisms (Copp et al., 2003; Deschamps and Duboule, 2017; Nikolopoulou et al., 2017).

Another important aspect of neural tube development is patterning along the DV axis. Our neural tube organoids have a default dorsal fate, as reported in previous studies (Meinhardt et al., 2014; Ranga et al., 2016). However, we have shown that the neural tube organoid responds to ventral morphogens and can differentiate into ventral cell types in response to the morphogens, in the same way that it does in vivo. This result highlights the potential of our organoid model for studying patterning and cell type derivation with morphogen gradients in vitro.

Finally, the accessibility of this system allows the emergence, delamination and migration of neural crest cells to be studied in real time. In particular, we could observe a phagocytosis event of migrating neural crest cells in real-time that has never been reported in a mammalian model due to its poor accessibility. As many neurocristopathies are caused by abnormalities in migration of neural crest cells, dissecting the migration of neural crest cells with this tractable system would allow us to better understand the mechanisms involved in these pathologies. We expect our organoids to open interesting avenues to the study of tissue mechanics, patterning, and neural crest and neural development in the mammalian neural tube at the cellular and molecular levels. Collectively, our system should contribute to understanding and facilitating the exploration of mammalian neural tube developmental processes and neurodevelopmental diseases in a tractable and potentially scalable manner.

Sox1^eGFP;Bra^mCherry double reporter (SBr; Deluz et al., 2016) mouse embryonic stem cells (mESCs) were used for the main experiments, except those shown in Fig. S2. Briefly, the mESC knock-in reporter line SBr was derived from CGR8 mouse ES cells with homology directed repair (HDR) and the fluorescent protein was cloned into Sox1 and Bra loci with 5′-P2A sequence. The cells were cultured in DMEM supplemented with 10% serum medium [10% ES-grade fetal bovine serum (FBS), 0.1 mM non-essential amino acids (NEAA), 1 mM sodium pyruvate, 3.44 mM GlutaMax, 0.1 mM β-mercaptoethanol (β-ME), 100 U/ml penicillin and 100 μg/ml streptomycin] with 2i (3 μM CHIR99021 and 1 μM PD025901) and 100 μg/ml LIF in a humidified incubator (5% CO₂, 37°C). The cells were split every other day, by washing with PBS (Gibco, 10010056) and dissociating with Accutase for 2 min at room temperature. The cell suspension was collected in a tube and centrifuged at 200 g for 5 min. After centrifugation, the cell pellet was resuspended in 1 ml of 10% serum medium and the cells were counted with a hemocytometer. 85,000 cells were plated in each well of a six-well plate (Corning) in 2 ml of 10% serum medium with 2i and LIF and cultured in the humidified incubator. The cells were passaged every 2 days.

mESCs were collected by washing with PBS and dissociating with Accutase for 2 min at room temperature. The cell suspension was collected in a centrifuge tube and was centrifuged at 200 g for 5 min. After centrifugation, the cell pellet was resuspended in 1 ml of 10% serum medium and the cells were counted with a hemocytometer. After cell counting, the cells were washed twice with PBS and were resuspended in N2B27 (47.4% Neurobasal medium, 47.4% DMEM/F-12 with 2.50 mM GlutaMAX, 1 mM GlutaMAX supplement, 100 U/ml penicillin, 100 μg/ml streptomycin, 0.1 mM NEAA, 1 mM sodium pyruvate, 0.1 mM β-ME, 1% B27 supplement and 0.5% N-2 supplement) and embedded in Matrigel to a final concentration of 1125 cells/10 μl drop of Matrigel for SBr. 10 μl of the resulting Matrigel-embedded cell suspension were deposited in each well of a 24-well plate (Corning). The cells in each well were cultured in 500 μl of Epi medium [N2B27 supplemented with 12 ng/ml basic fibroblast growth factor (bFGF, Thermofisher, PMG0031) and 20 ng/ml activin-A (R&D Systems, 338-AC-010)] for the first 3 days. On D3, the Epi medium of each well was substituted with 1 ml of N2B27. The medium of each well was substituted with 1 ml of N2B27 on D5 and the cells were cultured until D6. To test the reproducibility of self-elongating neural tube organoid formation, we tested three other mESC lines: (1) wild-type of the SBr cell line (Deluz et al., 2016) (SBr wild type in Fig. S2), (2) Sox1^GFP mESCs (Aubert et al., 2003; Meinhardt et al., 2014) and (3) Sox10^GFP mESCs (Kawaguchi et al., 2010) (Fig. S2). All the cell lines were maintained and cultured following the same protocol except for the initial cell embedding density in the Matrigel drop on D0 (Table S1).

To harvest the organoids, on D6, the Matrigel drops with organoids were washed once with PBS and then treated with 3 mg/ml collagenase IV and 1 mg/ml dispase in HBSS+/+ for 30 min on a shaker to dissolve the Matrigel. After Matrigel was dissolved, the organoids were collected with cut 1000 μl pipette tips coated with 10% serum and centrifuged at 150 g for 5 min. The supernatant was removed and the organoids were suspended in N2B27 medium. The organoids were collected and re-embedded in 10 μl of Matrigel (1 organoid/drop), plated on glass-bottomed 96-well plate and further cultured in 100 μl of N2B27 medium in 5% CO₂ and 37°C incubator, and the medium was changed every day.

The organoids were collected from Matrigel drop by washing with PBS and dissolving the Matrigel using 4 mg/ml collagenase-IV and 100 μg/ml dispase I in HBSS+/+ on a shaker in the incubator for 40 min. The organoids were fixed in 4% paraformaldehyde (PFA) in PBS for 30 min at room temperature. The organoids were washed three times with PBS for 20 min each, permeabilized in 0.2% Triton X-100 in PBS for 1 h and blocked in blocking solution [3% bovine serum albumin (BSA) and 0.01% Triton X-100 in PBS] for 3 h to overnight. Primary antibody (Tables S3, S5) solution was prepared in Blocking solution, with 2 μg/ml of DAPI. The organoids were stained overnight on a low-speed orbital shaker at 4°C. The organoids were thoroughly washed with PBS every 30 min for 4 h. Secondary antibody (Tables S4, S5) solution was prepared, treated and washed the same way as primary antibody solution. Here, the immunostaining was processed to avoid the green channel apart from SOX1 antibody because there are always SOX1-GFP⁺ cell populations in the organoid. Elsewhere, it was processed in the far-red channel or red channel, as there are no brachyury-mCherry⁺ cells in the organoid (Fig. S5A). Thus, the color of the immunostaining images in the figures do not necessarily match the color of the channel for the secondary antibody that was used for the staining. Confocal images of organoids were generated using LSM700 (Zeiss) on an inverted Zeiss AxioObserver Z1 (Bioimaging and Optics Core Facility, EPFL), equipped with 10×/0.30 air, 20×/0.80 air and 63×/1.40 oil objectives, 405 nm, 488 nm, 555 nm and 639 nm lasers, and controlled by ZEN 2010 imaging software (Zeiss). Images were analyzed using FIJI (NIH).

The organoids were collected, fixed and immunostained following the immunostaining method described above. CUBIC-mount [sucrose (250 g of 50%, w/v), 125 g urea (25%, w/v) and 125 g N,N,N′,N′-tetrakis (2-hydroxypropyl)ethylenediamine (25%, w/v) dissolved in 150 ml of distilled H₂O and brought up to 500 ml (Lee et al., 2016)] was diluted with water to reach a refractive index of 1.453. The organoids were embedded in 1% agarose solution in a 1.5 mm diameter capillary (701908) using a plunger and were cooled in 4°C to allow for agarose gelation. The glass capillaries with agarose-embedded organoids were dipped in CUBIC-mount overnight in a dark room at room temperature for tissue clearing. For imaging, the glass capillary with sample was mounted with sample holder and inserted into the Top-fed Zeiss Lightsheet Z1. The cell chamber was filled with CUBIC-mount solution and the organoids in the capillary were located in the field of view using the plunger. Images were acquired with a 20×/1.0 corr water W Plan Apochromat objective (for non-cleared samples) or a 20×/1.45 corr clearing EC Plan Neofluar objective (for cleared samples). DAPI, Alexa Fluor 488, Alexa Fluor 568 and Alexa Fluor 633 were excited with 405, 488, 555 and 639-nm lasers, respectively. Images were analyzed using IMARIS (Oxford Instruments) and FIJI (NIH).

To observe how the organoid and its lumen elongate between D4 and D6, the organoids of each well were treated with 1 ml of N2B27 with 500 nM Sir-actin (Table S5) for 4 h before imaging. Time-lapse images were acquired using either a Nikon Eclipse Ti inverted microscope system [equipped with a 10×/0.30 air objective, 632-nm filters, DS-Qi2 and Andor iXon Ultra DU888U (Oxford Instruments) cameras and controlled by NIS-Elements AR 5.11.02 (Nikon) software] or a Visitron CSU-W1 (equipped with a U PLAN S APO 60×/1.42 oil objective, sdc Cy5 filter, 50 μm Spinning Disk, ImagEMX2 camera and controlled by VisiView software) and a 647 nm laser for Sir-actin signal and bright-field illumination. The image was captured every hour between D4 and D6 (Nikon Eclipse Ti) or every 5 min between D5 and D6 (Visitron). Temperature was maintained at 37°C with 5% CO₂.

Time-lapse images were analyzed using FIJI. Bright-field images of organoids were segmented by thresholding and the Max Inscribed Circles function (FIJI function developed by O. Burri and R. Guiet, BIOP, EPFL) was used to fit circles in the segmented object. Axial length was determined by connecting the centers of the fitted circles. Elongation index was calculated by dividing axial length by the diameter of the largest inscribed circle determined by the Max Inscribed Circles function.

The collected organoids were washed with PBS and resuspended in lysis buffer. RNA was extracted immediately using the RNAeasy micro kit (Qiagen, 74004) according to the manufacturer instructions with on-column DNase I (Qiagen) digestion. The extracted RNA concentration was measured using NanoDrop. cDNA was synthesized using the iScript cDNA Synthesis Kit (BioRad, 170-8891) according to the manufacturer's instructions using 1 μg of RNA of each sample replicate in a BioRad thermocycler. Quantitative PCR analysis of mRNA levels for the genes on primer list was performed using Power SYBR Green PCR Master Mix (ThermoFisher Scientific, 4367659) on the Applied Biosystems Quantstudio 6 Flex Real-Time PCR system. PCR primers were designed using NCBI Primer-Blast software (Table S2). Expression values of each gene were normalized to the Gapdh expression value of each sample replicate using the ΔΔCT method. The graph was plotted using Prism 8 (GraphPad) and error bars represent the s.d. across three biological replicate samples.

The sample preparation was performed according to the instructions from 10X Genomics Single Cell Protocols (10X Genomics), with some modification. The collected organoids (100 organoids per condition) were washed twice with PBS and trypsinized with TrypLE Express (Gibco, 12605-028) for 7 min in 37°C. The organoids were dissociated into single cells by pipetting up and down with a 1000 μl pipette tip, strained using 40 μm cell strainer and collected in a 1.5 μl microcentrifuge tube. Those single cells were centrifuged at 300 g for 5 min and resuspended in 1 ml of PBS with 0.04% BSA after removing the supernatant. The cell concentration was then determined. After the same centrifugation step was repeated, the cells were resuspended in PBS with 0.04% BSA to achieve a final cell concentration of 1 million cells/ml. After the trypsinization, the cells were always placed on ice.

Single-cell transcriptomics profiling derived from the organoids from different timepoints (D6, D7, D8 and D10) was carried out using Chromium Single Cell Gene Expression (10X Genomics), according to the manufacturer's recommendations. Cellular suspension was loaded on a Chromium Controller instrument to generate Gel Beads-in-emulsion (GEMs), containing single cells. Reverse transcription and library preparation were performed using the Single Cell 3′ Reagent Kit v3 (10X Genomics), following the 10X Genomics protocol. The cells were sequenced at an average depth of 85,000 reads per cell.

Single cell RNA sequencing data analysis was performed using Seurat 4.1.0 (Stuart et al., 2019) on R, with raw feature-barcode matrices of the organoids from different timepoints (D6, D7, D8 and D10). Genes expressed in fewer than three cells and the cells with less than 200 features expressed were removed. Cells with unique feature counts between 3000 and 10,000, and less than 10% mitochondrial counts were used for analysis. Expression level was log normalized, and the top 2000 variable genes were calculated using the Variance Stabilizing Transformation (VST) algorithm. The datasets were integrated using ‘FindIntegrationAnchors’ and ‘IntegrateData’ function, and normalized by the ‘LogNormalize’ method with the dimension of 1:30 from CCA to specify the neighbor search space. After the integration, the effects of the batch ID, cell cycle, nUMI and mitochondrial counts were determined by regression analysis. The number of cells after filtering was 2076 (D6, self-elongating neural tube organoid replicate 1), 2402 (D6, self-elongating neural tube organoid replicate 2), 2225 (D6 NE, non-elongating neural tube organoid), 2565 (D7), 3337 (D8) and 3677 (D10). We scaled the data by applying a linear transformation and performed principal component analysis (PCA) on the scaled data. Clusters of cells were identified by a shared nearest neighbor (SNN) modularity optimization-based Louvain algorithm with a resolution of 1.0 (Waltman and van Eck, 2013). Clusters were visualized on a two-dimensional Uniform Manifold Approximation and Projection (UMAP) (Mcinnes et al., 2018 preprint). The cell type of each cluster was identified based on the markers of each cluster expressing more than 0.25 log₂ fold-change different over the others. To estimate RNA velocity of each cell, velocyto.py (La Manno et al., 2018) was used to generate loom files for each dataset from D6 to D10 to be used for Seurat analysis. To do so, the bam file of each dataset produced by Cell Ranger (10X Genomics) and the gtf file from 10X Genomics were used for annotation. To integrate the loom file and the Seurat object in R, metadata-like filtered cell IDs, UMAP coordinates and clusters were extracted. Once the loom files were generated, the loom files and the metadata were integrated, and RNA velocity could be estimated based on the same UMAP coordinates and clusters using scVelo (Bergen et al., 2020) on Python. For cell type identification, we have referenced mainly two single cell transcriptomic analysis reports in the developing mouse neural tube: spatial and temporal dynamics of gene expression in spinal cord (Delile et al., 2019), and neural crest development (Soldatov et al., 2019).

20 nt fluorescently labeled DNA probes for detecting Otx2, Gbx2, Hoxb1 and Hoxb9 transcripts were designed and ordered according to Stellaris (LGC Biosearch) instructions (Table S6). The probes were resuspended in 1× TE buffer [10 mM Tris-HCl and 1 mM EDTA (pH 8.0] to stock concentration 12.5 μM and aliquoted. The organoids were collected, washed once in 1× PBS and fixed in 4% PFA for 30 min. Next, whole-mount organoids were placed on SuperFrost slides (ThermoFisher Scientific) in a drop of 1× PBS. After the buffer was evaporated (∼30 min), organoids were attached to the slide for whole-mount staining. The organoids were permeabilized overnight in 70% ethanol at 4°C. An ImmEdge PAP pen (Vector laboratories) was used for creating a hydrophobic barrier around areas of interest. Next, the organoids were washed with 2× SSC twice. The tissue was incubated with smFISH hybridization buffer [2× SSC (Sigma Aldrich), 10% (w/v) dextran sulfate (Sigma Aldrich), 10% (v/v) Formamide (Ambion), 1 mg/ml E. coli tRNA (Roche) and 2 mM ribonucleoside vanadyl complexes (RVC) (Sigma Aldrich)] for 5 min without probes for buffer substitution. smFISH probes were resuspended in smFISH hybridization buffer to a final concentration of 250 nM to create the hybridization mix. Hybridization mix was incubated with the tissue for 4 h at 37°C in humidified environment. After hybridization, the tissue was washed three times for 10 min each with 20% formamide in 2× SSC at 37°C. Hoechst 33342 (ThermoFisher Scientific) staining was applied on the second wash to visualize nuclei. Next, the slides were washed twice with 2× SSC. The slides were mounted with ProLong Gold mounting media and imaged with an Upright Leica DM6 CS in wide-field imaging mode equipped with a DFC 7000 GT (CCD Grayscale) camera and a HC PL APO 40×/63× objective, and the following filters were used: DAPI (BP 350/50 excitation, BP 460/50 emission), RHOD (BP 546/10 excitation, BP 585/40 emission) and Y5 (BP 620/60 excitation, BP700/75 emission). To detect the smFISH spots in the image, an ImageJ plug-in RS-FISH (https://github.com/PreibischLab/RS-FISH) was used.

All data shown in column graphs are expressed as mean±s.d. Statistical analysis between two columns was performed using a two-tailed unpaired Student's t-test, whereas data containing more than two experimental groups were analyzed with one-way ANOVA followed by Tukey's multiple comparisons test or with two-way ANOVA followed by Šídák's multiple comparisons test. Statistical significance was calculated using the Graphpad Prism software, which was also used to generate all the graphs. *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (confidence intervals 95%; alpha level 0.05).

We thank R. Guiet and O. Burri of the BioImaging and Optics Core Facility (BIOP) for programming image-processing plug-in for elongation index quantification, and all the other staff for their help for the general support in imaging; all the staff from the Gene Expression Core Facility (GECF, EPFL) for their support on single-cell RNA sequencing (scRNAseq); S. Vianello, M. Leleu and N. Broguiere for the introduction to the scRNAseq analysis; A. Martinez Arias for help with defining cell types during scRNAseq analysis; G. Rossi for the introduction to Lightsheet fluorescence microscopy; and M. Girgin, S. Vianello, N. Brandenberg, A. Martinez Arias and E. Tanaka for the inputs on the manuscript. Some of the text and figures in this paper formed part of JiSoo Park's PhD thesis in the Doctoral Program in Biotechnology and Bioengineering (EDBB) at the Ecole Polytechnique Fédérale de Lausanne in 2021.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201052.reviewer-comments.pdf.