An in vitro model of human anterior neurulation.

(A) Micrographs of human embryo sections from the Virtual Human Embryo project showing sequential neural plate bending, neural fold elevation, and neural tube closure from Carnegie Stages (CS) 8, 10, and 11. Tissues are outlined according to epithelial thickness, denoting the approximate location of neural ectoderm (green), neural plate border (gray), and surface ectoderm (magenta). Scale bar: 300 µm.

(B) Epifluorescence images of seven hexagonally arranged neural tube organoids 1, 2, and 3 days after BMP4 addition, stained for nuclear marker DAPI and neural ectoderm marker N-cadherin (NCAD). Rightmost set of four images shows the organoids at day 4 of differentiation, stained for DAPI, NCAD, and surface ectoderm marker E-cadherin (ECAD). The six outer organoids show radial patterning, with inward-facing neural ectoderm (NCAD+, green) and outward-facing surface ectoderm (ECAD+, magenta). Scale bar: 300 µm.

(C) Confocal images of a representative single outer organoid after exposure to BMP4 for 1, 2, 3, and 4 days, fixed and stained for DAPI, NCAD, epithelial tight junction marker ZO-1. A neural plate is distinguishable by Day 2, which elevates and folds from a “U” shape to a “C” shape by Day 3 and finally fuses on Day 4 to enclose a single lumen (“O” shape). Scale bar: 100 µm.

(D) Quantification of neural tube closure in Day 4 organoids (left) and relative proportion of ECAD+ surface ectoderm to NCAD+ neural ectoderm (right) as measured by the projected area (see Methods). 100% of the outer organoids across biological replicates (N=4 replicates, n=24 organoids) exhibit closure based on observation of a continuous ring of NCAD (Figure S1F). The ratio of the surface to neural ectoderm has a coefficient of variation (CV) of 0.01.

E) Quantification of lumen (left) and organoid size (right) across Day 4 organoids (n=6), based on the projected area enclosed by NCAD and DAPI staining (Methods). Lumen and organoid size have a CV of 0.09 and 0.04, respectively.

F) Hierarchical clustering in the space of genes identified by sparse multimodal decomposition (Methods) of single-cell RNA-sequencing data from organoids on Day 2 after BMP4 addition (left) and Week 4 human embryos (right), showing the existence of neural ectoderm (green), neural plate border (gray) and surface ectoderm (magenta) cell types in vitro and in vivo marked by similar gene expression patterns. Distinct sets of transcription factors are enriched in neural and surface ectoderm, highlighted in green and magenta respectively, and a subset of these are co-expressed in border ectoderm. Neural ectoderm cells express high levels of forebrain-associated transcription factors (OTX1, OTX2, LHX2, SIX3), and low levels of midbrain/hindbrain-(EN1, EN2) or ventral forebrain-(NKX2-1, GSX2) associated transcription factors. Gene expression (color bar, bottom left) is normalized to median transcript count, followed by log- and min-max normalization.

Identifying transcription factor candidates regulating anterior neurulation.

A) Top left: Illustration denoting the mediolateral axis present in neural tube organoids, spanning the dorsal neural ectoderm, neural plate border, and surface ectoderm. Bottom left: Principal component analysis of single-cell RNA-sequencing (scRNA-seq) data from Day 2 organoids, colored by pseudo-spatial mediolateral coordinate. Right: Standardized mediolateral expression profiles of NCAD, ECAD, and 20 most medially expressed transcription factors from Day 2 organoids. Mediolateral coordinate is denoted by the green-magenta color bar (top) and expression is denoted by the blue-yellow color bar (bottom left).

B) Day 4 organoids transduced with a dual-guide RNA targeting ZIC2. Top left: Epifluorescence merged-channel image of organoids on the pattern stained for NCAD and ECAD. Scale bar: 300 µm. Right: Confocal image of one organoid stained for NCAD, ECAD, and DAPI. Scale bar: 100 µm. Bottom left: Quantification of neural tube closure in Day 4 control and ZIC2 knockdown (KD) organoids, showing that 0% of the control and 100% of ZIC2 knockdown organoids possess an open neural plate.

C) Top 40 up- and down-regulated genes in the neural ectoderm (NE) of Day 4 organoids upon ZIC2 knockdown as compared to scramble control, based on log2(fold change) (color bar, bottom). Transcription factors highlighted in bold.

D) Histogram of enrichment scores based on MEME SEA (Methods) of most overrepresented transcription factor family motifs in regulatory regions 10kb upstream and 100bp downstream of transcriptional start sites for the most down-regulated genes upon ZIC2 knockdown, using a threshold of log2(mean fold change relative to scramble control) < −0.4. Families are color-coded if a transcription factor from that family is present in our candidate selection (Figure 2E).

E) Plot of mean expression in neural ectoderm (transcripts per million, TPM) versus log2(fold change) in their expression between the neural and surface ectoderm in Day 2 organoids for all transcription factors (n=1482). Top 20 most medially expressed transcription factors (black or black outline) and additional candidates (gray, n=58) selected for knockdown (n=78, total). Data points corresponding to members of TF family with enrichment of binding sites near ZIC2-regulated genes (D, Methods) are colored by family membership. All 78 candidates are listed in Figure S2F.

High-throughput method for arrayed single-gene perturbations.

(A) Top: Schematic of pooled perturbations (as generated in previous work) that result in a mosaic organoid, in which a distinct sub-population of cells are perturbed by each guide RNA (shown in different colors). These mosaic organoids can be dissociated for single-cell transcriptomics. Bottom: Schematic of arrayed perturbations (generated in this work) that result in homogeneously perturbed organoids: all cells in each set of organoids are perturbed by the same guide (distinct guides shown in different colors). These organoids can be screened for morphological defects and subsequently dissociated and pooled for single cell transcriptomics.

(B) Efficiency measurements of transduction protocols. Representative images (left) and quantifications (right) correspond to the following protocols. Protocol 1: Virus added 24 hours after seeding and incubated for 24 hours. Protocol 2: Virus added 24 hours after seeding and incubated for 1 hour. Protocol 3: Virus added 1 hour after seeding and incubated for 24 hours. Protocol 4: Virus added 1 hour after seeding and incubated for 1 hour. Protocol 5: Virus added during seeding and incubated for 1 hour. See Figure S3D for transduction efficiency calculation.

(C) Correlation measurement of transduction and knockdown efficiencies. Left: Representative images of circular micropatterns show increasing transduction efficiency (mCerulean, mCer) and decreasing OCT4 expression with increasing viral volume of virus (volumes displayed left of each row) carrying a transfer plasmid with dual guides targeting OCT4, added to 150,000 hPSCs during seeding. Scale bar: 300 µm. Middle: Quantification of nuclear OCT4 fluorescence of individual cells in each condition. Dotted line (threshold) indicates 95th percentile value of OCT4 fluorescence in mCer-cells. Right top inset in each plot is the percentage of mCer+ cells with OCT4 fluorescence above the threshold, indicating OCT4 retention. Right: Plot of percentage of cells showing OCT4 knockdown vs. percentage of cells that express mCer expressed by the guide plasmid show a positive correlation (R2=0.99) demonstrating that live mCer fluorescence coverage of a micropattern is an accurate indicator of percentage of cells transduced.

(D) Schematic of small-volume viral growth and 24-well organoid chip seeding. Dual gRNA plasmids are first arrayed in a 96-well plate. Different colors represent unique lentiviral transfer plasmids targeting different genes (Figure S2B). HEK cells are transfected in a 96-well plate to form an arrayed viral library, then applied to a glass coverslip containing micropatterned hPSCs (depicted with blue dots) using a temporary well system (depicted in dark gray). The well system is then removed, resulting in a glass chip with 24 sets of differently perturbed organoids, each with multiple replicates, that can be cultured in a single media condition.

(E) Left: Quantifications of mCer fluorescence coverage of micropatterns in a mock-transduced well without virus (Well 1, 5% median mCer fluorescence coverage, n=16 micropatterns per well) and 23 scramble vector-transduced wells (Wells 2-24, 99% median mCer fluorescence coverage across wells, n=16 micropatterns per well). See Figure S3D for calculation. Right: Stitched images of viral transduction of micropatterns in Wells 1-6. Transduction efficiency across wells is high and consistent, with a low cross-contamination rate. Scale bar: 300 µm.

SOX11, ZIC2 and ZNF521 are necessary for neural tube closure.

(A) Phenotypic scores of organoids knocked down with each gene. Scoring was based on whether organoids had no closure defects (score = 0, white), a minor defect (score = 1, gray), or a major defect (score = 2, black) consisting of either a fully opened neural plate or multiple closure points.

(B) Epifluorescence microscopy of organoids at Day 4 (unless otherwise noted), showing examples of no defect (top row, scrambled control), minor defect scores (HESX1 knockdown, row 2), or major defect scores (ZIC2, SOX11 and ZNF521 knockdown, rows 3-5 respectively). Three biological replicates per condition are shown with n=6 outer organoids per replicate. Imaging of HESX1 knockdown organoids with a minor defect on Day 5 show that most of them eventually close and that minor openings represent delayed closure. Scale bar: 300 µm.

(C) Corresponding confocal images of single organoids for each of the phenotype classes in (B), showing a range of closure phenotypes. Scale bar: 100 µm.

SOX11 and ZIC2 co-regulate a shared set of genes in opposition to ZNF521.

(A) Fold changes of the 40 most up- or down-regulated genes in the neural ectoderm of each knockdown condition at Day 4, as compared to organoids transduced with a scramble control guide. Genes that are co-downregulated by SOX11 and ZIC2 knockdown are shown in bold.

(B) Cluster dendrogram of WGCNA analysis of scRNA-seq gene expression from neural cells (scramble control) showing modules of co-expressed genes. Module assignment (numbered) is indicated by the color bar; genes labeled in gray were not assigned.

(C) Module eigengene values (representative first principal component of expression of genes in a module) of Module 9 per cell in the neural cell population in each knockdown condition compared to scramble control data. SOX11 knockdown leads to significant downregulation of genes in Module 9 (p < 0.0001), which includes ZIC2 (Fig. S2D).

(D) Scatterplots showing fold changes of significantly co-downregulated genes (Anderson Darling test, p<0.05) in the neural cell populations from transcription factor knockdown pairs, as compared to organoids transduced with a scramble control guide. Genes are plotted if their fold change expression (FC) relative to scrambled control passes the threshold log2(FC) <-0.4 in both knockdown conditions. Genes are labeled if log2(FC) <-0.6 in both conditions.

(E) Scatterplots showing fold changes of genes that are significantly downregulated (Anderson Darling test, p<0.05) by ZIC2 or SOX11 knockdown and significantly upregulated by ZNF521 knockdown, as compared to organoids transduced with a scramble control guide.

(F) Gene regulatory model: SOX11 and ZIC2 are upstream a similar set of genes, including neural tube defect candidates from the literature like PAX2, HESX1, and CRABP1. ZNF521 upregulates a separate set of genes but downregulates several shared ZIC2- and SOX11-upregulated genes.

Development and validation of a hPSC-derived in vitro model of anterior neurulation.

A) Epifluorescence images of fixed and stained hPSCs 2 days after seeding (Day −2 according to scheme in S1D) onto hexagonally arranged 250-µm diameter micropatterns with 100-µm spacing (Methods), showing co-expression of pluripotency markers SOX2 and OCT4. Scale bar: 300 µm.

B) Live epifluorescence images of neural tube organoids expressing an OTX2-Citrine reporter to evaluate responses to different BMP4 treatment courses. The bottom right condition (Day 0-4 BMP4, 10 ng/ml) was selected for further experiments. Scale bar: 300 µm.

C) Live epifluorescence images of neural tube organoids expressing an OTX2-Citrine reporter to evaluate responses to different WNT activation and inhibition courses with agonist CHIR99021 (CHIR) and inhibitor IWR-1-endo (IWR-1). The top right condition (Day 0-4 IWR-1, 1 µM) was selected for further experiments. Scale bar: 300 µm.

D) Neural tube organoid differentiation protocol. hPSCs are seeded onto 2D micropatterns on Day −4 (2 days prior to differentiation) in mTeSR Plus (MT) with 10 µM ROCK inhibitor Y-27632 (Y). Media is changed to MT without Y on Day −3. On Day −2, organoids are folded into 3D cysts in high-concentration Matrigel and simultaneously differentiated toward ectodermal fates in N2B27 media (Methods) by addition of 10 µM TGFB inhibitor SB431542 (SB), 0.5 µM BMP inhibitor LDN193189 (LDN), and 1 µM WNT inhibitor IWP-3 on Day −2. On Day 0, addition of 10 ng/ml BMP4 induces mediolateral patterning of the ectoderm, while WNT inhibition is sustained by 1 µM IWR-1 addition. BMP4-containing media is replaced on Day 2. Each media replacement includes fresh high-concentration Matrigel.

E) Epifluorescence images of organoids after 1, 2, and 3 days of BMP4 treatment, fixed and stained for neural ectoderm marker NCAD, tight junction marker ZO1, and nuclear marker DAPI. Scale bar: 300 µm.

F) Epifluorescence images of three biological replicates of neural tube organoids on Day 4, showing complete closure in the form of a continuous ring of NCAD. 100% of organoids show closure (see Figure 1D). Scale bar: 300 µm.

Transcriptome-based identification of candidate transcriptional regulators of anterior neurulation.

A) Principal component analysis of single-cell RNA-sequencing data from Day 2 organoids, colored by expression of 20 most medially expressed transcription factors. Gene expression is transcript count-normalized, and units are in transcripts per million. Scale bar for each sub-panel is shown to the right of the PCA plot.

B) Cloning schematic for lentiviral dual-guide RNA vector (based on Replogle, et al., 2020), for which pLKO5 backbone was modified to carry an mCerulean (mCer, a cyan fluorescent protein) reporter, digested with BsmBI, and ligated with two pairs of annealed sgRNA oligos together with a BsmBI-digested gBlock. The final vector encodes gRNA1-CR3-CS1 under a human U6 (hU6) promoter, gRNA2-CR1-CS2 under a mouse U6 (mU6) promoter, and mCerulean under an EFS promoter (Methods).

C) dCas9-KRAB CRISPRi H1 (WA01) transduced with dual-guide RNAs consisting of either a control non-targeting scrambled sequence (“control”) or OCT4-targeting sequence (“OCT4”) under the hU6 or mU6 promoter in guide RNA positions 1 (“gRNA1”) and 2 (“gRNA2”), respectively, fixed and stained for DAPI and OCT4 48 hours post-transduction. Columns correspond to a control sequence at both positions (first column), an OCT4 sequence in position 1 and control sequence in position 2 (second column), a control sequence in position 1 and OCT4 sequence in position 2 (third column), and an OCT4 sequence at both positions (fourth column). Top and bottom rows correspond to merged channel images of the same fields of cells stained for OCT4 (red) and DAPI (gray), with or without the mCerulean channel (cyan) overlaid. The cells transduced with virus (mCer+) in the second, third, and fourth columns, but not those the first column, have low OCT4 levels, demonstrating that the OCT4-targeting guide is specific and functional when present in the first, second, or both positions. Scale bar: 50 µm.

D) scRNA-seq of Day 4 organoids transduced with a control non-targeting scramble dual-guide RNA construct. Control cells were clustered as neural ectoderm (NE, green), neural plate border (NPB, gray), or surface ectoderm (SE, magenta) based on expression of labeled transcription factors. Gene expression (color bar, bottom left) is normalized to median transcript count, followed by log- and min-max normalization.

E) Left: Cells from scRNA-seq of Day 4 scramble control and ZIC2 knockdown organoids plotted along the first two principal components for control cells. The data points corresponding to the three cell types are colored appropriately (NE, green; NPB, gray; SE, magenta). Right: Histogram of cell type proportions shows that ZIC2 knockdown organoids have less NE, more NPB, and a similar SE fraction compared to control organoids.

F) Full list of 78 transcription factor candidates selected for knockdown, with 20 neural transcription factors selected based on mediolateral expression (Figure 2A, S2A) in bold.

G) Standardized mediolateral expression profiles of top 200 most highly expressed transcription factors in Day 2 organoids, ordered from most medially expressed (top) to most laterally expressed (bottom). The 20 most medially expressed transcription factors (Figure 2A, right) are separated above the rest.

H) Plot of mean expression in neural ectoderm (transcripts per million, TPM) versus log fold change in expression between the neural and surface ectoderm in Day 2 organoids for the top 20 mediolaterally enriched transcription factors from Figure 2E. For all 20 transcription factors, mean expression > 70 TPM and log2(fold change) > 0.6.

Development of transfection and transduction methods.

(A) Gel electrophoresis of non-clonal plasmid stocks (see Figure S2B, Methods) after double digest with BamHI and DraIII. Expected band length for backbone (control lane): 7.4kb. Expected band lengths for backbone with ligated inserts (following lanes): 3.4kb and 4.5kb. Ladder: NEB 1kB Plus. Band intensities show that a majority of non-clonal plasmid stock contains the desired insert.

(B) Top: Epifluoresence images of mCerulean (mCer) expression in confluent hESC culture, 48 hours after lentiviral transduction with scramble control vector from one of 16 different HEK transfection protocols. Protocols altered 4 variables in a 96-well transfection of HEK cells. Bottom: Ratio of mCer to DAPI signal from each protocol above, rank ordered by transduction efficiency. Protocol 7 (*) is the standard protocol recommended by jetPrime. Scale bar: 300 µm.

(C) Ratio of mCer to DAPI signal from Figure S3B, paired by shared variables to show the effect of changing the variable in each controlled pair. Changes that resulted in increased transduction are shown in green. General trends show that decreasing viral growth volume and decreasing plasmid transfection volume result in higher viral production.

(D) Example of transduction efficiency analysis. Left: Colony diameter as measured by colony segmentation by Ilastik in the phase channel. Middle: Colony segmentation as applied to the mCer channel (white) and 2D background halo (area between gray circles). Right: Quantification of mCer signal from colony (white) and background halo (gray). Dotted line (threshold) represents 95th percentiles of mCer expression in the background halo. This approach is validated with nuclear analysis in Figure 3C. Scale bar: 300 µm.

(E) mCer expression in hESC colonies, 48 hours after transduction with scramble control vector made with Protocol 4 from Figure S3B. Top: hESCs were transduced after 48 hours of adherent growth on plates. Bottom: hESCs were transduced in single cell suspension directly after dissociation. Scale bar: 300 µm.

Perturbation of most transcription factors does not give rise to morphological phenotypes.

(A) 24 examples of live epifluorescence microscopy of lentiviral marker mCerulean (mCer), 48 hours after transduction in suspension and seeding on micropatterns as in Figure 3D. Scale bar: 300 µm.

(B) Quantification of transduction efficiencies from individual sgRNAs, 48 hours after transduction in suspension and seeding on micropatterns as in Figure 3D. The median transduction efficiency across all guides using the non-purified plasmid prep was 93%.

(C) Top: Epifluorescence image of RFX7 perturbation. Bottom: Maximum intensity projection of a confocal Z-stack (bottom). All organoids received a score of 0 (no phenotype) after confocal imaging. Scale bar: 300 µm

(D, E) Epifluorescence microscopy of terminal organoid morphologies, differentiated as in Figure S1D for 6 days, then fixed and immunostained. E: Examples of minor openings are indicated in white arrows. Scale bar: 300 µm.

(F) Distribution of mCer in organoids transduced with a scramble guide, measured 8 days post-transduction (2 days of mTeSR growth and 6 days of differentiation). Scale bar: 300 µm.

Analysis of scRNA-seq data from perturbed organoids.

(A) Principal component analysis of cell types in Day 4 organoids infected guides against SOX11 or ZNF521. Cells are plotted using the two major principal components from the scramble control data from Figure S2E. Cell type assignments to neural ectoderm (green), neural plate border (gray), and surface ectoderm (magenta) are based on expression of canonical transcription factors from scramble control data in Figure S2D.

(B) Quantification of cell type proportions from scRNA-seq of all control and major closure defect datasets (ZIC2, SOX11, and ZNF521 knockdown). Cell types are assigned as in (A).

(C) Top: Knockdown efficacy calculations of sgRNAs that led to no phenotype or minor closure defects. Seven out of eight guides showed significant knockdown of their target gene with a knockdown efficacy of >40%. Bottom: Knockdown efficacy calculation of sgRNAs that produced a major closure defect. All three guides showed significant knockdown of their target gene with a knockdown efficacy of >40%.

(D) Dendrogram module eigengene values from the remaining modules from Figure 5C. No other module shows significant changes in expression in any perturbation condition tested.

(E) Correlation values with module eigengene (module membership) of top 20 genes in Module 9, including ZIC2, from Figure 5C.

(F) Epifluorescence images of Day 4 organoids perturbed with guides against downstream targets of SOX11, ZIC2, and/or ZNF521. Scale bar: 300 µm