Transcriptome analysis reveals determinant stages controlling human embryonic stem cell commitment to neuronal cells

Proper neural commitment is essential for ensuring the appropriate development of the human brain and for preventing neurodevelopmental diseases such as autism spectrum disorders, schizophrenia, and intellectual disorders. However, the molecular mechanisms underlying the neural commitment in humans remain elusive. Here, we report the establishment of a neural differentiation system based on human embryonic stem cells (hESCs) and on comprehensive RNA sequencing analysis of transcriptome dynamics during early hESC differentiation. Using weighted gene co-expression network analysis, we reveal that the hESC neurodevelopmental trajectory has five stages: pluripotency (day 0); differentiation initiation (days 2, 4, and 6); neural commitment (days 8–10); neural progenitor cell proliferation (days 12, 14, and 16); and neuronal differentiation (days 18, 20, and 22). These stages were characterized by unique module genes, which may recapitulate the early human cortical development. Moreover, a comparison of our RNA-sequencing data with several other transcriptome profiling datasets from mice and humans indicated that Module 3 associated with the day 8–10 stage is a critical window of fate switch from the pluripotency to the neural lineage. Interestingly, at this stage, no key extrinsic signals were activated. In contrast, using CRISPR/Cas9–mediated gene knockouts, we also found that intrinsic hub transcription factors, including the schizophrenia-associated SIX3 gene and septo-optic dysplasia-related HESX1 gene, are required to program hESC neural determination. Our results improve the understanding of the mechanism of neural commitment in the human brain and may help elucidate the etiology of human mental disorders and advance therapies for managing these conditions.

Brain development is one of the most complicated and hierarchical events in mammals. A series of temporal processes, including neural commitment, patterning and sub-regionalization, neurogenesis, and neuronal network formation, are required to generate a functional brain (1)(2)(3)(4). Defects in any of these processes will lead to neurodevelopmental diseases such as autism spectrum disorders, schizophrenia, depression, and epilepsy. Neural commitment occurs after gastrulation and initiates the program to form neural epithelium (5)(6)(7). Many studies on amphibians and rodents reveal that the inhibition of BMP 3 signaling is an evolutionally conserved mechanism for neural induction; in addition, intrinsic factors play pivotal roles to direct fate transition from pluripotency to neural lineage (8 -11). Nonetheless, how the neural commitment is ensured in human brain has not been clearly elucidated. Systematically understanding the critical role of signaling pathways and transcriptional factors in modulating and shaping human brain development is hindered by limited accessibility to early embryos and inadequate amounts of stage-specific and cell type-specific materials. These problems may now be solved by the use of human embryonic stem cells (hESCs). Establishment of in vitro differentiation models that recapitulate normal development will facilitate the study in brain development and neurological disorders.
The establishment of neural differentiation protocols for hESCs makes it possible to investigate early events, including neural commitment in humans (12)(13)(14)(15). hESCs exhibit the restricted capacity to generate various subtypes of functional neurons by responding to extrinsic signals (16 -19), which reca-pitulate brain development in vivo. Global gene expression profiling technologies such as microarray and RNA sequencing (RNA-Seq) enable highly sensitive analysis of transcriptome associated with the neurodevelopment of hESCs. By using RNA-Seq data with nine samples from days 0 to 77 of the differentiating hESCs, van de Leemput et al. (20) establish a COR-TECON system to study human cerebral cortex development in vitro. Meanwhile, single-cell RNA-Seq is applied to identify various classes of neural progenitors (NPCs) and neurons, which is helpful for mapping the early human brain cells (21). Microarray analysis shows that the non-canonical WNT signaling pathway is important for partitioning neural versus epidermal fate during neural induction (22). It has been shown that the early neurodevelopment of hESCs advances much quicker than that in vivo (13,15,23). Therefore, the insufficient representation of differentiating time points analyzed by RNA-Seq or the low resolution of the microarray technique limits the outcome of systematic analysis on fast and transient cell fate transition such as neural induction.
In this study, we adapted and developed an in vitro hESC neural differentiation system, ending up with a high percentage of dorsal forebrain neurons. By specific co-expression gene assays of transcriptome data with 12 samples prepared every other day between differentiation day 0 and day 22, we show that the following five distinct stages exist during the early neural differentiation of hESCs: pluripotency (day 0); differentiation initiation (day 2/4/6); neural commitment (day 8/10); NPC proliferation (day 12/14/16); and neuronal differentiation stage (day 18/20/22). Expression profiling comparison of gene modules and transcription factor (TF) gene groups among several systems reveals that the Module 3-associated day 8/10 stage is a critical window for the fate transition from the pluripotency to the neural epithelium. Moreover, PAX6, SIX3, SIX6, HESX1, and ID3 are identified as key hub TF genes of this stage. The loss-of-function of either the SIX3 or HESX1 gene, mediated by the CRISPR/Cas9 gene editing system, leads to compromised neural commitment of hESCs.

Directed differentiation of hESCs mimics the early cortical development in vivo
To investigate the regulatory mechanisms of human neural commitment, we first adapted the previous protocols (12) and standardized an in vitro hESC (H9 line) neural differentiation system, with EB formation for 6 days, attached EB (aEB) for 10 days, sphere in N2 for 6 days, and then single cells replated in N2B27 for 4 weeks (Fig. 1A). Twelve RNA samples were collected every other day from day 0 to day 22 and subjected to RNA sequencing. Quantitative real-time PCR (RT-qPCR) assays with 12 RNA samples were performed to briefly characterize the differentiation process of hESCs. Following the differentiation, the expression of pluripotent genes POU5F1, SOX2, NANOG was decreased, and the expression of neuroectoderm genes POU3F1 and ZNF521, as well as neural epithelial marker genes such as PAX6 and SOX1, was increased and reached the peak at day 12. The expression of anterior forebrain progenitor marker genes FOXG1, OTX2, and SIX3 was up-reg-ulated at around day 16, followed by the elevation of neuronal marker genes TUJ1 (TUBB3) and MAP2 around days 16 -22 (Fig. 1B). The similar expression pattern could be obtained with another hESC line H1 (supplemental Fig. S1A), and the correlation coefficient between H9 and H1 cells is high (R 2 ϭ 0.95, supplemental Fig. S1B). To further confirm the synchronization of the neural differentiation protocol, we randomly collected 36 single cells at days 0, 4, 8, 14, and 20, respectively, and performed single cell qPCR for OCT6 (POU3F1), SIX3, and SOX1 genes. The results show that the majority of single cells show the comparable expression level for each gene, and the expression pattern of these genes is similar to the results from population cell samples (supplemental Fig. S1D). Consistent with RNA-expressing patterns, immunostaining assays revealed that many cells in day 10 aEB were double-positive for the NPC markers PAX6 and NESTIN (62.5%, Fig. 1C and supplemental Fig. S1C); in addition, PAX6 and OTX2 were co-expressed in those NPCs, indicating a cerebral cortical identity (71.3%, Fig.  1C supplemental Fig. S1C). Numerous NPCs in day 14 aEBs were PAX6 and SOX1 double-positive (85.1%, Fig. 1C supplemental Fig. S1C). At day 30, the majority of the cells not only displayed the prolonged neurites but also was positive for mature neuron markers MAP2 and NEUN (82.8%, Fig. 1C supplemental Fig. S1C). At differentiation day 50, the predominant neuronal subtype was TBR1ϩ (78.8%) and VGLUT1/2ϩ (88.6%) positive cortical glutamatergic projection neurons (Fig. 1, D and E). There were few THϩ dopaminergic neurons (5.1%) but no VAChTϩ and GAD67ϩ neurons (Fig. 1, D and E).
The results above demonstrate that our protocol can successfully program hESCs to differentiate into cortical NPCs and then cortical projection neurons, and this procedure is reliable and could be applied to different hESC lines.

Transcriptome analysis reveals five possible sub-stages of the hESC neural differentiation
Next, RNA-Seq assays with 12 RNA samples were conducted. Approximately 30 million sequencing reads of every sample were mapped to the human genome hg19. The average number of detectable genes was about 15,000 in each sample with fragments per kilobase per million (FPKM) more than 0.1 in at least one sample (supplemental Fig. S2, A and B). Correlation analysis confirmed that the normalized RNA-Seq tag counts of marker genes were consistent with their expression levels assessed by RT-qPCR ( Fig. 2A). Unsupervised hierarchical clustering analysis revealed that there were three main stages during the early hESC neural differentiation, designed as EB (days 0, 2, 4, and ), aEB (days 8, 10, 12, 14, and 16), and Sphere (days 18, 20, and 22) stage, respectively (Fig. 2B). Within the EB stage, days 2, 4, and 6 were closely clustered, compared with day 0. In the aEB stage, there were two sub-clusters, one with days 8 and 10 and the other with days 12, 14, and 16 (supplemental Fig.  S2C). Principal component analysis (PCA) revealed that the three major stages could be further divided into five sub-groups ( Fig. 2C and supplemental Fig. S2D). PC1 accounted for the temporal trajectory of differentiation, and PC2 separated the differentiation stages in more detail. As distributed on the transition line from the PC1-negative direction to the PC1-positive direction, day 8/10 probably is a critical period Determinant stages controlling hESC neural differentiation of fate transition. Consistently, clustering analysis using the highest PC-loading genes in the PC1 (24) (top 100 PC1-negative and -positive genes) suggested that day 8 in cluster 3 is a critical fate transition period during hESC neural differentiation (supplemental Fig. S2E). Especially, functional enrichment analysis showed that the genes with positive PC loadings in PC1 showed enrichment for forebrain development, neuron development, and axonogenesis. In contrast, the genes with negative PC loadings in PC1 were associated with RNA metabolic process, apoptosis, and cell cycle (Fig. 2D). Taken together, the dynamic gene expression profiles revealed a temporal developmental trajectory of hESC neural differentiation, and day 8/10 of cluster 3 is a critical fate transition period.

Gene-network modules define five temporal sub-stages of hESC neural differentiation
To systematically investigate the correlation of gene expression profiles, we performed weighted gene co-expression network analysis (WGCNA) to identify distinct co-expression modules among the clusters of associated transcripts without supervision and bias (25)(26)(27). We first focused on modules that can distinguish five sub-stages during the hESC neural differentiation. Among dozens of modules, five modules, marked in cyan, dark orange, light yellow, brown, and steel blue, were tightly associated with these five sub-stages, respectively (Fig.  3A). The expression of genes identified in each module was specifically enriched in the corresponding sub-stage (Fig. 3B).

Determinant stages controlling hESC neural differentiation
The expression of genes in Module 1 was down-regulated, whereas the expression of genes in Module 5 was up-regulated gradually along the differentiation of hESCs (Fig. 3, A-C).
To characterize five modules individually, functional enrichment analysis was conducted for genes in each module. Module 1 (day 0) was enriched for pluripotency related terms, e.g. regulation of cell motion, migration, and biological adhesion. Module 2 (day 2/4/6) was associated with stem cell differentiation and embryonic morphogenesis. The terms of chromatin assembly, forebrain development, and neural tube patterning were identified in Module 3 (day 8/10). Cell proliferation and cell growth were enriched in Module 4 (day 12/14/16). Module 5 (day 18/20/22) was associated with neuron differentiation, axonogenesis, and neuron development (Fig. 3B). To chart the activity of signaling pathways in the different modules, we examined the enrichment of target genes of important signaling pathways by using RNA-seq or microarray data from published perturbation experiments (29), and we found many signaling pathways were involved during hESC neural differentiation (Fig. 3D); some pathways are activated (red) and some are inhibited (green) in different modules, but there was barely any detectable signaling activity in Module 3, which is associated with day 8/10 ( Fig. 3D). The enrichment of neural related genes and the lack of extrinsic signaling genes in Module 3 support that day 8/10 is a critical period for the fate transition from pluripotency to neural lineage.
Based on the signature gene profiles in the defined modules, the early neural differentiation of hESC can be classified into five stages, including pluripotency stage with Module 1 at day 0, differentiation initiation stage with Module 2 at day 2/4/6, neu-

Determinant stages controlling hESC neural differentiation
Determinant stages controlling hESC neural differentiation ral commitment stage with Module 3 at day 8/10, NPC proliferation stage with Module 4 at day12/14/16, and neuronal differentiation stage with Module 5 at day18/20/22 (Fig. 3E). Taken together, we reveal that the hESC neurodevelopmental trajectory contains five stages, characterized by unique module genes, and day 8/10 is associated with neural commitment.

Comparative analysis validates Module 3 associated day 8/10 stage as a critical period for the neural commitment
In mouse early embryo development, the neural fate commitment occurs during gastrulation, and a portion of anterior ectoderm is specified to adopt the neural fate (28,29). One of our previous studies documented the spatial transcriptome of the epiblast at the mid-gastrulation (E7.0) (29). Pearson correlation coefficient (PCC) assays with RNA-Seq data of all the embryo regions and five hESC differentiation stages showed that Module 3 was specifically enriched with the prospective ectoderm in the anterior embryo at E7.0, where the neural epithelium originates later (Fig. 4A). In contrast, Module 2 was widely associated with the whole epiblast, indicating the status of the pluripotency (Fig. 4A). In mouse ESC neural differentiation, the neural fate commitment begins when NPC-specific gene Sox1 starts to express (30). By comparing with the microarray data of in vitro mouse ESC neural differentiation, we noticed that day 8 of Module 3 was highly correlated with the day 3 Sox1 ϩ cells, indicating the initiation of mouse neural fate (Fig. 4B). Meanwhile, Module 2 was associated with cells in pluripotent status (Fig. 4B). All the findings above demonstrated that neural fate commitment occurs at day 8/10 stage during hESC differentiation.
The CORTECON system is generated by analyzing RNA-Seq data with nine samples from day 0 to day 77 of the differentiating hESCs (20). PCC assays revealed that the CORTECON system could not reveal the conversion from the pluripotent epiblast to the neural lineage, which progresses quickly during the early development of the human nervous system (Fig. 4C). Transcriptome profiling data, generated by microarray from various laser-microdissected human fetal brain regions at different stages (31), are available in the BrainSpan Atlas of the Developing Human Brain (http://brainspan.org/). 4 The correlation analysis demonstrated that our hESC neural differentiation system resembles the development of early human brain at gestational week 8/9 (Fig. 4D).
Next, we asked whether there were novel markers that could be used to distinguish each stage. We performed Guilt-by-Association analysis (24,32) to identify putative sub-stagespecific markers (Fig. 4E). As shown in supplemental Fig. S3, the expression of putative markers in each stage was activated sequentially following the differentiation of hESCs. Then, RT-qPCR assays were conducted to confirm the expression profiles of selected genes in each module (supplemental Fig. S4A). GAL, GNA14, CALB1, and TUBA4A were novel marker genes in Module 1 for pluripotency; SALL3, PCDH1, and IRX2 in Module 2 for differentiation initiation; FGF8 and PTN in Module 3 for neural fate commitment; ENPP2, SEMA5A, and NEDD9 in Module 4 for NPC proliferation; and STMN2 and EMX2 in Module 5 for neuronal differentiation (supplemental Fig. S4, A  and B). The expression patterns of those novel marker genes with known stage-specific genes, e.g. OCT4, OCT6, PAX6, SOX1, and NCAM, during the hESC neural differentiation are summarized in supplemental Fig. S4B.
Overall, systematic comparison between our RNA-Seq data and other transcriptome profiling data suggests that day 8/10 stage is a critical window for the fate transition from the pluripotency to the neural lineage.

Regulation network of transcription factor genes during hESC neural differentiation
Given that intrinsic regulators, especially TFs, play essential roles in the neural commitment (20,33,34), we sought to identify key TF genes in each of the five modules. By overlapping the stagespecific genes with TF database (35) (www.transcriptionfactor. org), 4 distinct TF genes were enriched in individual modules (Fig. 5A). The gene ontology analysis showed that the stage transition could possibly be driven by these TFs during the hESC neural differentiation. For instance, TFs in Module 1 were enriched for the regulation of transcription, gene expression, cellular biosynthetic process, and cell proliferation (supplemental Fig. S5A). TFs in Module 2 were related to embryonic morphogenesis, tube development, tissue morphogenesis, and ectoderm development (supplemental Fig. S5B). TFs in Module 3 were associated with forebrain development, cell fate commitment, neural tube development, and midbrain-hindbrain development (supplemental Fig. S5C). TFs in Modules 4 and 5 were involved in cell proliferation and neuronal differentiation (supplemental Fig. S5, D and E).
To unravel the inter-group connection of TFs in each module, we analyzed the Connection Specificity Index (CSI) (36) of the five TF groups, and generated a TF co-expression network (CSI Ͼ0.9). The positive correlations between the different TFs are marked in red, and the negative correlations are marked in green (Fig. 5B). TFs in each module showed positive interactions with each other (Fig. 5B), indicating that TFs within the module may form a synergistic regulatory circuitry and act in a combinatory manner. TFs in Modules 1 and 2 mainly showed positive correlation, whereas TFs between Modules 1 and 5 displayed negative correlation. Interestingly, TFs in Module 3 that presented the critical period of the neural fate transition displayed unique correlations with other TF groups (Fig. 5B). They were mainly negatively correlated with TFs in Modules 1 and 2, inhibiting earlier processes; they were positively corre-

Determinant stages controlling hESC neural differentiation
lated with TFs in Modules 4 and 5, promoting later events, suggesting that they might behave like a gatekeeper to mediate the lineage transition from the pluripotency to the neural fate during the hESC differentiation (Fig. 5B).
TFs with the highest degree of connectivity in all modules are termed hub genes (31), which are expected to play imperative roles during the differentiating process. Hub genes of each module were identified and listed in Fig. 5C. PAX6, SIX3, SIX6,  HESX1, and ID3 were the top hub TFs with the highest CSI in Module 3 (Fig. 5C). In mice, both Pax6 and Six3 genes are essential for the appropriate forebrain patterning (37)(38)(39)(40). As expected, PAX6 and SIX3 centering in the hub generated the most connections with other modules (Fig. 5C). Both PAX6 and SIX3 were negatively correlated with hub TF genes in Module 1 and/or 2 but positively associated with hub TF genes in Modules 4 and 5 (Fig. 5C). SIX3 and SIX6 displayed similar regulation patterns with PAX6 (Fig. 5C). However, HESX1 and ID3 were negatively correlated with the hub genes in Module 5, and HESX1 positively correlated with the hub genes in Module 2. Thus all the data above support a notion that Module 3, which is associated with the neural fate commitment, inhibits Modules 1, 2, and 5, but it promotes Module 4 (Fig. 5D). The hub TF analysis also indicated that the cooperation of key TFs not only governs the progress in each sub-stage but also mediate the transition between temporal stage transition, which generates a dynamic and hierarchical TF regulatory network during hESC neural differentiation.

SIX3 and HESX1 play imperative roles in hESC neural differentiation
Among five hub genes of Module 3, the functions of PAX6 in human neurodevelopment have been well studied (41,42). To validate the function of other hub TF genes in Module 3, the CRISPR/Cas9 gene edition system was used to target SIX3, SIX6, HESX1, and ID3 in hESCs. Unfortunately, we failed to generate a stable hESC colony with either SIX6 or ID3 gene deletion (data not shown), indicating that SIX6 and ID3 genes might be essential for the maintenance of hESCs. The sgRNAs with complementary sequence to the coding sequence of the SIX3 or HESX1 gene were used to achieve gene edition (supplemental Table S2). Two individual colonies were screened for the editing of each gene. The sequencing data revealed that SIX3-sgRNA #1 (SIX3 #1 for short) and SIX3-sgRNA #2 (SIX3 #2 for short) had a 1-bp deletion at different loci of the SIX3 genome, and HESX1-sg RNA #1 (HESX1 no. 1 for short) and HESX1-sgRNA #2 (HESX1 #2 for short) had the same 12-bp

Determinant stages controlling hESC neural differentiation
deletion (data not shown). The Western blot assays confirmed that compared with the mock control at differentiation day 12, the expression of the SIX3 protein was barely detected in both colonies, and the expression of the HESX1 protein was significantly decreased in the two edited colonies (Fig. 6A). The two gene-edited colonies of either SIX3 or HESX1 gene could be passaged normally; in addition, the expression of OCT4, SOX2, and NANOG genes was comparable between those colonies and the control at both RNA and protein levels (supplemental Fig. S6, A and B), suggesting that these genes do not affect the pluripotency maintenance.
RT-qPCR assays revealed that at differentiation day 12, the expression of pluripotent genes OCT4 and NANOG was increased, but the neural epithelium marker genes OTX2, PAX6, SOX1, and ZIC1 was significantly reduced (Fig. 6B).
Consistently, Western blot assay confirmed that the expression of the PAX6 protein in the gene-edited colonies of either gene was also reduced (Fig. 6C).
Next, we asked what were the consequences of loss-of-function of the SIX3 or HESX1 gene during hESC neural differentiation. Immunocytochemical assays were performed to detect the expression of OCT4, SOX2, and OTX2 at day 6 EBs, and of PAX6 at day 12 attached EBs, and of NEUN at day 50 neurons in the control and SIX3 KO or HESX1 KO colony. In day 6 EBs, compared with the control, there were more OCT4 ϩ /SOX2 ϩ cells in either the SIX3 KO or HESX1 KO colony; but much fewer OTX2 ϩ /SOX2 ϩ cells (Fig. 6, D and E), suggesting that the early neural differentiation process was inhibit by loss-of-function of SIX3 or HESX1. In day 12 attached EBs, compared with the control, there were much fewer PAX6 ϩ cells in either SIX3

Determinant stages controlling hESC neural differentiation
KO or HESX1 KO colony (Fig. 6, F and G). In day 50 neurons, there were many NEUN ϩ cells in the control; however, in either the SIX3 KO or HESX1 KO colony, very few NEUN ϩ cells were detected (Fig. 6, H and I), suggesting that neuronal differentiation was blocked by SIX3 or HESX1 KO. Similar results were obtained in another clone of SIX3 and HESX1 deletion cells (supplemental Fig. S7). Taken together, these data suggest that the neural differentiation of hESCs is compromised in the lossof-function of either the SIX3 or HESX1 gene.

SIX3 and HESX1 promote neural differentiation by regulating downstream TF networks
According to the correlation analysis of hub genes among different modules (Fig. 5C), the potential downstream target genes of either SIX3 or HESX1 were identified (Fig. 7, A and D). If SIX3 or HESX1 promotes hESC neural differentiation through its downstream target genes, the loss-of-function of SIX3 or HESX1 gene should cause the up-regulation of its negatively correlated genes but the down-regulation of its positively correlated genes. Indeed, the expression of SIX3 negatively correlated genes NANOG, SMAD7, and LHX4 in Module 1 was significantly increased in both SIX3 KO colonies at day 0 (Fig. 7B); in contrast, the expression of positively correlated genes such as NR2F2, ZEB2, and LHX2 was reduced (Fig. 7C). In both HESX1 KO colonies, similar expressing changes were observed for the downstream target genes of HESX1, including HES3, ZIC5, ZIC2, and NR6A1 in Module 2 (Fig. 7E), as well as FOS, BACH2, CREBS, and INSM1 in Module 5 (Fig. 7F). These results suggest that SIX3 and HESX1 genes may intrinsically promote neural differentiation by regulating its downstream TF networks.

Discussion
In this study, we establish an in vitro hESC neural differentiation system, which mimics the early human brain develop-

Determinant stages controlling hESC neural differentiation
ment. Twelve RNA samples were collected from differentiating hESCs every other day between day 0 and day 22 for RNA-Seq. Transcriptome analysis reveals that five distinct co-expression gene modules can classify the hESCs neural differentiation into five sub-stages. The gene expression profiling and CSI analysis of transcription factor correlation suggest that Module 3 associated day 8/10 stage is a critical period for the neural fate commitment. PAX6, SIX3, SIX6, HESX1, and ID3 are identified as key hub TF genes of Module 3. CRISPR/Cas9 -mediated functional assays reveal that the neural differentiation of hESCs is impaired with the loss-of-function of either SIX3 or HESX1 gene.
Many studies with mouse systems have contributed significantly to the understanding of the development and the function of mammalian brain (43)(44)(45). Nevertheless, there are notable differences between the mouse and human brains, such as a 1,000-fold difference in size, more complicated human cortical layers, and advanced emotion-and stress-related functions of the human brain. Especially, mouse models of neural diseases cannot fully recapitulate the characteristics of human disorders (46). Neural differentiation of hESCs has been demonstrated to be a powerful system to study human brain development (47)(48)(49)(50)(51). In this study, we adapted and established an in vitro hESC neural differentiation system (Fig. 1). Similar to previous reports (14,15), RT-qPCR and immunofluorescent assays demonstrate that our hESC neural differentiation system is reliable and reproducible to resemble the early human cortical development in vivo ( Fig. 1 and supplemental Fig. S1). The early neural development of hESCs progresses very quick (13,15,23). The insufficient time points analyzed by RNA-Seq or the low resolution of the microarray technique in some previous studies limits the precision analysis of early neurodevelopmental events (20,22). Our transcriptome analysis, based on RNA-Seq data of 12 samples, reveals that during the early neural differentiation of hESCs, there are five temporal stages as follows: pluripotency (day 0); differentiation initiation (day 2/4/6); neural commitment (day 8/10); NPC proliferation (day 12/14/16); and neuronal differentiation stage (day 18/20/22) (Figs. 2 and 3). PCC assays with our RNA-Seq data and hESC CORTECON show that day 8/10 stage may mediate the fate switch from the pluripotent epiblast to the neural epithelium (Fig. 4). All the evidence above suggests that our hESC neural differentiation system is suitable to investigate the molecular and cellular mechanisms of early human brain development such as neural commitment.
In the past century, neural commitment has been well studied in amphibian and rodents (8,10). The inhibition of BMP signaling seems revolutionarily conserved during the neural fate initiation in different organisms (9,13,52). Inhibition of the TGF␤ pathway has now been demonstrated to be sufficient to directly induce neural fate in mammalian embryos as well as pluripotent mouse and human embryonic stem cells (53), and consistently, our data also showed that the TGF␤ signaling has been inhibited during hESC neural differentiation processing (Fig. 3D). However, the involvement of different extrinsic signaling pathways and intrinsic factors, including transcriptional and epigenetic factors, suggests that the regulatory machinery of neural determination is not that simple. So far, how the neu-ral fate is initiated in humans has not been clearly elucidated. In our study, the unique expression profiles of genes in PCA and WGCNA analysis suggest that among the five temporal substage, the Module 3 associated day 8/10 stage is a critical period for the fate transition (Figs. 2 and 3 and supplemental Fig. S2). Functional enrichment analysis reveals that Module 3 is associated with forebrain development (Fig. 3); in addition, TFs in Module 3 are related to cell fate commitment and neural tube development ( Fig. 5 and supplemental Fig. S5C). Moreover, systematic comparison between our RNA-Seq data and other transcriptome profiling data such as mouse Geo-Seq (29), mESC microarray (22), hESC CORTECON (20), and human fetal brain microarray (31) suggests that day 8/10 stage is essential for the neural determination (Fig. 4). Furthermore, no active extrinsic signals are present in Module 3 of day 8/10 ( Fig.  3), which supports the notion that the neural commitment is ensured by the intrinsic program. Finally, TFs in Module 3, as well as hub TFs, including PAX6, SIX3, SIX6, HESX1, and ID3, are identified ( Fig. 5 and supplemental Table S1). The CSI assays show that TFs in Module 3 are negatively correlated with TFs in Modules 1 and 2 but are positively associated with TFs in Modules 4 and 5 (Fig. 5). As expected, the loss-of-function of either SIX3 or HESX1 gene leads to the compromised neural differentiation of hESCs ( Fig. 6 and supplemental Fig. S7), which will be discussed in detail later. Thus, all our findings above demonstrate that day 8/10 stage is a critical window for the fate transition from the pluripotency to the neural lineage.
The intrinsic regulation is imperative for the neural induction. Nonetheless, how the key intrinsic factors, especially transcriptional factors, participate in the regulation of neural commitment is still largely unclear. The determination of day 8/10 stage as a neural fate conversion period provides an opportunity to explore the mechanism of human neural induction. PAX6, SIX3, SIX6, HESX1, and ID3 are hub TF genes in Module 3 related to day 8/10 ( Fig. 5). It has been known that PAX6 is not only a marker of human neural stem cells but is also a human neuroectoderm cell fate determinant (41,54,55). In humans, various mutations of the SIX3 gene, mapped to 2q21, are associated with holoprosencephaly (56 -60) and schizophrenia (61). Different mutations of the HESX1 gene, mapped to 3p14.3, are related to septooptic dysplasia (62)(63)(64)(65) and pituitary deficiency (63, 66 -68). In addition, the expression of Six3 and Hesx1 genes is highly enriched in the prospective E7.0 mouse anterior ectoderm (29), where the prospective neural epithelium originates, suggesting that SIX3 and HESX1 may play a crucial role in human neural commitment as well. Indeed, the early neural differentiation of hESCs is repressed with either SIX3 or HESX1 deficiency, achieved by the CRISPR/Cas9 gene edition system ( Fig. 6 and supplemental Fig. S7). Several possible downstream target genes of either SIX3 or HESX1 are identified by the correlation analysis. As expected, in the SIX3-deficient cells, the expression of the negatively associated genes such as NANOG, SMAD7, and LHX4 is increased; in contrast, the expression of the positive correlated genes, including NR2F2, ZEB2, and LHX2, is decreased (Fig. 7, A-C). Similar results are obtained for the HESX1-regulated genes in the HESX1-deficient cells (Fig. 7, D-F). Thus, our findings provide evidence that SIX3 and HESX1 are novel key intrinsic factors,

Determinant stages controlling hESC neural differentiation
which ensure human neural commitment by regulating the distinct transcriptional network. Our system can be used not only to identify the pivotal determinants in temporal events during human brain development, but also to unravel their mechanisms to program the appropriate neurodevelopment and to prevent the neural diseases.
Our in vitro hESC neural differentiation system recapitulates early human cortical development in vivo; especially the day 8/10 stage resembles the neural commitment in human brain. Many neural diseases such as autism spectrum disorders, schizophrenia, depression, epilepsy, and mental retardation are caused by the abnormal early neurodevelopment. Thus, our findings will benefit not only the understanding of the molecular mechanism of human brain development, but also the understanding of the etiology and the therapy of human mental diseases.

Human ESC culture
Human ESC lines H9 and H1 (passage 30 -45) were maintained and passaged every 6 -7 days on a feeder layer of irradiated embryonic mouse fibroblasts as described previously (12,69). The standard medium was supplemented with 4 ng/ml fibroblast growth factor (R&D Systems Inc.), and the differentiated colonies were physically removed before passaging.

Human ES cell neural differentiation
The procedure for hES cell differentiation is schematically summarized in Fig. 1A. Briefly, human ES cell colonies were detached from the feeder layer and floated in the ES cell growth medium without basic FGF. EB were formed after 4 days' suspension in the medium, and subsequently floated for 2 days in neural induction medium consisting of F-12/DMEM, N2 supplement, and non-essential amino acids (NEAA). The ES cell aggregates (EBs) were then adhered to a substrate in a neural induction medium. By about 10 days after ES cell differentiation, cells in the center of each colony differentiated into neuroectodermal cells, displaying small columnar morphology followed by organization of the columnar cells into neural tubelike rosettes after an additional 6 days. These neural epithelial cells were isolated from surrounding non-neural cells through differential response to dispase treatment. These aggregates of neural epithelial cells are floating and cultured in the same neural induction medium for another 6 days to form the neural sphere. Finally, these the neural progenitor spheres were digested with accutase (Innovation Cell Technology) into single cells and plated on PDL-coated dishes in the N2B27 medium for maturation (12).

Immunocytochemistry and cellular quantification
Immunocytochemistry staining on dish cultures was performed as described previously (15). The primary antibodies used in this study are listed in supplemental  Table S3). Images were taken with Nikon Eclipse FN1 confocal laser-scanning microscope. To quantify the differentiation efficiency, three to five fields were randomly selected, and the neuronal cells and total cells (DAPIstained) were counted using ImageJ software. All the experiments were performed in triplicate. A two-tailed Student's t test was used for statistical analysis.

Quantitative reverse-transcription PCR (qRT-PCR)
Total RNA was extracted using TRIzol reagent (Invitrogen), and cDNA was reverse-transcribed using the SuperScript III First-strand cDNA synthesis kit (Invitrogen). Quantitative RT-PCR was performed using the MyiQ real-time PCR detection system (Bio-Rad), as described (14). A minimum of three biological replicates from two separate experiments were examined.

Single-cell qPCR
The single-cell qPCR follows the protocol as published previously (70). Briefly, the single cell is collected into lysis buffer for the first cDNA synthesis, and a modified Smart2-seq protocol is performed to amplify the cDNA. cDNAs that pass quality control are then used for qPCR. The procedure has been optimized for each step: tissue collection; cell lysis; RNA isolation; and single cell-based PCR amplification. Quantitative RT-PCR was performed using the MyiQ real-time PCR detection system (Bio-Rad), as described (14). A minimum of three biological replicates from two separate experiments were examined.

CRISPR/Cas9 gene edition
The single-cell qPCR follows the protocol as published previously (71). First, pX330 plasmid, including the T7 promoter, was linearized by NotI. Linearized templates were purified and transcribed in vitro with mMESSAGE mMACHINE T7 ULTRA kit (Life Technologies, Inc.). sgRNAs with the T7 promoter were amplified by PCR and transcribed in vitro using MEGAshortscript T7 kit (Life Technologies, Inc.). After transcription, the Cas9 sgRNAs were purified with MEGAclear kit (Life Technologies, Inc.) according to the manufacturer's instructions. The sequences for preparation of template for in vitro transcription of sgRNA are shown in supplemental Table S2.

Determinant stages controlling hESC neural differentiation
Western blotting was carried out using horseradish peroxidaseconjugated IgG as a secondary antibody and enhanced chemiluminescence system for detection (ABclonal).

RNA-seq library construction
The RNA-seq library construction follows the method published previously (73). cDNA samples were sheared by ultrasonication on a Covaris S2 for 80 s with the following parameter settings: duty cycle ϭ 10; cycles per burst ϭ 200; intensity ϭ 4. Samples were ethanol-precipitated and resuspended in Tris buffer for subsequent enzymatic modifications. End repair was carried out with T4 DNA polymerase and T4 polynucleotide kinase (New England Biolabs), and samples were purified with Ampure XP paramagnetic beads (Beckman Coulter). Bluntend, 3Ј-phosphorylated products were 3Ј-adenylated with exo-Klenow fragment in the presence of dATPs (New England Biolabs), purified with Ampure XP beads, and ligated to sequencing adapters (Illumina) by T4 DNA ligase at 20°C for 30 min. PCR amplification of library constructs was carried out with AccuPrime DNA polymerase (Invitrogen) for 13 cycles. Molarity and size distribution of sequencing libraries were assessed by HS-DNA microfluidic chips on the 2100 Bioanalyzer (Agilent). Sequencing was performed in 100-bp pairedend format on the Illumina HiSeq 2000.

RNA-Seq data preprocessing
Raw reads were mapped to the hg19 genome using TopHat2 version 2.0.4 program (74). We calculated FPKM as expression level using Cufflinks version 2.0.2 with default parameters (75). Then, we discarded genes that do not have FPKM Ͼ0.1 in at least one sample within all 12 time points during the hESC neural differentiation, and we next transformed expression levels to log-space by taking the log 2 (FPKM ϩ 1).

Weighted gene co-expression network analysis
A signed and weighted correlation network was constructed by first creating a matrix of pairwise correlations between all pairs of genes across the measured samples (26). Next, the adjacency matrix was constructed. Soft power parameter was estimated and is interpreted as a soft threshold of the correlation matrix. Based on the resulting adjacency matrix, we calculated the topological overlap, which is a robust and biologically meaningful measure of network interconnectedness (i.e. the strength of two genes' co-expression relationship with respect to all other genes in the network). Genes with highly similar co-expression relationships were grouped together by performing average linkage hierarchical clustering on the topological overlap. We used the Dynamic Hybrid Tree Cut algorithm to cut the hierarchical clustering tree, and we defined modules as branches from the tree cutting. We summarized the expression profile of each module by representing it as the first principal component (referred to as module eigengene). Modules whose eigengenes were highly correlated (correlation above 0.75) were merged.

Functional enrichment analysis
Functional enrichment of gene sets with different expression patterns was performed using the Database for Annotation, Visualization, and Integrated Discovery version 6.7 (76).