ResourceEpigenomic Landscapes of hESC-Derived Neural Rosettes: Modeling Neural Tube Formation and DiseasesGraphical AbstractHighlightsd Neural rosettes have transcriptional and epigenomic landscapes distinct from 2D NSCs d Lowmethylated regions in NRs predict enhancers active later in neurodevelopment d NR regulatory maps provide insights into mechanisms of neural tube defectsValensisi et al., 2017, Cell Reports 20, 1448–1462 August 8, 2017 ª 2017 The Author(s). http://dx.doi.org/10.1016/j.celrep.2017.07.036Authors Cristina Valensisi, Colin Andrus, Sam Buckberry, Naresh Doni Jayavelu, Riikka J. Lund, Ryan Lister, R. David Hawkins Correspondence rdhawk@uw.edu In Brief Neural rosettes offer a 3Dmodel of neural tube development. Valensisi et al. show that these cells have distinct epigenomic landscapes compared to 2D neural stem cells. cis-regulatory elements identified in NRs harbor folate-associated CpGs, which may be important during neural tube development, and neurodevelopmental-disorder- associated variants. Cell Reports ResourceEpigenomic Landscapes of hESC-Derived Neural Rosettes: Modeling Neural Tube Formation and Diseases Cristina Valensisi,1,2 Colin Andrus,1,2 Sam Buckberry,3,4 Naresh Doni Jayavelu,1,2,5 Riikka J. Lund,5,6 Ryan Lister,3,4 and R. David Hawkins1,2,5,7,* 1Division of Medical Genetics, Department of Medicine and Department of Genome Sciences, University of Washington School of Medicine, Seattle, WA, USA 2Institute for Stem Cell and Regenerative Medicine, University of Washington School of Medicine, Seattle, WA, USA 3Australian Research Council Centre of Excellence in Plant Energy Biology, School of Molecular Sciences, The University of Western Australia, Perth, WA, Australia 4Harry Perkins Institute of Medical Research, Perth, WA, Australia 5Turku Centre for Biotechnology, University of Turku, Turku, Finland 6A˚bo Akademi University, Turku, Finland 7Lead Contact *Correspondence: rdhawk@uw.edu http://dx.doi.org/10.1016/j.celrep.2017.07.036SUMMARY We currently lack a comprehensive understanding of the mechanisms underlying neural tube formation and their contributions to neural tube defects (NTDs). Developing a model to study such a complex morphogenetic process, especially one that models human-specific aspects, is critical. Three-dimen- sional, human embryonic stem cell (hESC)-derived neural rosettes (NRs) provide a powerful resource for in vitro modeling of human neural tube formation. Epigenomic maps reveal enhancer elements unique to NRs relative to 2D systems. A master regulatory network illustrates that key NR properties are related to their epigenomic landscapes. We found that folate-associated DNA methylation changes were enriched within NR regulatory elements near genes involved in neural tube formation and metabolism. Our comprehensive regulatory maps offer insights into the mechanisms by which folate may prevent NTDs. Lastly, our distal regulatory maps provide a better understanding of the potential role of neuro- logical-disorder-associated SNPs. INTRODUCTION Advances in 3D culture systems have opened the door to un- precedented opportunities for recapitulating the in vivo features of human CNS development (Kelava and Lancaster, 2016). Neural rosettes (NRs) derived from human embryonic stem cells (hESCs) can recapitulate the molecular and morphoge- netic sequence of events from gastrulation to neural tube (NT) formation (Deglincerti et al., 2016; Kelava and Lancaster, 2016; Zhang et al., 2001). These 3D tubular structures emerge from embryoid bodies (EBs) with a differentiation timing com-1448 Cell Reports 20, 1448–1462, August 8, 2017 ª 2017 The Author This is an open access article under the CC BY-NC-ND license (http://parable to that seen in vivo (Pankratz et al., 2007). NRs display apical-basal polarity typical of neuroepithelium and radial orga- nization similar to that of radial glial cells, a primary neural pro- genitor cell population present during development and in the adult brain (Conti and Cattaneo, 2010). Previous studies have explored molecular markers and signaling pathways establish- ing the identity and differentiation potential of NRs and showed that NRs better recapitulate the in vivo properties of bona fide neural stem cells (NSCs) (Elkabetz et al., 2008; Koch et al., 2009). Recently, 3D culture systems were used to investigate molecular mechanisms and signaling pathways governing self-organization of NT formation and NT closure in both mouse and non-human primates (Meinhardt et al., 2014; Zhu et al., 2016). These studies serve as proof of concept of the power of 3D culture systems for modeling NT formation and its disruption. Disruption in NT formation and closure leads to NT defects (NTDs) (Copp and Greene, 2013), the second most common birth defects and a global public health burden (Zaganjor et al., 2016). An important player in NTD etiology is folate (Nazki et al., 2014), a critical component in the one-carbon metabolism pathway providing methyl groups for a range of biochemical re- actions, including methylation of DNA and histones. How folate affects NT formation is still poorly understood, especially for human-specific aspects (Wallingford et al., 2013), but major evidence points to changes in DNA methylation and gene expression (Copp and Greene, 2013). Epigenomic remodeling plays a major role in development and cell fate determination, and alterations to DNA methyl- ation and histone modifications, here referred to as epige- netics, have been associated with disease (Romanoski et al., 2015). To date, we are still missing a comprehensive under- standing of epigenetic mechanisms for NT formation and NTDs. Recently, epigenomic dynamics of neural differentiation from hESCs were investigated using monolayer culture sys- tems, namely 2D NSCs (Xie et al., 2013; Ziller et al., 2015). However, 2D culture systems lack the complexity to(s). creativecommons.org/licenses/by-nc-nd/4.0/). A B C Figure 1. Distinct Neural Rosette Transcriptional Profiles Compared to 2D Neural Stem Cells (A) Top: schematic of the neural rosette differentiation protocol. Bottom: immunostaining for neural markers in NRs at day 11 of differentiation. (B) UCSC Genome Browser snapshot of the PAX6 locus (chr11:31,800,000–31,850,000) showing normalized signals of DNA methylation level (mCG/CG), RNA- seq (read count expressed as TPM), and ChIP-seq reads (RPKM, input normalized). ChIP-seq peak calls are shown as bars above each track. (C) Heatmap representation of expression levels (TPM) of entropy-based cell-type-specific protein-coding (left) and lncRNA (right) genes. See also Supplemental Experimental Procedures, Figure S1, and Table S1.recapitulate morphogenetic properties of the NT, a key aspect to understand how its disruption leads to NTDs. Here, we leverage the 3D nature of hESC-derived NRs to explore the epigenomic regulation of NT formation and possible applica- tions for studying NTDs and other neurological disorders con- nected with early development. RESULTS Distinct NR Transcriptional Profiles Compared to 2D NSCs NRswere generated from hESCs via formation of EBs (Figure 1A; Supplemental Experimental Procedures). To gain insight into the role of epigenomic regulation in defining NR identity and NT for- mation, we generated comprehensive maps of gene expression by strand-specific RNA sequencing (RNA-seq); histone H3mod- ifications via chromatin immunoprecipitation sequencing (ChIP- seq) that are commonly used to predict promoter and enhancerelements and their functional status, including lysine 4 trimethy- lation (H3K4me3) and monomethylation (H3K4me1) and lysine 27 acetylation (H3K27ac) and trimethylation (H3K27me3); and genome-wide DNA methylation by whole-genome bisulfite sequencing (WGBS) (Figure 1B). To explore molecular features that set apart NRs from 2D NSCs, we compared the NR tran- scriptome and epigenome to publicly available data for three 2D NSCs that match the developmental window of NRs, namely neural progenitor cells (NPCs) (Xie et al., 2013), neuroepithelia (NE), and early radial glia (ERG) (Ziller et al., 2015). We first focused on the NR transcriptome. Given that all cells will express key genes for NSCs and neurodevelopment, we used an entropy-based method (Schug et al., 2005) to define cell-type-specific gene expression from our global transcrip- tional analysis. We identified 415 protein-coding genes and 122 long noncoding RNAs (lncRNAs) with expression restricted to NRs (Figure 1C; Table S1A). To validate our findings, we also assessed differential expression in a pairwise comparisonCell Reports 20, 1448–1462, August 8, 2017 1449 fashion using rank product (Figure S1A; Table S1B) and observed high concordance with the entropy-based analysis, with 96% of NR-restricted genes recovered in the top 10% for rank product for NR versus all 2D NSCs (Figure S1B; Table S1C). Among NR-restricted protein-coding genes, we found several genes whose function during neurulation can potentially influence NT formation. One such gene was T-box 3 (TBX3), known for its role in maintaining pluripotency in mouse ESCs (Niwa et al., 2009). Interestingly, recent evidence showed that TBX3 knockdown in human ESCs decreases NR formation (Esmailpour and Huang, 2012), suggesting that it may also be important during neurulation. The expression of betaine-homo- cysteine methyltransferase (BHMT), a member of the one-car- bon metabolism pathway, was also found to be restricted to NRs. Genetic variants in BHMT have been reported to influence the risk for spina bifida (Boyles et al., 2006; Morin et al., 2003). Another such gene was ephrin receptor A5 (EPHA5), which is involved in regulation of cell adhesion during cranial develop- ment. EphA5/mice display anencephaly, likely due to the cra- nial neural folds failing to fuse in the dorsal midline (Holmberg et al., 2000). This collection of genes is a strong indication of a NT phenotype. Epigenomic Landscapes at NR Promoters A large fraction of chromatin-enriched regions were unique to NRs with respect to the other cell populations (Figure S2A; Table S2), with the exception of H3K4me3. H3K4me3 marked 13,880 transcription start sites (TSSs; gene-level annotations), com- pared to nearly 20,000 TSSs marked by H3K4me3 in 2D NSCs (Figure S2B). However, NR promoters were not devoid of chromatin dynamics. H3K27ac showed between 39% and 68% of peaks to be NR unique in pairwise comparisons (Fig- ure S2A; Table S2), of which 2,780 were found at promoters (Figure S2C). Gene ontology (GO) term analysis of these puta- tive NR-specific active promoters suggested these genes are involved in general developmental functions (Figure S2D). H3K27me3 showed between 65% and 67% of peaks to be NR unique (Figure S2A; Table S2). We identified 830 bivalent pro- moters (H3K4me3 and H3K27me3) in NRs (6% of all H3K4me3 marked regions), a signature that is a key process during lineage and cell-type specification (Voigt et al., 2013), and we observed that the NR bivalent signature was highly cell-type specific (Fig- ure S2E). GO term analysis of NR bivalent genes showed enrich- ment for processes such as glycoprotein metabolism, metal ion binding, and voltage-gated channel activity (Figure S2F). The poised nature of the promoter chromatin state suggests a de- layed developmental expression of these genes, potentially at a time when distinct neurons are specified. Next, we showed how different classes of H3K4me3-marked promoters, namely active (co-enriched for H3K27ac), bivalent (co-enriched for H3K27me3), and poised (no co-enrichment for H3K27ac or H3K27me3) correlated with the expression level of the genes (Figure 2A). DNA methylation remodeling is an essential component of epigenetic regulation during development (Smith and Meissner, 2013). Recent studies showed that unmethylated and low-meth- ylated regions (UMR/LMRs) can be used to predict regulatory elements (Burger et al., 2013; Stadler et al., 2011). UMRs,1450 Cell Reports 20, 1448–1462, August 8, 2017defined as regions with average methylation lower than 10%, are mostly found at promoters and CpG islands (CGIs). We iden- tified 19,049 UMRs with an average size of 2 kb (Figure S2G) and confirmed their expected genomic locations (Figure S2H). Consistent with their genomic location, a larger fraction of UMRs colocalized with H3K4me3-marked promoters and en- compassed chromatin-marked promoters in different states (Figure 2B). Contrasting UMRs with the chromatin maps, we observed that almost all UMRs were located at chromatin- marker regulatory region (Figure 2C). As for the UMRs that overlapped H3K4me1-marked regions, those are predominantly promoter-proximal enhancers. Consistently with their chro- matin signature, H3K4me3-marked promoters were largely de- methylated, while their genes appeared to be expressed (Fig- ure 2D). This was consistent when we considered methylation at promoters genome-wide, where genes with demethylated promoters were expressed at a high level, while methylation in promoters (>75%) corresponded to genes being undetected or expressed at a low level (transcripts per million [TPM] < 1) (Figure 2E). To compare NRs to 2D NSCmethylomes, we defined differen- tially methylated regions (DMRs) for each pairwise comparison (Figure S2I). When considering regions that were differentially methylated in at least one of the three pairwise comparisons be- tween NRs and 2D NSCs, we found 28,102 DMRs, whereas only 15,769 DMRs were found in pairwise comparisons among 2D NSCs, suggesting that DNA methylation signature in NRs is more different with respect to any 2D NSCs than 2D NSCs are among each other. Furthermore, NRs tend to harbor hypomethy- lated DMRswhen compared to any 2DNSC (Figures 2F and S2L) (25,195 hypomethylated DMRs out of 28,102 DMRs). Of the 1,143 NR-specific promoter DMRs, 459 were hypomethylated, and their genes showed significantly higher expression than genes of hypermethylated promoters (Figure 2G). NRRegulatory Enhancer Networks andDynamics during Nervous System Development Enhancers play a major role in cell-type-specific transcriptional regulation (Shlyueva et al., 2014). The H3K4me1 enhancer signa- ture (Heintzman et al., 2009) was the most cell-type-specific epigenomic modification when compared against all 2D NSC population, with 62% (38,181) peaks unique to NRs (Figure 3A). Co-enrichment of the activating H3K27ac modification with H3K4me1 reveals enhancers in an active state, while elements marked by H3K4me1 alone tend to be in a poised state (Creyghton et al., 2010; Hawkins et al., 2011; Rada-Iglesias et al., 2011). Using these criteria, we identified 22,475 active and 41,797 poised enhancers in NRs. Collectively, these results show that NRs have a very distinct transcriptional and epige- nomic landscape compared to 2D NSCs and highlight a role for distal regulatory elements in defining NR epigenomic signature. To predict the regulatory network shaping the enhancer land- scape in NRs, we performed transcription factor (TF) binding site (TFBS) motif enrichment analysis on both active and poised enhancer regions (Figures 3B and S3A; Tables S3A and S3B) and constructed in silico protein-protein interaction maps for TFs with detectable transcripts in NRs (Figures 3C and S3B). -1kb peak 1kb center U M R (1 9, 04 9) H3K27ac H3K27me3 H3K4me1 H3K4me3 H 3K 4m e3 -m ar ke d pr om ot er s H3K27ac H3K4me3 DNAme Expression -2 14 Methylation Level TPM -1kb peak 1kb center -2kb peak 2kb center 0 18 0 1 RPKM RPKM 0 50 1 0.5 0 M ethylation Level NR NPC ERG NE 33.9% 15.3% 13.2% 5.3% 8.0% 24.2% 0.5 0 D ensity E xp re ss io n (lo g2 (T P M )) Methylation Level 8 4 0 0 0.25 0.75 1 h h E xp re ss io n (lo g2 (T P M )) 10 5 0 Hyper DMR Hypo DMR ** active bivalent poised 9,253 5,311 793 2,840 3,799 798 active promoter (K4me3/K27ac) poised promoter (K4me3) bivalent promoter (K4me3/K27me3) active enhancer (K4me1/K27ac) poised enhancer (K4me1) no overlap E xp re ss io n (lo g2 (T P M )) 10 5 0 *** *** * A F G E D C B Figure 2. Epigenomic Landscapes at Neural Rosette Promoters (A) Violin plot of distribution of the expression values (log2(TPM)) of genes for active (H3K4me3/H3K27ac), poised (H3K4me3 only), and bivalent (H3K4me3/ H3K27me3) promoters. ***p < 0.0001, *p% 0.01 (t test). (B) Pie chart showing overlap of UMRs with chromatin-marked promoters and enhancers. (C) Heatmap representation of ChIP-seq signal intensity (RPKM, input normalized) for H3K4me1, H3K4me3, H3K27me3, and H3K27ac over 2-kb windows centered on UMRs. (D) Heatmap representation of ChIP-seq signal intensity (RPKM, input normalized) for H3K4me3 and H3K27ac over 1-kb windows centered on TSSs, DNA methylation levels over 2-kb windows centered on TSSs, and gene expression (TPM) for H3K4me3-marked promoters. (E) Density plot showing gene expression (log2(TPM)) and DNA methylation levels for all promoters (defined as ±1 kb around TSSs; Gencode 19). (F) Heatmap and clustering visualization of aggregate mC levels (mCG/CG) in DMRs (delta-mC > 0.2) from all comparisons for positions with >53 coverage in all cell types. 31,395 DMRs are shown in the heatmap. (G) Violin plot of distribution of the expression values (log2(TPM)) of genes which promoters harbor a DMR. **p < 0.001 (t test). See also Figure S2 and Table S2. Cell Reports 20, 1448–1462, August 8, 2017 1451 AC D E B Figure 3. NR Regulatory Enhancer Networks and Dynamics during Nervous System Development (A) Venn diagram shows the overlap for all H3K4me1 peaks among NRs, NPCs, NE, and ERGs. Fraction of the cell-type-unique peaks is also given as percentage of the total number of peaks. (legend continued on next page) 1452 Cell Reports 20, 1448–1462, August 8, 2017 Our analysis suggests that SOX2 plays a central role in TF networks anchored at either active or poised enhancers. SOX2 is expressed throughout development and is involved in main- taining pluripotency as well as differentiation of several lineages (Sarkar and Hochedlinger, 2013). In NSCs, SOX2 is critical for maintaining neural progenitor identity by inhibiting premature dif- ferentiation toward the neuronal lineage (Graham et al., 2003; Lee et al., 2014). Another major hub in both networks is TCF3. TCF3 is a key player in controlling the balance between differen- tiation and maintenance of neural progenitors (Kuwahara et al., 2014) by antagonizing the Wnt/b-catenin signaling pathway (At- lasi et al., 2013). TCF3 is also necessary for establishing the ante- rior-posterior axis (Merrill et al., 2004). Among TFmotifs uniquely enriched at active enhancers, we observed E2F4, a member of a family of TFs involved in regulating the cell cycle. E2f4/ mice have brain malformations and abnormal fetal expression of key neural markers such as Pax6, Nkx2.1, and Dlx (Swiss and Ca- saccia, 2010). Among TF motifs uniquely enriched at poised en- hancers is the transcriptional repressor GLI3, a mediator of the Shh signaling pathway. Gli3 mediates the graded Shh signal (Stamataki et al., 2005) and is essential for establishing dorsal- ventral patterning in the NT (Persson et al., 2002). Given GLI3 is expressed in NRs, its enrichment at poised enhancers may shed light on the dynamic state of these elements. Next, we asked how the NR enhancer landscape might change during nervous system development. To address this question, we obtained DNase I hypersensitivity site (DHS) data for human fetal spinal cord (15 weeks of gestation) and brain (18 weeks of gestation) from the Roadmap Epigenome Project (Thurman et al., 2012). We found that 13% of both active and poised enhancers are hypersensitive at these later stages in neu- rodevelopment, suggesting they remain or become active (Fig- ure 3D, cluster 1). We assigned hypersensitive enhancers to the nearest neighboring genes (NNGs) (Tables S3C and S3D) and performed functional annotation clustering analysis of associated GO terms (Figure 3E). For those active NR enhancers that are hypersensitive during brain and spinal cord develop- ment, top enriched clusters suggested these genes may play a role in more general neurodevelopmental processes, such as TF activity and neuron differentiation. Top GO term clusters for poised NR enhancers that potentially become active as the brain and spinal cord form were instead enriched for terms linked to synapse signaling and axonogenesis. These latter functions, associated with more mature neuronal activities, are consistent with the poised for activation state of the enhancers. To validate these findings, we also used an alternative method for assigning regulatory elements to putative gene targets, GREAT (McLean et al., 2010), and observed more than 72%(B) Heatmap representation of number (expressed as percentage of the total) of e motifs whose TF was expressed (TPM > 1) in NRs are shown. See Table S3A for (C) In silico protein-protein interaction maps for TFs whose motifs are enriched a number of interactions; color represents expression value (TPM) of the gene as m (D) Distribution of the average enrichment of DHS signals (RPKM) over 6-kb windo (K = 2) was used to identify enhancers enriched for DHSs (cluster 1, in blue). Clu (E) Bar plots show the most significantly enriched clusters of GO terms associated poised enhancers, respectively. See also Figure S3 and Table S3.concordance with NNGs (Figure S3C); GO term analysis of GREAT-assigned putative targets showed high similarity as well (Figure S3D). Broad Distal Regulatory Domains Reveal Master Regulators of NT Formation Recently, a great deal of interest has emerged around broad distal regulatory domains—comprising sites stitched together, namely super-enhancers (Hnisz et al., 2013; Whyte et al., 2013) or stretch enhancers—defined largely by broad acetyla- tion peaks (Parker et al., 2013). Besides sharing all the prop- erties of average sized enhancers (commonly considered less than 1 kb), super-enhancers are often located in close proximity to known master regulators (e.g., pluripotent genes and oncogenes) and are deemed to act as switches to deter- mine cell fate. To identify broad enhancers and their putative target genes based on genomic proximity, we first defined broad H3K4me1 and H3K27ac domains as peaks longer than 3 kb and identified 4,245 (7%) and 1,964 (6%) peaks, respectively (Figure S4A). While most broad H3K27ac peaks were co-enriched for H3K4me1, indicating enhancers in an active state, a large fraction of broad H3K4me1 domains ap- peared to exist in a poised state (Figure 4A). Consistent with their functional state, active broad enhancers were largely found near highly expressed genes, while expression values of nearest genes to poised broad enhancers were significantly lower (Figure 4B). All expressed genes were compared as a control. Among putative targets of active broad enhancers were SKI and HES5 (Figure 4C; see Tables S4A and S4B for com- plete list). HES5 is a primary effector of the Notch signaling pathway and plays a critical role in NT development (Louvi and Artavanis-Tsakonas, 2006; Yoon and Gaiano, 2005). In NRs, HES5 prevents the switch from neural to neuronal iden- tity (Abranches et al., 2009). Similarly, SKI functions as a repressor of transforming growth factor (TGF) beta signaling and may play a role in neurulation and patterning of the NT (Berk et al., 1997). Overall, GO term analysis showed that pu- tative targets are primarily involved in processes specific to NT formation, such as establishment of rostral-caudal axis and regionalization of the nervous system (Figure 4D). Using GREAT to assign broad enhancers to putative gene targets, we observed over 75% concordance with NNG method (Figure S4B). Next, we explored which TFs were more likely to be found at broad enhancers than at standard enhancers. We performed TFBS motif enrichment analysis of both active and poised broad enhancers using standard peaks as background (Figure 4E;nhancer regions were a motif was significantly enriched (q-value < 0.001). Only complete list. t active (left) and poised (right) enhancers. Size of node is proportional to the easured by RNA-seq. ws centered on active (top) and poised (bottom) enhancers. K-mean clustering ster 2, in green, represents enhancers lacking enrichment for DHSs. with NNGs in cluster 1 of (D). See Tables S3C and S3D for NNGs of active and Cell Reports 20, 1448–1462, August 8, 2017 1453 12 8 4 0 -5kb peak 5kb center R P K M R P K M br oa d H 3K 27 ac (1 ,9 65 ) br oa d H 3K 4m e1 (4 ,2 45 ) H3K4me1 H3K27ac 12 8 4 0 H3K27ac H3K4me1 all genes 9 6 3 0 E xp re ss io n (lo g2 (T P M )) -6 -4 -2 0 nervous system development bone development eye development pattern specification process limb morphogenesis forebrain development neural tube development regionalization neurogenesis neg. regulation of gliogenesis DLL1 SKI FOXP1 SMAD3 FZD3 SOX11 GLI2 SOX2 HES5 WNT1 ID1 ZIC5 LHX2 FZD5 PAX6 HES3 PBX2 IRX2 RFX3 NOTCH1 SOX1 OTX1 WNT7B OTX2 log10 enrichment score H3K4me1 H3K27ac SKI 50 kb HES5 A D E C B *** *** Fo ld E nr ic hm en t 0.0 1.0 2.0 3.0 4.0 5.0 FO XA 1 HO XA 2 HO XA 9 HO XB 4 LH X2 LH X3 PA X7 RF X2 SO X2 broad H3K4me1 broad H3K27ac Figure 4. Broad Distal Regulatory Domains Reveal Master Regulators of Neural Tube Formation (A) Heatmap representation of the ChIP-seq signal intensity (RPKM, input normalized) over a 10-kb window centered on H3K4me1 and H3K27ac peaks across broad H3K27ac peaks (top) and broad H3K4me1 peaks (bottom). (B) Violin plot of distribution of the expression values (log2(TPM)) of the NNGs for broad H3K4me1 (blue) and H3K27ac (orange) peaks. Expression of all genes detected by RNA-seq in NRs is shown as a control. ***p < 0.0001 (t test). (C) UCSC Genome Browser snapshots of SKI and HES5 loci. ChIP-seq signals (RPKM, input normalized) for H3K4me1 (blue) and H3K27ac (orange) are dis- played. Solid bars indicate peaks called by MACSv1.4. (D) Bar plot shows the most significantly enriched GO terms associated with NNGs of broad enhancers. (E) Bar plot of TFBS motifs significantly enriched (q-value < 0.001) at broad enhancers versus standard enhancers. Enrichment at broad H3K4me1 (blue) and H3K27ac (orange) regions is expressed as ratio between fraction of broad regions and fraction of standard regions enriched for the given TFBS motif. See also Figure S4 and Table S4.Tables S4C and S4D). We found that LHX2 and HOXA9 motifs are enriched at both broad enhancer classes. Both factors are known to co-localized with Sox2 during early neurodevelopment in mice and are critical determinants of regional-specific Sox2 occupancy. For example, Lhx2 is more enriched in rostral regions (anterior fate), and HoxA9 is more enriched in caudal re-1454 Cell Reports 20, 1448–1462, August 8, 2017gions (posterior fate) (Hagey et al., 2016). While SOX2 and LHX2 genes are highly expressed in NRs, HOXA9 was not detected in NRs. The co-presence of motifs for both anterior and posterior factors may contribute to a key feature of NRs, which is a largely anterior identity while still possessing the capability of res- ponding to posterior patterning cues (Conti et al., 2005; Elkabetz et al., 2008; Pankratz et al., 2007). Taken together, these results depict the use of broad enhancers near master regulators that underlie key properties of NRs, such as dependency on Notch signaling for maintaining NSC identity and unrestricted differen- tiative potentials. At the same time, NRs, like the NT, may be primed to respond to imminent regionalization signaling by the likely accessibility of key TFs, such as HOXA9 at broad enhancers. LMRs Predict Distal Regulatory Elements and Their Developmental Dynamics LMRs have an average methylation ranging from 10% to 50%; are regions of low CG density; tend to be enriched for H3K4me1, DHSs, and p300/CBP; and are primarily located distal to promoters in intergenic or intronic regions (Burger et al., 2013; Stadler et al., 2011). We identified 56,558 LMRs (500 bp average size) in the NR methylome (Figure S2G) and determined their genomic annotations (Figure S5A). We con- trasted LMRs with NR-specific hypomethylated DMRs. We found that NR hypomethylated DMRs largely occur at LMRs (48% of 2,995 hypo-DMRs common across all three pairwise comparisons against 2D NSCs; Figure S5B), whereas only a small fraction (6%) were found to be located at UMRs. This was consistent with the observation that NR hypomethylated DMRs largely occur in intergenic and intronic regions (Fig- ure S5C) and are enriched for H3K4me1 and, to a lesser extent, for H3K27ac, but not for H3K4me3 (Figure S5D). We then inte- grated all LMRs with our chromatin maps (Figures 5A and 5B). As expected based on their genomic localization, LMRs were depleted of the promoter mark H3K4me3 but enriched for the enhancer mark H3K4me1. The larger fraction of LMRs was devoid of chromatin marks. It has been suggested that TF binding may be both necessary and sufficient to induce local demethylation and establish LMRs (Sta- dler et al., 2011). We asked whether unmarked LMRs might represent distal regulatory regions that are poised to become active at a later stage in development but have yet to acquire indicative chromatin marks. To explore this hypothesis, we first sought to explore whether these elements might be able to bind TFs relevant for neural differentiation. Thus, we performed TFBS motif enrichment analysis on the unmarked LMRs (Table S5A) and compared to active and poised enhancers (Tables S3A and S3B). We found that unmarked LMRs were enriched for a similar set of motifs to those found at enhancers (Figure 5C; Table S5B), with 73% of motifs shared with enhancers (Fig- ure S5E), suggesting that these regions have the ability to func- tion as distal regulatory elements. Among the motifs most enriched in unmarked LMRs but absent in H3K4me1-marked enhancers was the motif for the repressor REST, which is known for its role in maintaining NSCs and preventing premature neuronal differentiation (Huang et al., 2011). The presence of repressive TF motifs at LMRs may prevent these putative poised enhancers from activating and also explains why some LMRs lack an enhancer chromatin signature. Other motifs enriched only in unmarked LMRs were those for EOMES/TBR2 and members of the ETS family. EOMES/TBR2 is well characterized for its role in controlling differentiation inradial glia toward neurons during early brain development (Arnold et al., 2008; Sessa et al., 2008). Eomes/Tbr2 also medi- ates Pax6 control of the balance of NSC self-renewal and neuro- genesis (Sansom et al., 2009). Members of the ETS family of TFs are known for their role in regulating a spectrum of develop- mental processes. For example, Etv1 is involved in differentiation of dopaminergic neurons (Flames and Hobert, 2009) and cere- bellum granule cells (Abe et al., 2011), while Ets1 plays a role in neural crest differentiation (Wang et al., 2015). Next, we assigned unmarked LMRs to NNGs (Table S5C). As validation, 85% of putative target genes identified by GREAT were also identified by NNG approach (Figure S5F). Putative tar- gets of unmarked LMRs were highly involved in development of the nervous system (Figure 5D; Figure S5G). Interestingly, the GO terms associated with these genes (e.g., transmission of the nerve impulse and axonogenesis) are indicative of a function at a later stage of development relative to NT formation, such as when neurons are formed, which would be consistent with a poised regulatory element. Furthermore, we contrasted unmarked LMRs with fetal brain and spinal cord DHSs and found that over half (23,448 of 39,078) are hypersensitive at later stages in development (Fig- ure 5E). To further supporting the poised for activation state of unmarked LMRs, we compared the expression of putative target genes of unmarked LMRs containing the REST motif between NRs and fetal spinal cord and found that the expression of these genes significantly increased at a later stage in development (Figure 5F). Upregulation of these genes suggest the elements switch to an active state. Collectively, these results support our hypothesis that un- marked LMRs may predict distal regulatory elements prior to acquiring chromatin features more commonly used to identify enhancers. Given that NRs and other 2D NSCs express a large repertoire of neural TFs that are present throughout devel- opment (see above and Ziller et al., 2015), these factors may establish poised enhancers (both H3K4me1-only and LMR- only elements), which activate upon localization of additional cell-specific factors. Folate-Associated CpGs Overlap NR Regulatory Elements High levels of folate during the periconceptional period are the most effective means to reduce the risk of NTDs (Zaganjor et al., 2016). Recently, a large study correlated a wide range of maternal folate levels during early stages of pregnancy with DNA methylation in cord blood at birth in 1,988 healthy mother-newborn pairs (Joubert et al., 2016). Methylation levels at 443 individual CpGs were correlated with folate levels. We sought to explore how our model of NT formation could unravel the regulatory mechanism whereby folate-associated methyl- ation influences NT formation. To do so, we assayed the genomic context for the 443 folate-associated CpGs and found that most were located in intergenic or intronic regions, where enhancers are more likely to be located (Figure S6A). Approxi- mately 40% of the folate-associated CpGs are at a regulatory element, with three-quarters of those overlapping H3K4me1- marked enhancers (Figure 6A). By contrast, overlapping folate- associated CpGs with 2D NSC H3K4me1-marked enhancersCell Reports 20, 1448–1462, August 8, 2017 1455 A B C D E F Figure 5. LMR Predictions of Distal Regulatory Domains and Their Developmental Dynamics (A) Pie chart showing the overlap of LMRs with chromatin-marked promoters and enhancers. (B) Heatmap representation of the ChIP-seq signal intensity (RPKM, input normalized) for H3K4me1, H3K4me3, H3K27me3, and H3K27ac over 2-kb windows centered on LMRs. (C) Heatmap representation of negative log10(p value) for significantly (q-value < 0.01) enriched TF binding site motifs at active enhancers, poised enhancers and unmarked LMRs. Only motifs whose TF was expressed (TPM > 1) in NR are shown. See Table S5B. (D) Bar plot of the most significantly enriched GO terms associated with NNGs of unmarked LMRs. See also Table S5C. (E) Heatmap representation of the DHS signal (RPKM) for human fetal brain and spinal cord (18 and 15 weeks of gestation, respectively) over 2-kb windows centered at unmarked LMRs. (F) Violin plot of expression (log2(TPM)) of NNGs in NRs and fetal spinal cord (15 weeks of gestation) for unmarked LMRs where the REST binding site motif was found to be significantly enriched. ***p < 0.0001 (t test). See also Figure S5 and Table S5.showed little to no overlap (Figure 6B), supporting a potential functional relevance of these CpGs in NT formation. Analysis of the folate-associate CpGs at NR regulatory elements showed that about half of the folate-associate CpGs are indeed de- methylated (methylation level <50%) in the NR methylome (Figure 6C).1456 Cell Reports 20, 1448–1462, August 8, 2017In linewith these results,GO termanalysis (Table S6A) of NNGs of the folate-associatedCpGsatNR regulatory elements showed enrichment for neurodevelopment andmetabolic processes (Fig- ure 6D; Tables S6B andS6C). 73%ofNNGswere validated using GREAT, andGO termanalysis ofGREATputative targets showed highconcordance (FiguresS6BandS6C).We found that over half -4.00 -2.00 0.00 regulation of metabolic process pattern specification process regulation of metabolic process anterior/posterior pattern specification ion transmembrane transport cell-cell signaling G-protein coupled receptor signaling pathway regulation of biosynthetic process generation of neurons regionalization cyclic nucleotide biosynthetic process positive regulation of gastrulation A B D HES3 APC2 H3K4me1 H3K27ac DNAme CpG 20 kb 20 kb E 20 7 88 32 6 9 12 LMR H3K4me1/LMR H3K4me1 H3K4me1/UMR UMR H3K4me3/UMR H3K4me3 C log10 p-value 0.00 0.25 0.50 0.75 1.00 M ethylation Level cell type % folate- associated CpGs Fisher's exact test p- value NR 29% 3.4E-93 ERG 1% 6.3E-01 NE 0% 1.0E+00 NPC 7% 8.6E-02 TBX3 10 kb Figure 6. Folate-Associated CpGs Overlap NR Regulatory Elements (A) Pie chart showing overlap of 174 folate-asso- ciated CpGs from Joubert et al. (2016) and NR regulatory elements. Fisher’s exact test log10 p < 5 for all overlaps shown. (B) Percentages of the 443 folate-associated CpGs identified by Joubert et al. (2016) over- lapping enhancers from NRs and 2D NSCs. Sta- tistical significance of the overlap was calculated using Fisher’s exact test. (C) Heatmap representation of methylation level at the 174 folate-associated CpGs that overlap NR regulatory elements. (D) Bar plot of the most significantly enriched GO terms associated with NNGs of regulatory ele- ments that overlap folate-associated CpGs. (E) UCSC Genome Browser snapshots of ChIP-seq signals (RPKM, input normalized) for H3K4me1, H3K4me3, H3K27me3, and H3K27ac at some loci that contain folate-associated CpGs. Peaks calls are shown as bars above each track. See also Figure S6 and Table S6.of the putative target genes for folate-associated CpGs at NR regulatory elements are expressed in NRs (TPM> 1; Figure S6D). Among putative target genes (Figure 6E), we found TBX3, a gene with expression restricted to NRs. TBX3 is associated with several GO terms related to NT patterning and axis establish- ment. Another related gene was FOXA2, a transcriptional acti- vator and effector of SHH signaling. Several genetic models showed that this pathway is important for proper NT closure and disruption of Shh signaling leads to NTDs (Murdoch and Copp, 2010). In mice, Foxa2 mediates Shh signaling in the floor plate, a critical organizing center for NT patterning (Cho et al., 2014; Mansour et al., 2011). Given the poised state of putative FOXA2 enhancers and that the gene is not expressed in NRs (TPM < 1), it may be possible that this gene is activated at a later stage of NT development. Notably, the APC2 locus harbors nine individual folate-associate CpGs, eight of which are within NR enhancers. APC2 is involved in the regulation of WNT signaling and is expressed in both the fetal and adult human brain. A study in mice showed an association between folate and DNA methylation in its homologous gene, APC (Sie et al., 2013). It is worth noting that several NNGs of enhancers that contain folate-associate CpGs are implicated in signaling pathways thatCell Repregulate NT axis establishment, that is dorsal-ventral and rostral-caudal axes. This process likely plays an important role in NT closure and NTDs. NT closure is a continuous process along the tube, and key signaling pathways, such as the SHH and Notch pathways, drive this process and play a role in regulating timing of closure along the NT (Greene and Copp, 2009). Interestingly, we observed enhancers that overlapped folate-associate CpGs are highly en- riched for broad enhancers (28 out of127; Fisher’s test p value = 3.8e-32), over three times more than expected, given that 6%–7% of all NR enhancers are broad. Among the putative target genes of broad enhancers that contain a folate-associate CpG is HES3, a regulator of the Notch signaling pathway, required to maintain NSCs population by inhibiting premature differentiation into neurons (Imayoshi et al., 2010) (Figure 6E). Mouse studies showed that Hes3 mutants develop NTDs (Hirata et al., 2001), likely as results of premature neural differentiation before NT closure is complete (Copp and Greene, 2013). These results reveal insights into mechanisms by which folate may influence neurodevelopment by altering methylation levels at distal regulatory elements such as enhancers. Furthermore, it demonstrates the functional relevance of our in vitro 3D model of human NT formation, which failed to be captured when analyzing regulatory maps in 2D NSCs (Figure 6B). NR Enhancers Are Enriched for Genetic Variants Associated with Neurological Disorders Recent studies have shown that non-coding single-nucleotide polymorphisms (SNPs) associated with risk for disease may be enriched at enhancer elements (Corradin and Scacheri, 2014;orts 20, 1448–1462, August 8, 2017 1457 Hawkins et al., 2010b, 2013). Given the limited genetic studies for NTDs that have resulted in very few candidates for this con- dition, a large collection of NTD-associated SNPs is not available for analysis at NR enhancers. Nonetheless, we asked whether we could observe enrichment for other neurological disorders or conditions potentially founded in neurodevelopment. To do so, we collected SNPs associated with a variety of neurological traits and other conditions from public data (573 traits from the NHGRI-EGI genome-wide association studies (GWAS) catalog; https://www.ebi.ac.uk/gwas/) and examined their overlap within NR enhancers (Tables S7A and S7B). We screened both lead SNPs and SNPs in linkage disequilibrium (r2R 0.8) for localiza- tion at NR enhancers. Our analysis showed that for 20% of these traits (123 out of 573), a significant enrichment (q-value < 0.05) of the associated SNPs was found at NR enhancers, of which 10 were neurological disorders and traits, including schizophrenia and autism (Figure 7A). For both these neurolog- ical disorders, an increasing body of evidence is emerging sup- porting a neurodevelopmental etiology (Donovan and Basson, 2017; Stachowiak et al., 2013). From the 10 enriched neurological disorders and traits, we found their overlapping SNPs to reside in 120NR enhancers (Fig- ure 7B), 59 of which were active (Figure 7C) and 101 were in a poised state (Figure 7D). Next, we sought to identify the putative target genes for the 120 NR enhancers. We leveraged a recent study that provides prediction of target genes for over 500,000 putative regulatory elements using correlations between DHS signal and gene expression levels across 72 human cell types (Sheffield et al., 2013). This allowed us to predict gene targets for 68 NR enhancers (Figure 7B). Among both active and poised enhancers, we observed that vast majority of putative target genes were expressed in NRs (Table S7B; Figures 7C and 7D). For example, among predicted gene targets expressed in NRs was autism susceptibility candidate 2 gene (AUTS2), a gene implicated in neurodevelopment and neurological disorders, including autism spectrum disorders, intellectual disability, and developmental delay (Oksenberg and Ahituv, 2013). Notably, SNPs associated with subcortical brain region volumes are more abundant at poised enhancers, and their target genes are largely off in NRs. Altered activation of these enhancers and target genes may provide insight into encephalopathies. Overall, our results suggested that our model can be useful for evaluating the impact of genetic risk factors for neurological dis- orders associatedwith dysfunction of early development. For var- iants at active enhancers, our results likely reveal insights into the developmental origins of these disorders. Given that just over half (63%) of the enhancer SNPs are within poised enhancers, the impact of these variants is likely realized later in development. DISCUSSION We generated an initial set of transcriptomic and epigenomic maps in a 3D human NSC culture system to build a comprehen- sive regulatory framework of hESC-derived NRs, which largely recapitulate NT development unlike 2D systems. A key feature of NRs is their unrestricted differentiative potential within the neural lineage (Elkabetz et al., 2008; Koch et al., 2009). Maintain- ing neural stemness and inhibiting premature progression to-1458 Cell Reports 20, 1448–1462, August 8, 2017ward neuronal fate is critical. Our results show that these key properties may be mediated by epigenomic features. Examina- tion of in silico TF-TF interaction maps derived from putative enhancer binders revealed two important nodes, which were occupied by SOX2 and TCF3, both of which are known to play a role in maintaining NSC identity (Atlasi et al., 2013; Graham et al., 2003; Kuwahara et al., 2014). Similarly, repression of TGF-beta signaling and activation of Notch signaling (Louvi and Artavanis-Tsakonas, 2006; Yoon and Gaiano, 2005) are two fundamental events in early neurodevelopment and estab- lishing NR identity. Our study suggests that a repressor of TGF-beta signaling, SKI, and an effector of Notch signaling, HES5, may be transcriptionally regulated via broad enhancers, which resemble super-enhancers/stretch enhancers typically enriched near master regulators. Another key feature of NRs is an anterior identity while still pos- sessing the capability of responding to patterning cues (Conti et al., 2005; Elkabetz et al., 2008; Pankratz et al., 2007). This developmental potential and ability to respond to signaling cues may be mediated by the largely poised enhancer land- scape. We found that broad enhancers, predominantly poised, were highly enriched in proximity to genes involved in NT region- alization. Furthermore, we identified motifs for TFs associated with the anterior and posterior phenotype at enhancers. The regulatory networks anchored in chromatin-marked regulatory elements provides a better understanding of the epigenetic mechanisms guiding gene regulation, whereby patterning of the NT is established and orchestrated. Little is known about the molecular mechanisms that link maternal folate levels to mechanisms of neurodevelopment and malformations such as NTDs during the early stages of hu- man embryogenesis (Nazki et al., 2014). Our study demonstrates the power of human in vitro modeling to fill this gap. This study highlights several major points. Our findings clearly demonstrate that NRs are a more relevant model for studying epigenetic mechanisms of folate, as we found little to no overlap between folate-associated CpGs and 2D NSC enhancer maps. This pro- vides a framework to explain how folate-associated effects on DNA methylation may in turn influence NT formation and NTD. Moreover, enhancers are likely a critical mediator of mecha- nisms that link folate to NT formation. Not only are the majority of reported folate-associated CpGs enriched outside of pro- moters and genic regions, but also, more specifically, our findings show that three-quarters of those overlapping NR regu- latory map are at enhancers, as defined by histone modifications or a DNA methylation signature. This is in line with the ever- increasing understanding of the central role that enhancers play in development and predisposition to disease. Future studies that make using of genome-wide DNA methylation maps are likely to reveal more extensive folate-associated CpGs, as studies conducted to date have utilized array technol- ogies with minimal CpG representation. Our study also provides insights into other neurodevelopment related disorders and traits. Enhancer variants offer an unprece- dented opportunity to unravel how genetic variants in the non- coding genome can influence risk for many diseases and traits. Cell-type-specific contexts are key for assessment of the effects of enhancer variants. Representing an early stage B D C –log10 p-value 52 36 20 4 1 0.1 0.01 A P er ce nt ag e O ve rla pp in g S N P s In te rs ec tio n S iz e (lo g1 0) 5 3 1 Set Size Enhancers Corrected Target Genes DHSs Enhancers Corrected Target Genes DHSs Figure 7. NR Enhancers Are Enriched for Genetic Variants Associated with Neurological Disorders (A) Bubble chart showing traits and conditions from public data (https://www.ebi.ac.uk/gwas/) for which both lead SNPs and SNPs in linkage disequilibrium (r2R 0.8) were found inside NR enhancers. Circle size represents the number of SNPs overlapped by enhancers for each trait. Colored circles represent traits associated with neurological disorders. Grey circles represent traits not associated with neurological disorders. Green dashed line indicates p value = 0.05 (hypergeometric test with Bonferroni correction). (legend continued on next page) Cell Reports 20, 1448–1462, August 8, 2017 1459 of neurodevelopment, NRs provide a useful tool for investigating neurological disorders and traits that may have a developmental origin, such as schizophrenia and autism. Furthermore, we were able to identify their putative target genes and found among them genes involved in autism spectrum disorders, intellectual disability, and developmental delay (Oksenberg and Ahituv, 2013). In conclusion, this study provides comprehensive epige- nomic maps of regulatory elements in NRs and supports the use of this NSC population as a 3D human model for studying NT formation and NTDs. EXPERIMENTAL PROCEDURES For detailed descriptions, see Supplemental Experimental Procedures. NR Differentiation hESCs (WA09, WiCell) were cultured on Matrigel in mTeRS1 (STEMCELL Technologies, Inc.). NRs were differentiated for 5 days in STEMdiff Neural Induction Medium according to the manufacturer’s guidelines and on Aggre- Well 800 plates (STEMCELL Technologies, Inc.), transferred to polyornithine- and laminin-coated plates, and cultured in STEMdiff Neural Induction Medium for 7 days. NR 3D tubular-structures were harvested at day 12 of differentiation. Genome-wide Assays, Library Construction, and Sequencing ChIP-seqwas performed as previously described (Hawkins et al., 2010a) using 2 million cells per ChIP. RNA-seq was performed using ScriptSeq v2 RNA-Seq Library Preparation Kit (Epicenter). WGBS was performed as previously described (Lister et al., 2009) with minor modifications using 500 ng genomic DNA per library. Libraries were sequenced either on NextSeq500 in single-end 75-cycle runs or on HiSeq2500 in paired-end 100-cycle runs. All experiments were performed in duplicate. Computational Analysis Sequencing raw reads were trimmed for low quality (q score < 30) and adapters,mapped to human genome (hg19) using Bowtie2 (ChIP-seq), Kallisto (RNA-seq), and BBMap (WGBS). Duplicates were merged. Gene expression values were calculated as TPM. Alternative pairwise differential expression analysis of RNA-seq was performed with RankProd R package. ChIP-seq peaks were called using MACS v1.4 using the nomodel mode. Methylation was quantified using BS-Seeker2; DMRs were identified using DSS and UMRs and LMRs using MethylSeekR. Publically Available Data Used in This Study DHS data for human fetal brain (18 weeks gestation) and fetal spinal cord (15 weeks gestation) were obtained from the Roadmap Epigenome Project (accession numbers GEO: GSM595913 and GEO: GSM878661). Folate-asso- ciated CpGs coordinates were obtained from Joubert et al. (2016). SNPs for all available traits were downloaded from the NHGRI-EGI GWAS catalog (https:// www.ebi.ac.uk/gwas/). RNA-seq, ChIP-seq, and WGBS raw data for NPCs were obtained from Xie et al. (2013); data for NE and ERG were obtained from Ziller et al. (2015). Statistical Analysis Statistical significance of differences in the distribution of expression values for active versus poised or bivalent genes as well as for promoters with hyper- versus hypo-DMRs was assessed using the t test. Statistical significance of differences between expression of NNGs for broad H3K4me1 or H3K27ac do- mains and all genes expressed was assessed using t test. Significance of(B) Bar chart shows number of enhancers overlapping neural related SNPs and D conditions are shown by linked dots below bar. (C and D) Alluvial plots showing the number of active (C) and poised (D) enhancers the count of target genes. Colored boxes show the number of target genes expr 1460 Cell Reports 20, 1448–1462, August 8, 2017the overlap of folate-associated CpGs at broad enhancers was determined using a Fisher’s test. Statistically significant overlap of trait-associated SNPs with NR enhancers was assessed using a hypergeometric test with Bonferroni correction.ACCESSION NUMBERS The accession number for the data generated on neural rosettes is SRA: SRP094721. All data generated on NRs as well as data for 2D NSCs used in this study can be visualized in a UCSC genome browser session at http:// depts.washington.edu/hawklab/ucsc_sessions/public/Neural_Rosette/site/nr_ page.html.SUPPLEMENTAL INFORMATION Supplemental Information includes Supplemental Experimental Procedures, six figures, and seven tables and can be found with this article online at http://dx.doi.org/10.1016/j.celrep.2017.07.036. AUTHOR CONTRIBUTIONS C.V. and R.D.H. conceived the study, designed the experiments and compu- tational analysis, and wrote the manuscript. C.V. performed the experiments and analyzed data. C.A. contributed to data analysis and writing of the manu- script. S.B. and R.L. contributed to DNA methylation data analysis. N.D. contributed to data analysis. R.J.L. provided assistance with cell culture. All authors read, commented on, and edited the manuscript. ACKNOWLEDGMENTS We thank Marjo Hakkarainen, Pa¨ivi Junni, and Bogata Fezazi for technical assistance with hESC culture. This study was supported by the National Insti- tutes of Health (NIH: R21ES021959), the Academy of Finland (259913), the Washington Life Sciences Discovery Fund (265508), and Biocenter Finland. S.B. is supported by a National Health and Medical Research Council of Australia and Australia Research Council (NHMRC-ARC) Dementia Research Development Fellowship grant (APP1111206). R.L. was supported by a Sylvia and Charles Viertel Senior Medical Research Fellowship and an ARC Future Fellowship (FT120100862). Received: December 10, 2016 Revised: May 31, 2017 Accepted: July 13, 2017 Published: August 8, 2017 REFERENCES Abe, H., Okazawa, M., and Nakanishi, S. (2011). The Etv1/Er81 transcription factor orchestrates activity-dependent gene regulation in the terminal matura- tion program of cerebellar granule cells. Proc. Natl. Acad. Sci. USA 108, 12497–12502. Abranches, E., Silva, M., Pradier, L., Schulz, H., Hummel, O., Henrique, D., and Bekman, E. (2009). Neural differentiation of embryonic stem cells in vitro: a road map to neurogenesis in the embryo. PLoS One 4, e6286. Arnold, S.J., Huang, G.-J., Cheung, A.F.P., Era, T., Nishikawa, S., Bikoff, E.K., Molna´r, Z., Robertson, E.J., and Groszer, M. (2008). The T-box transcription factor Eomes/Tbr2 regulates neurogenesis in the cortical subventricular zone. Genes Dev. 22, 2479–2484.HS (with associated target genes) as represented by bar height. Overlapping overlapped by neural condition SNPs, overlapping DHS with gene targets, and essed (on = TPM > 1) or off. Atlasi, Y., Noori, R., Gaspar, C., Franken, P., Sacchetti, A., Rafati, H., Mah- moudi, T., Decraene, C., Calin, G.A., Merrill, B.J., and Fodde, R. (2013). Wnt signaling regulates the lineage differentiation potential of mouse embryonic stem cells through Tcf3 down-regulation. PLoS Genet. 9, e1003424. Berk, M., Desai, S.Y., Heyman, H.C., and Colmenares, C. (1997). Mice lacking the ski proto-oncogene have defects in neurulation, craniofacial, patterning, and skeletal muscle development. Genes Dev. 11, 2029–2039. Boyles, A.L., Billups, A.V., Deak, K.L., Siegel, D.G., Mehltretter, L., Slifer, S.H., Bassuk, A.G., Kessler, J.A., Reed,M.C., Nijhout, H.F., et al.; NTDCollaborative Group (2006). Neural tube defects and folate pathway genes: family-based as- sociation tests of gene-gene and gene-environment interactions. Environ. Health Perspect. 114, 1547–1552. Burger, L., Gaidatzis, D., Schuebeler, D., and Stadler, M.B. (2013). Identifica- tion of active regulatory regions from DNA methylation data. Nucleic Acids Res. 41, e155. Cho, G., Lim, Y., Cho, I.-T., Simonet, J.C., and Golden, J.A. (2014). Arx together with FoxA2, regulates Shh floor plate expression. Dev. Biol. 393, 137–148. Conti, L., and Cattaneo, E. (2010). Neural stem cell systems: physiological players or in vitro entities? Nat. Rev. Neurosci. 11, 176–187. Conti, L., Pollard, S.M., Gorba, T., Reitano, E., Toselli, M., Biella, G., Sun, Y., Sanzone, S., Ying, Q.-L., Cattaneo, E., and Smith, A. (2005). Niche-indepen- dent symmetrical self-renewal of a mammalian tissue stem cell. PLoS Biol. 3, e283. Copp, A.J., and Greene, N.D.E. (2013). Neural tube defects–disorders of neurulation and related embryonic processes. Wiley Interdiscip. Rev. Dev. Biol. 2, 213–227. Corradin, O., and Scacheri, P.C. (2014). Enhancer variants: evaluating func- tions in common disease. Genome Med. 6, 85. Creyghton, M.P., Cheng, A.W., Welstead, G.G., Kooistra, T., Carey, B.W., Steine, E.J., Hanna, J., Lodato, M.A., Frampton, G.M., Sharp, P.A., et al. (2010). Histone H3K27ac separates active from poised enhancers and pre- dicts developmental state. Proc. Natl. Acad. Sci. USA 107, 21931–21936. Deglincerti, A., Etoc, F., Ozair, M.Z., and Brivanlou, A.H. (2016). Self-organiza- tion of spatial patterning in human embryonic stem cells. In Essays on Devel- opmental Biology Part A, P. Wassarman, ed. (Elsevier), pp. 99–113. Donovan, A.P.A., and Basson, M.A. (2017). The neuroanatomy of autism: a developmental perspective. J. Anat. 230, 4–15. Elkabetz, Y., Panagiotakos, G., Al Shamy, G., Socci, N.D., Tabar, V., and Studer, L. (2008). Human ES cell-derived neural rosettes reveal a functionally distinct early neural stem cell stage. Genes Dev. 22, 152–165. Esmailpour, T., and Huang, T. (2012). TBX3 promotes human embryonic stem cell proliferation and neuroepithelial differentiation in a differentiation stage- dependent manner. Stem Cells 30, 2152–2163. Flames, N., and Hobert, O. (2009). Gene regulatory logic of dopamine neuron differentiation. Nature 458, 885–889. Graham, V., Khudyakov, J., Ellis, P., and Pevny, L. (2003). SOX2 functions to maintain neural progenitor identity. Neuron 39, 749–765. Greene, N.D.E., and Copp, A.J. (2009). Development of the vertebrate central nervous system: formation of the neural tube. Prenat. Diagn. 29, 303–311. Hagey, D.W., Zaouter, C., Combeau, G., Lendahl, M.A., Andersson, O., Huss, M., and Muhr, J. (2016). Distinct transcription factor complexes act on a permissive chromatin landscape to establish regionalized gene expression in CNS stem cells. Genome Res. 26, 908–917. Hawkins, R.D., Hon, G.C., Lee, L.K., Ngo, Q., Lister, R., Pelizzola, M., Edsall, L.E., Kuan, S., Luu, Y., Klugman, S., et al. (2010a). Distinct epigenomic land- scapes of pluripotent and lineage-committed human cells. Cell Stem Cell 6, 479–491. Hawkins, R.D., Hon, G.C., and Ren, B. (2010b). Next-generation genomics: an integrative approach. Nat. Rev. Genet. 11, 476–486. Hawkins, R.D., Hon, G.C., Yang, C., Antosiewicz-Bourget, J.E., Lee, L.K., Ngo, Q.-M., Klugman, S., Ching, K.A., Edsall, L.E., Ye, Z., et al. (2011). Dynamicchromatin states in human ES cells reveal potential regulatory sequences and genes involved in pluripotency. Cell Res. 21, 1393–1409. Hawkins, R.D., Larjo, A., Tripathi, S.K., Wagner, U., Luu, Y., Lo¨nnberg, T., Ra- ghav, S.K., Lee, L.K., Lund, R., Ren, B., et al. (2013). Global chromatin state analysis reveals lineage-specific enhancers during the initiation of human T helper 1 and T helper 2 cell polarization. Immunity 38, 1271–1284. Heintzman, N.D., Hon, G.C., Hawkins, R.D., Kheradpour, P., Stark, A., Harp, L.F., Ye, Z., Lee, L.K., Stuart, R.K., Ching, C.W., et al. (2009). Histone modifi- cations at human enhancers reflect global cell-type-specific gene expression. Nature 459, 108–112. Hirata, H., Tomita, K., Bessho, Y., and Kageyama, R. (2001). Hes1 and Hes3 regulate maintenance of the isthmic organizer and development of the mid/ hindbrain. EMBO J. 20, 4454–4466. Hnisz, D., Abraham, B.J., Lee, T.I., Lau, A., Saint-Andre´, V., Sigova, A.A., Hoke, H.A., and Young, R.A. (2013). Super-enhancers in the control of cell identity and disease. Cell 155, 934–947. Holmberg, J., Clarke, D.L., and Frise´n, J. (2000). Regulation of repulsion versus adhesion by different splice forms of an Eph receptor. Nature 408, 203–206. Huang, Z., Wu, Q., Guryanova, O.A., Cheng, L., Shou, W., Rich, J.N., and Bao, S. (2011). Deubiquitylase HAUSP stabilizes REST and promotes maintenance of neural progenitor cells. Nat. Cell Biol. 13, 142–152. Imayoshi, I., Sakamoto, M., Yamaguchi, M., Mori, K., and Kageyama, R. (2010). Essential roles of Notch signaling in maintenance of neural stem cells in developing and adult brains. J. Neurosci. 30, 3489–3498. Joubert, B.R., den Dekker, H.T., Felix, J.F., Bohlin, J., Ligthart, S., Beckett, E., Tiemeier, H., van Meurs, J.B., Uitterlinden, A.G., Hofman, A., et al. (2016). Maternal plasma folate impacts differential DNA methylation in an epige- nome-wide meta-analysis of newborns. Nat. Commun. 7, 10577. Kelava, I., and Lancaster, M.A. (2016). Protocol review. Stem Cells 18, 736–748. Koch, P., Opitz, T., Steinbeck, J.A., Ladewig, J., and Br€ustle, O. (2009). A rosette-type, self-renewing human ES cell-derived neural stem cell with poten- tial for in vitro instruction and synaptic integration. Proc. Natl. Acad. Sci. USA 106, 3225–3230. Kuwahara, A., Sakai, H., Xu, Y., Itoh, Y., Hirabayashi, Y., and Gotoh, Y. (2014). Tcf3 represses Wnt-beta-catenin signaling and maintains neural stem cell population during neocortical development. PLoS One 9, e94408. Lee, K.E., Seo, J., Shin, J., Ji, E.H., Roh, J., Kim, J.Y., Sun, W., Muhr, J., Lee, S., and Kim, J. (2014). Positive feedback loop between Sox2 and Sox6 inhibits neuronal differentiation in the developing central nervous system. Proc. Natl. Acad. Sci. USA 111, 2794–2799. Lister, R., Pelizzola, M., Dowen, R.H., Hawkins, R.D., Hon, G., Tonti-Filippini, J., Nery, J.R., Lee, L., Ye, Z., Ngo, Q.-M., et al. (2009). Human DNA methyl- omes at base resolution show widespread epigenomic differences. Nature 462, 315–322. Louvi, A., and Artavanis-Tsakonas, S. (2006). Notch signalling in vertebrate neural development. Nat. Rev. Neurosci. 7, 93–102. Mansour, A.A., Nissim-Eliraz, E., Zisman, S., Golan-Lev, T., Schatz, O., Klar, A., and Ben-Arie, N. (2011). Foxa2 regulates the expression of Nato3 in the floor plate by a novel evolutionarily conserved promoter. Mol. Cell. Neurosci. 46, 187–199. McLean, C.Y., Bristor, D., Hiller, M., Clarke, S.L., Schaar, B.T., Lowe, C.B., Wenger, A.M., and Bejerano, G. (2010). GREAT improves functional interpre- tation of cis-regulatory regions. Nat. Biotechnol. 28, 495–501. Meinhardt, A., Eberle, D., Tazaki, A., Ranga, A., Niesche, M., Wilsch-Bra¨u- ninger, M., Stec, A., Schackert, G., Lutolf, M., and Tanaka, E.M. (2014). 3D reconstitution of the patterned neural tube from embryonic stem cells. Stem Cell Reports 3, 987–999. Merrill, B.J., Pasolli, H.A., Polak, L., Rendl, M., Garcı´a-Garcı´a, M.J., Anderson, K.V., and Fuchs, E. (2004). Tcf3: a transcriptional regulator of axis induction in the early embryo. Development 131, 263–274.Cell Reports 20, 1448–1462, August 8, 2017 1461 Morin, I., Platt, R., Weisberg, I., Sabbaghian, N., Wu, Q., Garrow, T.A., and Rozen, R. (2003). Common variant in betaine-homocysteinemethyltransferase (BHMT) and risk for spina bifida. Am. J. Med. Genet. A. 119A, 172–176. Murdoch, J.N., and Copp, A.J. (2010). The relationship between sonic Hedgehog signaling, cilia, and neural tube defects. Birth Defects Res. A Clin. Mol. Teratol. 88, 633–652. Nazki, F.H., Sameer, A.S., and Ganaie, B.A. (2014). Folate: metabolism, genes, polymorphisms and the associated diseases. Gene 533, 11–20. Niwa, H., Ogawa, K., Shimosato, D., and Adachi, K. (2009). A parallel circuit of LIF signalling pathways maintains pluripotency of mouse ES cells. Nature 460, 118–122. Oksenberg, N., and Ahituv, N. (2013). The role of AUTS2 in neurodevelopment and human evolution. Trends Genet. 29, 600–608. Pankratz,M.T., Li, X.-J., Lavaute, T.M., Lyons, E.A., Chen, X., and Zhang, S.-C. (2007). Directed neural differentiation of human embryonic stem cells via an obligated primitive anterior stage. Stem Cells 25, 1511–1520. Parker, S.C.J., Stitzel, M.L., Taylor, D.L., Orozco, J.M., Erdos, M.R., Akiyama, J.A., van Bueren, K.L., Chines, P.S., Narisu, N., Black, B.L., et al.; NISC Comparative Sequencing Program; National Institutes of Health Intramural Sequencing Center Comparative Sequencing Program Authors; NISC Comparative Sequencing Program Authors (2013). Chromatin stretch enhancer states drive cell-specific gene regulation and harbor human disease risk variants. Proc. Natl. Acad. Sci. USA 110, 17921–17926. Persson, M., Stamataki, D., te Welscher, P., Andersson, E., Bo¨se, J., R€uther, U., Ericson, J., and Briscoe, J. (2002). Dorsal-ventral patterning of the spinal cord requires Gli3 transcriptional repressor activity. Genes Dev. 16, 2865– 2878. Rada-Iglesias, A., Bajpai, R., Swigut, T., Brugmann, S.A., Flynn, R.A., and Wysocka, J. (2011). A unique chromatin signature uncovers early develop- mental enhancers in humans. Nature 470, 279–283. Romanoski, C.E., Glass, C.K., Stunnenberg, H.G., Wilson, L., and Almouzni, G. (2015). Epigenomics: Roadmap for regulation. Nature 518, 314–316. Sansom, S.N., Griffiths, D.S., Faedo, A., Kleinjan, D.-J., Ruan, Y., Smith, J., van Heyningen, V., Rubenstein, J.L., and Livesey, F.J. (2009). The level of the tran- scription factor Pax6 is essential for controlling the balance between neural stem cell self-renewal and neurogenesis. PLoS Genet. 5, e1000511. Sarkar, A., and Hochedlinger, K. (2013). The sox family of transcription factors: versatile regulators of stem and progenitor cell fate. Cell Stem Cell 12, 15–30. Schug, J., Schuller, W.-P., Kappen, C., Salbaum, J.M., Bucan,M., and Stoeck- ert, C.J., Jr. (2005). Promoter features related to tissue specificity as measured by Shannon entropy. Genome Biol. 6, R33. Sessa, A., Mao, C.-A., Hadjantonakis, A.-K., Klein, W.H., and Broccoli, V. (2008). Tbr2 directs conversion of radial glia into basal precursors and guides neuronal amplification by indirect neurogenesis in the developing neocortex. Neuron 60, 56–69. Sheffield, N.C., Thurman, R.E., Song, L., Safi, A., Stamatoyannopoulos, J.A., Lenhard, B., Crawford, G.E., and Furey, T.S. (2013). Patterns of regulatory ac- tivity across diverse human cell types predict tissue identity, transcription fac- tor binding, and long-range interactions. Genome Res. 23, 777–788. Shlyueva, D., Stampfel, G., and Stark, A. (2014). Transcriptional enhancers: from properties to genome-wide predictions. Nat. Rev. Genet. 15, 272–286. Sie, K.K.Y., Li, J., Ly, A., Sohn, K.-J., Croxford, R., and Kim, Y.-I. (2013). Effect of maternal and postweaning folic acid supplementation on global and gene- specific DNA methylation in the liver of the rat offspring. Mol. Nutr. Food Res. 57, 677–685.1462 Cell Reports 20, 1448–1462, August 8, 2017Smith, Z.D., and Meissner, A. (2013). DNA methylation: roles in mammalian development. Nat. Rev. Genet. 14, 204–220. Stachowiak, M.K., Kucinski, A., Curl, R., Syposs, C., Yang, Y., Narla, S., Ter- ranova, C., Prokop, D., Klejbor, I., Bencherif, M., et al. (2013). Schizophrenia: a neurodevelopmental disorder–integrative genomic hypothesis and thera- peutic implications from a transgenic mouse model. Schizophr. Res. 143, 367–376. Stadler, M.B., Murr, R., Burger, L., Ivanek, R., Lienert, F., Scho¨ler, A., van Nimwegen, E., Wirbelauer, C., Oakeley, E.J., Gaidatzis, D., et al. (2011). DNA-binding factors shape the mouse methylome at distal regulatory regions. Nature 480, 490–495. Stamataki, D., Ulloa, F., Tsoni, S.V., Mynett, A., and Briscoe, J. (2005). A gradient of Gli activity mediates graded Sonic Hedgehog signaling in the neural tube. Genes Dev. 19, 626–641. Swiss, V.A., and Casaccia, P. (2010). Cell-context specific role of the E2F/Rb pathway in development and disease. Glia 58, 377–390. Thurman, R.E., Rynes, E., Humbert, R., Vierstra, J., Maurano, M.T., Haugen, E., Sheffield, N.C., Stergachis, A.B., Wang, H., Vernot, B., et al. (2012). The accessible chromatin landscape of the human genome. Nature 489, 75–82. Voigt, P., Tee, W.-W., and Reinberg, D. (2013). A double take on bivalent pro- moters. Genes Dev. 27, 1318–1338. Wallingford, J.B., Niswander, L.A., Shaw, G.M., and Finnell, R.H. (2013). The continuing challenge of understanding, preventing, and treating neural tube defects. Science 339, 1222002. Wang, C., Kam, R.K.T., Shi, W., Xia, Y., Chen, X., Cao, Y., Sun, J., Du, Y., Lu, G., Chen, Z., et al. (2015). The proto-oncogene transcription factor Ets1 regu- lates neural crest development through histone deacetylase 1 to mediate output of bone morphogenetic protein signaling. J. Biol. Chem. 290, 21925– 21938. Whyte, W.A., Orlando, D.A., Hnisz, D., Abraham, B.J., Lin, C.Y., Kagey, M.H., Rahl, P.B., Lee, T.I., and Young, R.A. (2013). Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell 153, 307–319. Xie, W., Schultz, M.D., Lister, R., Hou, Z., Rajagopal, N., Ray, P., Whitaker, J.W., Tian, S., Hawkins, R.D., Leung, D., et al. (2013). Epigenomic analysis of multilineage differentiation of human embryonic stem cells. Cell 153, 1134–1148. Yoon, K., and Gaiano, N. (2005). Notch signaling in themammalian central ner- vous system: insights from mouse mutants. Nat. Neurosci. 8, 709–715. Zaganjor, I., Sekkarie, A., Tsang, B.L., Williams, J., Razzaghi, H., Mulinare, J., Sniezek, J.E., Cannon, M.J., and Rosenthal, J. (2016). Describing the preva- lence of neural tube defects worldwide: a systematic literature review. PLoS ONE 11, e0151586. Zhang, S.-C., Wernig, M., Duncan, I.D., Br€ustle, O., and Thomson, J.A. (2001). In vitro differentiation of transplantable neural precursors from human embry- onic stem cells. Nat. Biotechnol. 19, 1129–1133. Zhu, X., Li, B., Ai, Z., Xiang, Z., Zhang, K., Qiu, X., Chen, Y., Li, Y., Rizak, J.D., Niu, Y., et al. (2016). A robust single primate neuroepithelial cell clonal expan- sion system for neural tube development and disease studies. Stem Cell Reports 6, 228–242. Ziller, M.J., Edri, R., Yaffe, Y., Donaghey, J., Pop, R., Mallard, W., Issner, R., Gifford, C.A., Goren, A., Xing, J., et al. (2015). Dissecting neural differentiation regulatory networks through epigenetic footprinting. Nature 518, 355–359.