CRISPR-Cas9 Screens Reveal Genes Regulating a G0-like State in Human Neural Progenitors

The coordination of developmental potential and proliferation in stem and progenitor cells is essential for mammalian development and tissue homeostasis. We performed CRISPR-Cas9 screens in human neural progenitor cells (hNPCs) and identified genes, including CREBBP, NF2, PTPN14, TAOK1, or TP53, that limit expansion. Knockout of these genes causes increased hNPC proliferation via skipping of a transient G0-like state, characterized by expression of genes associated with quiescent neural stem cells and neural development and molecular features of quiescent cells (e.g., hypophosphorylated Rb, low CDK2 activity, and p27 stabilization). Single-cell RNA-sequencing of hNPCs revealed distinct G0/G1 populations, altered in G0-skip mutants through both distinct and convergent downstream effectors, including cell cycle, Hippo-YAP, and novel targets. Our results provide a molecular and phenotypic portrait of expanding hNPCs including a gene expression map of their cell cycle and characterization of antiproliferative factors that regulate cell cycle exit with likely roles in maintaining developmental potential.


Introduction
In mammals, most developing and adult tissues are hierarchically organized such that tissue growth and maintenance is driven by the production of lineage-committed cells from populations of tissue-resident stem and progenitor cells (Reya, Morrison et al., 2001). In adult tissues, stem cells are typically found in a quiescent or reversible G0 state, but must re-enter the cell cycle and divide to promote lineage commitment (Doetsch, 2003, Obernier, Cebrian-Silla et al., 2018. Their lineage committed progeny, e.g., amplifying progenitors, further balance lineage potential with proliferation to produce enough terminally differentiated cells to keep pace with demand (Lin, 2008). This balance of proliferation and self-preservation is particularly crucial for neural progenitor cells which cannot compensate for differences in cell number by altering cell size due to the functional requirements of the resulting neural networks (Homem, Repic et al., 2015). Thus disruption of the balance between self-renewal and proliferation during development may lead to severe defects in brain architecture, such as microcephaly or megacephaly (Homem et al., 2015).
Neural progenitor cells (NPCs) derived from the developing mammalian telencephalon (Davis & Temple, 1994, Johe, Hazel et al., 1996 can be cultured ex vivo yet recapitulate the expansion, specification, and maturation of each of the major cell types in the mammalian central nervous system. When grown in defined, serum-free conditions in vitro in the presence of epidermal and fibroblast growth factors (EGF and FGF), these cells retain progenitor-like molecular features, including Nestin and Sox2 expression, and capacity for multi-lineage differentiation (Pollard, Conti et al., 2006, Sun, Pollard et al., 2008. As such, these cells represent a model for both understanding basic stem/progenitor cell biology and a tractable system for functional genomics to uncover the regulatory pathways that maintain the precise regulation of cell cycle exit.
Here, we identify genes that govern the self-renewal and proliferation of hNPCs in vitro by CRISPR-Cas9 genome-wide screening for genes rate-limiting for expansion.
We found that knock-out (KO) of five genes, CREBBP, NF2, PTPN14, TAOK1, or TP53, all known or candidate tumor suppressors, robustly increased hNPC expansion as a result of shortening G1 phase and, more specifically, skipping a transient G0-like state of variable length.
In-depth characterization and analysis of hNPCs using single-cell RNAsequencing revealed four definable G0/G1 populations, which are altered in G0-skip mutants. In particular, we observe that the G0-skip mutants promote expansion of a late G1 population enriched in key cell cycle regulators (e.g., MYC, CCND1) and Hippo-YAP targets at the expense of a G0-like subpopulation enriched for expression of neurodevelopment genes and those found in quiescent NSCs (e.g., CLU, HOPX).
Further, characterization of G0/G1-specific gene expression for each G0-skip gene KO yielded a rich data resource for analysis of the p53 transcriptional network, Hippo-YAP targets, cell cycle gene regulation, and many novel targets and pathways, including those downstream of CREBBP and TAOK1. The results suggest that hNPC expansion is constrained by a transient G0-like state, regulated by multiple pathways, that facilitates retention of neurodevelopmental identity.
These screens revealed dozens of candidate antiproliferative screen hits, which were enriched for known tumor suppressor genes (Futreal, Coin et al., 2004) (Fig S1B) and genes found mutated across 35 different cancer types ( Fig 1D). In contrast to NPCs, there was little evidence of expansion enhancement in CRISPR-Cas9 screens for patient-derived glioblastoma stem cell-like (GSCs) isolates (Joo, Kim et al., 2013) grown in the same conditions (Fig S1C), consistent with pre-existing, cancer-associated disruptions of key proliferative and tumor suppression pathways.
Examining the intersection of all of our screen data revealed five reproducible and robust proliferation-enhancing screen hits: CREBBP, NF2, PTPN14, TAOK1, and TP53 ( Fig 1E), which we validate, examine, and discuss below.

Validation of CREBBP, NF2, PTPN14, TAOK1, and TP53 as antiproliferative genes in hNPCs
To control for off-target effects, we performed small guide RNA (sgRNA) tiling screens for each of these genes, whereby each gene is targeted with 138 to 466 sgRNAs tiled across most exons. At least 70% of tiling sgRNAs for each candidate proliferation limiting gene were significantly enriched in NPCs, whereas an sgRNA tiled gene essential for DNA replication, MCM2, showed significant depletion over time and control, non-targeting (NTC) sgRNAs were largely inert (Figs 2A and S2A; Table S1).
KO of target genes was confirmed by Western blots (Fig 2B) and by deep sequencing analysis for CREBBP, NF2, and TP53 (Toledo et al., 2015). KO of CREBBP,NF2,PTPN14,TAOK1,and TP53 in hNPCs caused a significant proliferative advantage over control cells in a 23-day outgrowth competition assay, while KO of the essential gene KIF11 showed the opposite result ( Fig 2C).
However, the competitive advantage did not appear to be based on differences in survival since no changes in Annexin-V staining were observed following normal culturing or in co-cultures, where apoptosis remained <2% regardless of the experimental condition.
Using cell proliferation assays (Fig 2D-G), we found that each KO significantly increased cell accumulation in 48-96 hour outgrowth assays. Importantly, this effect was independent of cell density, as KO cells showed increased proliferation at both low and high densities (e.g., Fig 2E,F). Further, the doubling time significantly decreased for each KO, shortening from ~50 hours to 30-40 hours (Fig 2H), similar to two GSC isolates used in the same assay.
We also measured transit time through G0/G1 and S/G2/M in individual NPCs using time-lapse microscopy ( Fig S2B-E). For G0/G1 transit times, we found that our control hNPCs exhibit variable G1 transit times, consistent with observations of a transient G0-state seen in other non-transformed human and rodent cells (Spencer, Cappell et al., 2013). We observed a wide distribution of G0/G1 transit times in control hNPCs, from fast (4.3 hrs), medium, and extremely slow (95 hrs) (averaging 32.5 hrs) (Fig 3d). By contrast, S/G2/M transit times were much more uniform (~12.4 hrs).
In mammals, cell cycle ingress is governed by progressive phosphorylation of Rb by CDK4/6 and CDK2 as cells pass through the restriction point in late G1, causing derepression of E2F transcription factors (Sherr & McCormick, 2002, Weinberg, 1995, Yao, Lee et al., 2008, Zetterberg, Larsson et al., 1995. We observed that KO of CREBBP, NF2, PTPN14, TAOK1, or TP53 in U5-NPCs results in a pronounced increase in the intensity of phosphorylated Rb during G1, consistent with an enrichment for a late G1 state (Fig 4B).
CDK2 activity correlates with cell cycle progression; if CDK2 activity levels are low during G1, cells enter G0 (Spencer et al., 2013). If CDK2 activity is intermediate (relative to its peak during G2/M) they progress past the restriction point and into Sphase (Spencer et al., 2013). Using the steady-state cytoplasmic to nuclear ratios of a DNA helicase B (DHB)-mVenus reporter as a readout of CDK2 activity (Hahn, Jones et al., 2009, Spencer et al., 2013, we observed significant increases in CDK2 activity in each KO in G0/G1 cells (Fig 4C,D). This was true either by total intensity or the proportion of cells with a reporter ratio greater than 1, which corresponds with S-phase entry in mammary epithelium (Spencer et al., 2013). Control cells averaged ~8% of G1 cells with >1 cytoplasmic:nuclear reporter ratios CDK2 activity, while KOs were 20-27% ( Fig 4D).
Collectively, the above data demonstrate that KO of antiproliferative genes in U5-NPCs causes a cell autonomous decrease in cell cycle length with less distributed and faster G0/G1 transit times, an increase in the molecular features associated with late G1, and a reduction in the molecular features associated with G0 ( Fig 4G). These data are consistent with KOs either blocking entry of cells into a transient G0 state or causing failure to maintain cells in G0. Therefore, we call these G0-skip genes.
G0-skip genes have overlapping but distinct roles in regulating G0/G1 gene expression To further characterize G0-skip genes, we performed gene expression analysis of KO cells specifically in G0/G1 phase. Among our G0-skip genes, it has been previously reported that TP53/p53 and its transcriptional target CDKN1A/p21 have causal roles in modulating entry into G0 (Arora, Moser et al., 2017, Barr, Cooper et al., 2017, Spencer et al., 2013, Yang, Chung et al., 2017. Consistent with this notion, CDKN1A KO scored as a proliferation enhancer of U5-NPCs using the 3rd generation CRISPR-Cas9 library ( Fig S1A). Thus, one key question is whether each of the other G0-skip genes merely acts through p53 and/or p21 to regulate G0. CREBBP/CBP, for example, has been implicated in positively regulating p53 turnover and activity via acetylation (Ito, Lai et al., 2001).
Other G0-skip genes have known or suspected roles in negatively regulating Hippo-YAP pathway signaling which regulates cell and tissue growth and response to contact inhibition (Yu, Zhao et al., 2015). For the latter, repression of Hippo signaling is sufficient to arrest cells in G0/G1 at high densities (Zhao, Wei et al., 2007); however, our phenotypes occur at sub-confluent densities, suggesting that precocious Hippo-YAP activation rather than bypass of contact inhibition would be driving Hippo-YAP-related effects. NF2 (Neurofibromin 2, Merlin), PTPN14, and TAOK1 have been demonstrated or proposed to negatively regulate YAP, the major downstream effector of HIPPO signaling, through modulating LATS1/2 kinase activity (Plouffe, Meng et al., 2016, Wilson, Li et al., 2014, Zhang, Bai et al., 2010 or, in the case of PTPN14, through preventing nuclear localization of YAP (Lin, Poon et al., 2013).
To analyze gene expression, RNA-seq was performed on mCherry-CDT1+ sorted NPCs after KO, which captures both G0 and G1 subpopulations (Fig 5A and   Table S2). In control NPCs, as expected, comparing G0/G1 sorted cells to unsorted populations revealed down-regulation of genes involved in cell cycle regulation, DNA replication, and mitosis (Fig 5A and Table S3). Overall comparisons between the KOs and NTC U5-NPCs showed that KO of NF2 and PTPN14 were most similar by unsupervised clustering as well as having the most overall gene changes, while TAOK1 KO was most similar to the controls ( Fig 5B). However, comparison of the overlapping up-or down-regulated genes showed that TAOK1 KO up-regulated genes were more similar to NF2 and PTPN14 KO than the other KOs ( Fig S3A).
First evaluating p53 target genes, we found that only TP53 KO significantly down-regulated the expression of high confidence p53 targets including BAX, CDKN1A/p21, RRM2B, and ZMAT3 (Fischer, 2017) (Figs 5C and S3B). None of the other KOs showed inhibition of p53 targets or p53 itself, strongly suggesting that the G0-skip genes are not acting through p53-dependent transcriptional activity.
For evaluation of the transcriptional targets of the HIPPO-YAP pathway, we examined the expression of 55 YAP conserved gene targets (Cordenonsi, Zanconato et al., 2011). Each KO, except for CREBBP, showed significant enrichment for YAP targets with NF2 KO having increased expression of the largest subset (Figs 5C and S3C). Interestingly, NF2 KO activated one subset of YAP targets important in the biological process of extracellular matrix (ECM) organization, while TAOK1 KO activated a different subset of YAP targets important in nuclear chromosome segregation, such as during mitosis (Fig S3C-E). NF2 and PTPN14 KO shared the most overlap in YAP target activation, including targets considered universal Hippo-YAP targets (e.g., CTGF, CYR61, and SERPINE1).
We next performed gene ontology (GO) analysis for genes up-or down-regulated after KO (The Gene Ontology, 2017) (Figs 5D,E and S4A,B; Supplementary Table 4).
For up-regulated genes, the only common biological processes related to extracellular matrix and structure organization. Among up-regulated gene sets, NF2 and PTPN14 KO again showed the most overlap with 127 GO categories, while up-regulation of various Hippo-YAP pathway members, including LATS2, TEAD1, and YAP1, suggests a possible feedback regulation of the pathway unique to NF2 and PTPN14 KO (Fig S4C).
There were also a number of notable categories unique to some KOs among upregulated genes. Intriguingly, TAOK1 KO triggered up-regulation of >40 key regulators of mitosis (e.g., AURKA, BUB1, CCNB1/2, CDK1, KIF11,etc.), suggesting it may act to inhibit their precocious activation in G0/G1 or expression after mitosis has completed ( Fig S4D). CREBBP KO caused up-regulation of key nuclear-encoded mitochondrial genes, including members of the NADH dehydrogenase complex, the succinate dehydrogenase complex, and mitochondrial DNA polymerase ( Fig S4E). Since many of these genes are direct transcriptional regulatory targets of nuclear respiratory factors 1 and 2 (NRF1 and NRF2) (Kelly & Scarpulla, 2004), CREBBP likely acts to oppose NRF1/2 function in G0/G1 and to attenuate cellular respiration and mitochondrial function.
Analysis of down-regulated gene sets shared among all KOs revealed significant enrichment for genes involved in neural development with categories such as nervous system development, regulation of neurogenesis, and glial cell differentiation (Figs 5E and S4B). In addition, many other components of neural differentiation are downregulated among subsets of the G0-skip genes, including regulation of synaptic plasticity, membrane potential, oligodendrocyte differentiation, and dendrite development, suggesting that G0-skip genes help maintain neurodevelopmental gene expression in G0/G1.
The S, S/G2, and G2/M clusters were enriched for genes whose expression peaks in these phases (Cyclebase (Santos, Wernersson et al., 2015)) including CCNE2 in S-phase and CCNB1/2 in G2/M (Figs 6C and S5C) and produced tightly interconnected networks of key cell cycle genes ( Fig S5D). This enrichment included the DNA replication genes PCNA, MCM3/4/5/6/7/10, GMNN, and RPA2/3 for S-phase and the mitosis genes CDC20, AURKA, and BUB1 during G2/M. There were four definable G0/G1 clusters (Figs 6B and S5B). Despite being the largest cluster, the G1 cluster had the smallest number of enriched identifying genes which included IGFR1 signaling genes (e.g., IGFBP3 and IGFBP5), and significant reductions of genes expressed in S, S/G2, and G2/M clusters ( Fig S5C). The M/early G1 cluster showed low but significant residual expression of M phase genes and enrichment for splicing factor genes. The Late G1 cluster was defined by genes important in G1 cell cycle progression, including CCND1 and MYC, and enriched for cholesterol biosynthesis, cell adhesion genes, and the subset of YAP target genes activated following NF2/PTPN14 KO, such as CTGF and SERPINE1 (Figs 6C and S6A).
Finally, the Neural G0 cluster also showed significant repression of genes peaking in other phases of cell cycle, including suppression of CCND1 expression, which is an indicator of cell cycle exit (Sherr, 1995). Moreover, this population uniquely showed significant enrichment for many types of neural development gene categories, including glial cell differentiation, neurogenesis, neuron differentiation, and oligodendrocyte differentiation (Table S5). These genes included transcription factors with known roles in balancing stem cell identity and differentiation, including OLIG2, HEY1, HOPX, BEX1, SOX2, SOX4, and SOX9 (Bergsland, Werme et al., 2006, Sakamoto, Hirata et al., 2003, Scott, Wynn et al., 2010. Network analysis of mean and overlapping cluster gene expression resolved the seven clusters into a pattern that fits well with cell cycle progression and predicted transit through G0/G1 (Fig 6D). The gene expression of the candidate G0 population only significantly overlapped with G1 cluster, consistent with a cell cycle exit from G1.
Interestingly, for our S, S/G2, G2/M, and M/early G1 cluster genes, we found significant overlap with the activated NSCs of Llorens-Bobadilla et al. (2015) and the NPCs of Artegiani et al. (2017) (Fig S6D), which are no longer quiescent.
Lastly, we applied our cluster criteria to a scRNA-seq dataset from mouse embryonic stem cells (mESCs) sorted by live Hoechst staining into G1, S-phase, and G2/M populations (Buettner, Natarajan et al., 2015). We found that our G1 category captured 83% of their Hoechst G1 cells, our G2/M category captured 89% of their G2/M, and their S-phase cells were split between G1, S, and G2/M, which is consistent with their Hoechst S-phase gate overlapping portions of these populations ( Fig S6C).
However, the mESCs failed to score in our Neural G0, Late G1, or M/early G1 categories. This is consistent with the shorter G1 of ESCs compared to somatic cells (Coronado, Godet et al., 2013) and suggests that these G0/G1 subpopulations may be unique to adult or tissue-specific stem cells.
These results demonstrate that hNPCs have definable G0 and G1 subpopulations and that the candidate G0 population is uniquely enriched for expression of neurodevelopment genes, including those associated with quiescent NSCs.

KO of G0-skip genes results in transcriptional reprogramming of genes associated with G0/G1 NPC subpopulations
We next examined if genes defining each of the scRNA-defined clusters were significantly changed in the RNA-seq profiles of KO NPCs (Fig 7A,B). The analysis revealed that genes associated with the Neural G0 cluster were significantly downregulated after KO of CREBBP, NF2, PTPN14, TAOK1, or TP53. These genes included those expressed in quiescent NSCs and others cited above with key roles in neural development (e.g., CLU, HOPX, ID3, PTN, PTPRZ1, SOX2, and SOX4) ( Fig S7A).
By contrast, we found that genes associated with the late G1 cluster, including CCND1 and MYC, are significantly up-regulated in each KO ( Fig S7B). Further, we find KO TAOK1, and to a lesser degree, PTPN14, caused up-regulation of S/G2, G2/M, and M/early G1 associated genes.
To further confirm these results, scRNA-seq was performed on G0/G1-sorted hNPCs with KO of TAOK1, since its KO produced the largest degree of cross-cluster G0/G1 gene expression changes ( Fig 7A). The Neural G0 population and, to a lesser degree, G1 in TAOK1 KO cells is significantly diminished from 21.3% to 10.3% and 58.9% to 53.3%, respectively (Fig 7C-E). However, the late G1 population is increased (from 3.0% to 9.8%) as are cells in the M/early G1 (from 7.8% to 15.3%) and G2/M phase (from 1.5% to 4.4%). The expansion of the M/early G1 in TAOK1 KO cells explains and cross-validates the increase in mitotic genes observed in the bulk G0/G1 RNA-seq data in TAOK1 KO cells (Fig S4D), suggesting that TAOK1 helps attenuate expression of mitotic genes from the previous cell cycle.
These results strongly suggest that G0-skip mutants lose a significant fraction of the hNPC Neural G0 subpopulation and reprogram G1 transcription networks to promote entry into G1-S. We also examined the G0/G1 gene expression profiles from two patient GSC isolates to see if we could identify changes in similar G0/G1 cell cycle subpopulations in cancer cells (Fig S7C,D and Table S2). We found that genes associated with G1 or Neural G0 were depleted in the sorted GSCs, consistent with the decreased G0/G1 transit time seen through time-lapse FUCCI analysis ( Fig 3D). There was a greater increase in genes related to S and S/G2 in G0/G1 than in the NPCs and limited changes in genes associated with Late G1, suggesting that there may be other G1 subpopulations altered in the GSCs.

Loss of G0-skip genes attenuates developmental programs during hNPC differentiation
Since G0-skip mutants affect subpopulations of NPCs enriched for expression of neural genes, we further examined if KO of G0-skip genes affects developmental transcriptional programs triggered upon hNPC differentiation. To this end, we differentiated hNPCs by removal of mitogenic growth factors, FGF2 and EGF. These conditions trigger specification of multiple neural cell types over the course of 8 weeks, including TUJ1+ neurons, GFAP+ astrocytes, and a small number of O4+ oligodendrocytes. After only 7 days of mitogen-withdrawal, U5-NPCs show increased markers of the astrocytic and neuronal lineage (GFAP and TUJ1, respectively) and radial glia markers (GFAP, VIM, GLI3, PAX6, SOX2), reduction of neural stemnessrelated genes (e.g., NES/Nestin, PROM1/CD133), and an elongated, bipolar morphology ( Fig S7E and Table S6). There are also increases in 93 brainspecific/restricted transcription factors (e.g., BCL11B/CTIP2, NEUROD1, POU3F2/BRN2, RFX4, ZBTB18, ZIC4, etc.) (Gray, Fu et al., 2004) (Table S6) and genes associated with cell-cell adhesion, extracellular matrix organization, synapse assembly, and axoneme assembly. Consistent with exit from the cell cycle during this differentiation, we find genes involved in mitosis, DNA replication, and G1/S transition are significantly reduced (Table S6).
Among G0-skip genes, NF2 and PTPN14 KO more closely resembled each other and undifferentiated cells ( Supplementary Fig 7f). Further, their KO also prevented upregulation of most of the 93 neural transcription factors (59 total) up-regulated during early differentiation, while other KOs overlapped these 59 to varying degrees ( Fig   S7G). These included GATA2, which drives NPCs out of the cell cycle concomitant with down-regulation of cyclin D1 expression and up-regulation of p27/Kip1 (El Wakil, Francius et al., 2006), and FOXJ1, which controls NPC fate determination (Li, Floriddia et al., 2018). Moreover, we found that G0-skip mutants disrupt the early neural differentiation transcriptional programs, including reduced GFAP expression, reduced cilia-related genes (including DRC1), and reduced synapse-assembly genes (Figs 7F and S7H).
CREBBP and PTPN14 KO cells were more mitogen-independent, maintaining enrichment for cell cycle genes even after the withdrawal of FGF2 and EGF to induce differentiation (Table S7 and Fig 7G). In contrast, TP53 KO cells showed the most mitogen-dependence, with almost all cell cycle genes examined down-regulated similar to and even beyond values in control cells. NF2 KO cells displayed a partial effect with Late G1/S genes up-regulated (e.g., CCND1, CCNE1, and MYC), but no observable enrichment for key mitotic genes.
Thus, KO of G0-skip genes affect neurodevelopmental gene expression during conditions favoring lineage commitment, suggesting that the Neural G0 subpopulation may play an important role in maintaining fidelity of neurodevelopmental transcription programs.

DISCUSSION
Here we report the identification and characterization of five genes, CREBBP, NF2, PTPN14, TAOK1, and TP53, that function in hNPCs to limit proliferation through the promotion of a definable G0-like state ( Fig 7H). KO of these genes causes However, a second restriction point was recently proposed for cell cycle exit, which is governed by the amount of residual CDK2 activity from the previous cell cycle (Spencer et al., 2013). In mammalian cell lines, this restriction point is modulated by the presence of endogenous DNA damage and mitogen-driven Ras signaling during the previous cell cycle, which determine entry into a transient G0-like state in the next cell cycle (Arora et al., 2017, Barr et al., 2017, Naetar, Soundarapandian et al., 2014, Yang et al., 2017. Consistent with this second restriction point operating in cultured hNPCs, we observe a variable length G0/G1 with a portion of hNPCs entering a transient G0-like state each cell cycle. We further demonstrate that TP53, which regulates CDKN1A/p21 expression in hNPCs, is required for this transient G0-arrest. Moreover, we identify p53independent effectors that likely act through Hippo-YAP signaling pathway (NF2 and PTPN14) or via novel mechanisms (e.g., CREBBP and TAOK1) to affect expression of genes that allow bypass or escape of the G0-like state (e.g., MYC, HRAS, CCND1/2, CCNA2, CDKN1B/p27) (Figs 5C and S3).
Further, using scRNA-seq, network, and molecular analysis, we defined the gene expression signature associated with this G0-like state. This signature prominently featured neurodevelopmental genes and genes specifically expressed in quiescent NSCs in post-natal mice (Fig 6F), which were not found significantly expressed in other hNPC cell cycle populations. Since CREBBP, NF2, PTPN14, TAOK1, and TP53 KO both decrease residence in a G0-like state and disrupt differentiation specification, the results are consistent with Neural G0 being important for hNPC development potential, at least in vitro (Figs 7F,G and S7E-H).
Changes in the length of G1 can affect cell-fate specification. During cerebral development, fast dividing multipotent stem/progenitor cells that first populate the neocortex have a short G1, and, as lineage-commitment occurs, neuro-progenitors' cell cycle is slowed by lengthening G1 phase, eventually giving rise to post-mitotic, terminally differentiated cells (Miyama, Takahashi et al., 1997, Takahashi, Nowakowski et al., 1995. Experimentally modulating G1  Lastly, our results also suggest that G0-skipping could play unexpected roles in developing neoplasms, as the G0-skip genes are among the most highly mutated across 35 different tumor types ( Fig 1D). In addition to faster proliferation, if the transient G0-like state has dual roles in maintaining developmental potential and allowing time to repair DNA damage, then G0-skipping could simultaneously reduce the barrier to adopt alternative developmental states and promote genomic instability. In glioblastoma patient isolates, we observe phenotypes consistent with G0-skipping, including uniform, fast G1 transit times and diminished Neural G0 gene expression in G0/G1 populations (Figs 3D and S7B).
Collectively, this work reveals novel aspects of the proliferation of cultured hNPC through the use of functional genetics and single-and subpopulation-resolved genomic profiling, which may further aid in understanding NPC biology. These data are also a novel resource for the study of CREBBP, NF2, PTPN14, TAOK1, and TP53 and the dynamics of stem cell proliferation, including a detailed transcriptional model of the hNPC cell cycle.

Lentiviral production
For virus production, lentiCRISPR v2 plasmids (Sanjana et al., 2014) were transfected using polyethylenimine (Polysciences) into 293T cells along with psPAX and pMD2.G packaging plasmids (Addgene). To produce lentivirus for the whole-genome CRISPR-Cas9 libraries, 25x150mm plates of 293T cells were seeded at ~15 million cells per plate. Fresh media was added 24 hours later and viral supernatant harvested 24 and 48 hours after that. For screening, virus was concentrated 1000x following ultracentrifugation at 6800xg for 20 hours. For validation, lentivirus was used unconcentrated at an MOI<1.

CRISPR-Cas9 screening
For large-scale transduction, NPC cells were plated into T225 flasks at an appropriate density such that each replicate had 250-500-fold representation, using the two previously published CRISPR-Cas9 libraries (Doench, Fusi et al., 2016, Shalem, Sanjana et al., 2014 (Addgene) or a custom synthesized sgRNA library (Twist Biosciences) targeting 1377 genes derived from (Toledo et al., 2015). Full lists of the guides are in Table S3. NPCs and GSCs were infected at MOI <1 for all cell lines. Cells were infected for 48 hours followed by selection with 1-2 µg/mL (depending on the target cell type) of puromycin for 3 days. Post-selection, a portion of cells were harvested as Day 0 time point. The remaining cells were then passaged in T225 flasks maintaining 250-500-fold representation and cultured for an additional 21-23 days (~10-15 cell doublings) or 10 days. Genomic DNA was extracted using QiaAmp Blood Purification Mini or Midi kit (Qiagen).
A two-step PCR procedure was performed to amplify sgRNA sequence. For the first PCR, DNA was extracted from the number of cells equivalent to 250-500-fold representation (screen-dependent) for each replicate (2-4 replicates) and the entire sample was amplified for the guide region. For each sample, ~100 separate PCR reactions (library and representation dependent) were performed with 1 μg genomic DNA in each reaction using Herculase II Fusion DNA Polymerase (Agilent) or Phusion High-Fidelity DNA Polymerase (Thermo Fisher). Afterwards, a set of second PCRs was performed to add on Illumina adaptors and to barcode samples, using 10-20ul of the product from the first PCR. Primer sequences are in Table S8. We used a primer set to include both a variable 1-6 bp sequence to increase library complexity and 6 bp Illumina barcodes for multiplexing of different biological samples. The whole amplification was carried out with 12 cycles for the first PCR and 18 cycles for the second PCR to maintain linear amplification. Resulting amplicons from the second PCR were column purified using Monarch PCR & DNA Cleanup Kit (New England Biolabs; NEB) to remove genomic DNA and first round PCR product. Primers sequences are available in Table S8. Purified products were quantified (Qubit 2.0 Fluorometer; Fisher), mixed, and sequenced using HiSeq 2500 (Illumina).
Bowtie was used to align the sequenced reads to the guides (Langmead, Trapnell et al., 2009). The R/Bioconductor package edgeR was used to assess changes across various groups (Robinson, McCarthy et al., 2010). For the tiling library, only guides that mapped once to the genome and are within the gene's coding region were considered for further analysis. Raw and mapped data files are available at the Gene Expression Omnibus database (GSE117004).

Individual lentiviral-sgRNA assembly for validation
For retests, individual or pooled sgRNA were cloned into lentiCRISPR v2 plasmid.
Briefly, DNA oligonucleotides were synthesized with sgRNA sequence flanked by the following: 5': tatatcttGTGGAAAGGACGAAACACCg 3': gttttagagctaGAAAtagcaagttaa PCR was then performed with the ArrayF and ArrayR primers (Table S8). The PCR product was gel purified using the ZymoClean Gel DNA recovery kit (Zymo Research).
Gibson Assembly Master Mix (NEB) was used to clone the PCR product into lentiCRISPR v2 plasmid (Sanjana, Shalem et al., 2014). The ligated plasmid was then transformed into Stellar Competent cells (Clontech), and streaked onto LB agar plates.
The resulting clones were grown up and sequence verified (GeneWiz).

Viability and Proliferation Assays
Cells were infected with lentiviral gene pools containing 3-4 sgRNAs per gene or with lentivirus containing a single sgRNA to the respective gene (Table S8). Initial cell density was carefully controlled for in each experiment by counting cells using a Nucleocounter NC-100 (Eppendorf) and cells were always grown in subconfluent conditions. For viability assays, following selection, cells were outgrown for 7-10 days, then harvested, counted, and plated in triplicate onto 96-well plates coated with laminin in dilution format starting at 1,000 cells to 3,750 cells per well (cell density depended on cell isolate and duration of assay). Cells were fed with fresh medium every 3-4 days.
After 7-12 days under standard growth conditions, cell proliferative rates were measured using Alamar blue reagent according to manufacturer's instructions (Invitrogen). For analysis, sgRNA-containing samples were normalized to their respective nontargeting control (NTC) samples. For doubling time assays, cells infected with individual sgRNAs or NTC were routinely cultured (split every 3-5 days), and counted at each split (Nucleocounter NC-100; Eppendorf). The overall growth of each well containing an individual sgRNA was calculated and compared to the NTC well.
Comparisons between multiple experiments were normalized.

Competition experiment
NPCs were infected with lentiviral gene pools containing 3-4 sgRNAs per gene, puromycin selected, and mixed with NPCs infected with lentiviruses containing turboGFP at an approximate 1:9 ratio, respectively. Cultures were outgrown for 23 to 31 days and flow analysis (FACS Canto; Becton Dickinson) was conducted every 7-8 days for GFP expression. Flow analysis data was analyzed using FlowJo software. For each sample, the GFP-population for each time point was normalized to its respective Day 0 GFP-population and the NTC (competition index).

Time-lapse microscopy
NPCs were infected with lentiviral gene pools containing 3-4 sgRNAs per gene or with individual sgRNAs, puromycin selected, outgrown for >13 days, and plated onto 96-well plates or 24-well plates. Plates were then inserted into the IncuCyte ZOOM (Essen BioScience), which was in an incubator set to normal culture conditions (37° and 5% CO 2 ), and analyzed with its software. For the cell confluency experiment, phase images were taken every hour for 72 hours. For the FUCCI cell cycle experiment, images were taken every 10-15 minutes for 72-120 hours. Cell cycle transit time for G0/G1 (mCherry-CDT1(aa30-120)+) and S/G2/M (mAG-Geminin(aa1-110)+) was manually scored by three different observers in actively dividing cells (those that could be followed from mitosis to mitosis). Each KO was scored by at least 2 independent observers and consistency between scorers was checked through shared analysis of a standard.

Western blotting
Cells were harvested, washed with PBS, and either immediately lysed or snap-frozen and stored at -80˚C until lysis. Cells were lysed with modified RIPA buffer (150mM NaCl, 50mM Tris, pH 7.5, 2mM MgCl 2 , 0.1% SDS, 2mM DDT, 0.4% deoxycholate, 0.4% Triton X-100, 1X complete protease inhibitor cocktail (complete Mini EDTA-free, Roche) and 1U/µL benzonase nuclease (Novagen) at room temperature for 15 minutes. Cell lysates were quantified using Pierce 660nm protein assay reagent and proteins were loaded onto SDS-PAGE for western blot. The Trans-Blot Turbo transfer system (Bio-Rad) was used according to the manufacturer's instructions. See Table S8 for antibodies and dilutions. An Odyssey infrared imaging system was used to visualize blots (LI-COR) following the manufacturer's instructions.
Cells were grown out for 21 days with splitting every 3-4 days and maintaining equivalent densities. Cells were counted (Nucleocounter NC-100; Eppendorf) and plated 3 days before analysis on an LSR II (BD). Controls cultured in the same conditions included cells transduced with guides against 3 non-growth limiting genes, including GNAS1, and showed equivalent FUCCI ratios. Results were analyzed using FlowJo software.

Immunofluorescence and CDK2 Activity
U5-NPCs were plated on acid-washed glass coverslips (phosphorylated Rb and CDK2 activity) or 96-well imaging plates (differentiation; Corning). They were fixed overnight in 2% paraformaldehyde (USB) at 4˚C, washed with DPBS (with calcium and magnesium) (Fisher), and blocked and permeabilized with 5% goat serum (Millipore), 1% bovine serum albumin (Sigma), and 0.1% triton X-100 (Fisher) in DPBS for 45 minutes at room temperature. Samples were stained with primary antibody diluted in 5% goat serum in DPBS overnight at 4˚C, washed with DPBS, and stained with secondary antibody (diluted 1:200 in 5% goat serum in DPBS) at 37˚C for 45 minutes. See Table S8 for antibodies and dilutions. Samples were washed with DPBS, dyed with 100 ng/mL 4',6-diamidino-2-phenylindole (DAPI) diluted in DPBS for 20 minutes at room temperature, and washed with DPBS. Coverslips were preserved using ProLong Gold Antifade Mountant (Thermo Fisher) and inverted on glass slides. For differentiation, images were acquired on Nikon Eclipse Ti using NIS-Elements software (Nikon). Global/Otsu thresholding, and distinguishing clumped objects by shape). CDK2 activity was defined by the cytoplasmic to nuclear ratio of the mVenus-DHB reporter, with the cytoplasmic intensity of the DHB reporter defined as the upper quartile intensity of a 2pixel ring around the CDT1-defined nucleus due to the irregular shape of the U5-NPCs.

p27 reporter
The p27 reporter was constructed after (Oki et al., 2014), using a p27 allele that harbors two amino acid substitutions (F62A and F64A) that block binding to Cyclin/CDK complexes but do not interfere with its cell cycle-dependent proteolysis. This p27Kallele was fused to mVenus to create p27K --mVenus. To this end, the p27 allele and mVenus were synthesized as gBlocks (IDT) and cloned via Gibson assembly (NEB) into a modified pGIPz lentiviral expression vector (Open Biosystems). Lentivirally transduced cells were puromycin selected and validated using mCherry-CDT1 FUCCI and HDAC inhibitor treatment (48 hours of 5 µM apicidin (Cayman)) to induce G0/G1 arrest using FACS (LSR II from Becton Dickinson and FlowJo software).

Bulk RNA sequencing expression analysis
For G0/G1 NPC, cells singly positive for mCherry-CDT1 FUCCI were sorted on a FACSAria II (BD) directly into TRIzol reagent (Life Technologies). For differentiating cells, cells were sparsely plated and cultured with growth medium without EGF or FGF-2 for 7 days before being lysed with TRIzol reagent. For both, 2 replicates per condition were harvested. RNA was extracted using Direct-zol RNA MiniPrep Plus (Zymo Research). Total RNA integrity was checked and quantified using a 2200 TapeStation (Agilent). RNA-seq libraries were prepared using the KAPA Stranded mRNA-seq Kit with mRNA capture beads (KAPA Biosystems) according to the manufacturer's guidelines. Library size distributions were validated using a 2200 TapeStation (Agilent).
Additional library QC, blending of pooled indexed libraries, and cluster optimization was performed using the Qubit 2.0 Fluorometer (Thermo Fisher). RNA-seq libraries were pooled and sequencing was performed using an Illumina HiSeq 2500 in Rapid Run mode employing a paired-end, 50 base read length sequencing strategy.

Bulk RNA sequencing data analysis
RNA-seq reads were aligned to the UCSC mm10 assembly using Tophat2 (Trapnell, Roberts et al., 2012) and counted for gene associations against the UCSC genes database with HTSeq (Anders, Pyl et al., 2015). Differential expression analysis was performed using R/Bioconductor package edgeR (Robinson et al., 2010). Samples for G0/G1 bulk RNA-seq were collected in two batches, so batch-dependent genes were removed before analysis (inter-batch p-value<0.01 by Wilcoxon-Mann-Whitney). To ensure that no genes were eliminated that may be regulated specific to a particular knockout, genes with a CPM variability greater than 2-fold compared to the internal batch control and an expression greater than 1 CPM in at least one sample were retained. Differentially expressed genes (DEG) at the transcription level were found using a statistical cutoff of FDR < 0.05 and visualized using R/Bioconductor package pheatmap. Kolmogorov-Smirnov test were conducted in R using the function ks.test from stats package. Raw sequencing data and read count per gene data can be accessed at the NCBI Gene Expression Omnibus (GSE117004).

Gene ontology analysis
Gene Ontology (GO)-based enrichment tests were implemented using GOseq (v 1.23.0) (Young, Wakefield et al., 2010), which corrects for gene length bias. Gene lists were also analyzed for pathways using the R/Bioconductor package ReactomePA (v 1.15.4) (Yu & He, 2016). Analysis used all genes either up or down-regulated with a FDR<0.05 compared to NTC. GO terms with adjusted P-values<0.05 were considered significantly enriched. Venn diagrams were generated on http://bioinformatics.psb.ugent.be/webtools/Venn/.

Single cell RNA-sequencing Sample Preparation
Single cell RNA-sequencing was performed using 10x Genomics' reagents, instruments, and protocols. Single cell RNA-Seq libraries were prepared using GemCode Single Cell 3' Gel Bead and Library Kit. FUCCI U5-NPCs (both with and without lentiviral TAOK1 KO, >14 days outgrowth) were harvested and half the cells were sorted using the

scRNA-seq Analysis
CellRanger (10x Genomics) was used to align, quantify, and provide basic quality control metrics for the scRNA-seq data. Using Seurat version 2.2.1, the scRNA-seq data from wild-type U5 cells and sgTAOK1 knock-out cells were merged and analyzed.
Both scRNA-seq data were loaded as counts, normalized, and then scaled while taking into account both percent of mitochondria and the number of UMIs per cell as covariates. The union of the top 1,000 most variant genes from each dataset were used in canonical correlation analysis (CCA) to merge the two datasets via alignment of their subspace. We then identified clusters of cells using a shared nearest neighbor (SNN) modularity optimization-based clustering algorithm. Marker genes for each cluster were identified as differentially expressed genes, and the determination of 8 clusters was based on the discovery of strong markers for 6 of the eight clusters (both the G1 and low RNA clusters did not have significantly upregulated marker genes). Identity of clusters was determined primarily through the expression of cyclins and cyclindependent kinases, and secondarily through the function of other marker genes. A tSNE visualization was generated with a perplexity setting of 23.
Network analysis was used to determine the trajectories of cells through the cell cycle. First, the cluster centroids (mean expression for each gene across all the cells from a cluster) were used to compute the Canberra distance measure. In a cycle like a cell cycle, it is expected that on average there would be 2 edges between each cell cycle state. A distance cutoff of 240 led to 2.28 connections per cluster was used to turn the distance matrix into a network.
Network analysis of the clusters was performed using the STRING database (Szklarczyk, Morris et al., 2017) and visualized using Cytoscape software. Transcription factors were identified according to TFcheckpoint (Chawla, Tripathi et al., 2013).

Hypergeometric Analysis and Representation Factor Calculations
Hypergeometric tests (Johnson, Kemp et al., 2005) were carried out in R using function phyper (stat.ethz.ch/R-manual/R-devel/library/stats/html/Hypergeometric.html). Gene lists were prefiltered for the shared genes in each analysis to get the total gene population size, (i.e., 2739 genes for single cell analysis that had greater than 3 counts per cell in at least 10 cells and removing batch-effected genes for G0/G1 bulk RNAsequencing).
Representation factors were calculated according to (Kim, Lund et al., 2001). The representation factor shows whether genes from one list (list A) are enriched in another list (list B), assuming that genes behave independently.

Author contributions
Project conception and design was carried out by P.J.P., H.M.F., and C.M.T., CRISPR-    The data are presented as the mean ± standard deviation (SD). Significance was assessed using a two-tailed student's t test.

Figure 2
N T C T P 5 3 -2 N F 2 -2 T P 5 3 -3 T A O K 1 -3 T P 5 3 -1 P T P N 1 4 -2 The data are presented as the mean ± SD. Significance was assessed using a twotailed student's t test (C) or Mann-Whitney test (D).    White dots indicate FDR<0.05.

D,E Overlap of a selection of biological process enriched in up (D) and down (E)
regulated genes in G0/G1 following KO of G0-skip genes using gene ontology analysis.

Figure S1 -CRISPR-Cas9 Gene Knockout Screens Identify hNPC Candidate Antiproliferative Genes; Effect is Diminished in Glioblastoma Tumor Isolates.
A Overlap of screen hits between the four CRISPR-Cas9 screens. For all screens, hits were defined as genes having multiple guides with FDR<.01 and z-score>2. For Brunello library, there were additional requirements of at least one guide with z-score>3 and at least 20 reads at Day 0. CREBBP, NF2, PTPN14, TAOK1, and TP53 were defined as hits in every screen. Specificity to one of the NPC isolates was only defined if the same library was tested in both isolates. B Overlap of the candidate antiproliferative genes with validated tumor suppressors (Futreal et al., 2004) (Tier 1 and 2; Cancer Gene Census v84). Significance assessed using hypergeometric analysis. C Scatterplot of guides in U5-NPCs, GSC-0131 (n=5), GSC-0827 (n=5), and GSC-1502 (n=5) with a positive log 2 FC (logFC) following ~3 weeks outgrowth post-selection following transduction with a custom CRISPR-Cas9 retest library. Dots for individual guides were colored only if significant (sig) (FDR<0.05 and logFC>1). Loss of the proliferation-limiting genes has a more limited effect in GSCs compared to NPCs in both magnitude and number of significant guides. Intersection Size   A Hierarchical clustering and percent overlap of differentially expressed up-and down-regulated genes between G0-skip KO and NTC G0/G1 U5-NPCs. Percent gene overlap is relative to the smaller gene set. B

A
Number of unique molecular identifiers (nUMI) per cell in each single cell RNA-seq cluster. The Low RNA cluster has dramatically less nUMI than any other cluster, making it difficult to distinguish actual cells belonging to this cluster from background. Therefore, this cluster was excluded from further analysis. The nUMI increases as the cell cycle progresses, with G2/M cells having the greatest nUMI. Significance assessed using a two-tailed student's t-test.

B
Comparison of cluster size from scRNA-seq for the entire cell population of U5-NPCs compared to a subpopulation sorted for G0/G1 (mCherry-CDT1(aa30-120)-positive). The clusters with an increase in representation were characterized as G0/G1. Clusters that remained the same or moderately decreased in representation are considered transition states with some of the cells retaining or initiating CDT1 accumulation. Clusters that decreased in representation greater than 2-fold were considered S/G2/M and mCherry-CDT1-negative. C Dot plot for some of the key markers for each cluster progressing through cell cycle. Color is according to average expression (Avg. Exp.) per cell, while size of the dots corresponds to the proportion of cells within the cluster in which the gene was sequenced (% Exp.). Red text indicates cyclin and cyclin-dependent kinase genes that are differentially expressed in the different clusters. The genes that are depleted in each cluster supports the G0/G1 classification of the cells. For example, CDK1 is significantly reduced in Neural G0, G1, Late G1, and early G1/M, but enriched in S, S/G2, and G2/M. Similar depletion of other cell cycle genes is seen, including MKI67 (depleted in Neural G0, late G1, and enriched in G2/M and S/G2, AURKA (depleted in Neural G0, G1, S, and M/early G1 and enriched in G2/M), CCNA2 (depleted in Neural G0, G1, late G1, and M/early G1 and enriched in S/G2 and G2/M), and TOP2A (depleted in Neural G0, G1, late G1, and M/early G1 and enriched in S/G2 and G2/M). In addition, Neural G0 showed the most significant repression of these genes peaking in other phases of cell cycle, including suppression of CCND1 expression.

D
Protein-protein interaction networks based on the STRING database (Szklarczyk et al., 2017) for the enriched genes that define the G2/M, S, and S/G2 cluster after single-cell RNA-sequencing in unsorted U5-hNPCs. The thickness of the connecting lines corresponds to the strength of the predicted interaction based on experimental data, curated databases, co-expression analysis, and gene orthology. The network visualization only includes genes that connect with the rest of the network. Genes from enriched gene ontology groups are indicated by colored nodes, including mitotic cell division (red), DNA replication (purple), DNA conformation change (yellow), and RNA-related genes from various gene ontology biological processes (green).    A Protein-protein interaction networks based on the STRING database (Szklarczyk et al., 2017) for the enriched genes that define the G0/G1 clusters: Neural G0, Early G1/M and Late G1. The network visualization only includes genes that connect with the rest of the network. G1 only showed enriched interactions between IGFR1 signaling genes IGFBP3 and IGFBP5 (not visualized). Neural G0 has a more diffuse network with only a few interactions per gene and includes a number of transcription factors that may be important in mediating this cell cycle state (pink). The late G1 cluster has two gene groups: the cholesterol synthesis group (orange) and YAP-NF2 targets (yellow). M/early G1 retains residual expression or non-degradation of a subset of genes enriched in G2/M, including the mitosis-peaking genes CCNB2, CDKN3, and BIRC5 (red) as well as spliceosomal genes (green). B Overlap of conserved YAP target genes up-regulated following G0-skip knockout of NF2 (yellow), TAOK1 (red), or both (purple) (see Fig S3C) with genes defining cell cycle transcriptional clusters. The NF2-dependent YAP targets significantly overlap with the Late G1 cluster, while the TAOK1-dependent genes overlap with the S/G2 and G2/M clusters. Significance assessed using hypergeometric analysis. C Proportion of mouse embryonic stem cells (mESCs) (Buettner et al., 2015) classified to each hNPC transcriptional cluster following sorting for cell cycle by Hoechst staining for amount of DNA. D Significance of overlap of the neurogenesis cell subpopulation-defining genes in the early neurogenic lineage from two murine single cell RNA-sequencing studies (Artegiani et al., 2017, Llorens-Bobadilla et al., 2015 compared to single cell cluster definitions (up-regulated genes) from unsorted U5-hNPCs grown in culture. Clusters presented in order of increasing activation with quiescent neural stem cells on the left and proliferating progenitors on the right. Significance assessed though hypergeometric analysis. RF = representation factor. NSC = neural stem cell. OPC = oligodendrocyte progenitor. q = quiescent. a = activated.

YAP-targets (NF2) Transcription Factors
Significant, but not in network:

DDAH1 GADD45B
Significant, but not in network: SOX6 TSC22D1 HEY1 A,B Overlap of genes that are significantly down (A) or up (B) following G0-skip KO in G0/G1 U5-NPCs compared to the Neural G0 (A) or Late G1 (B) single cell RNA-sequencing clusters. C,D Significance of overlap of the down (C) and up (D) regulated genes from bulk RNAsequencing of G0/G1 sorted 0131 and 0827 glioblastoma stem cell (GSC) isolates with the single cell cluster definitions (up-regulated genes). Significance assessed though hypergeometric analysis. RF = representation factor. E Immunocytochemistry of U5-NPCs in normal NPC culture conditions and after 1 week removal of growth factors (GF) (EGF and FGF withdrawal). Staining is for various neural lineage markers including the neural stem/progenitor markers Nestin (green) and Sox2 (red), the glial marker GFAP (red), and the neuronal marker beta-tubulin III (TUJ1; green). Scale bar = 50 µm. F Dendrogram of unbiased hierarchical clustering of gene expression from U5-NPCs before and after 1 week of growth factor withdrawal (controls) and after KO of G0-skip genes. G Venn diagram for brain-specific transcription factors (Gray et al., 2004) that are significantly reduced (FDR<0.05) in G0-skip mutants compared to NTC following 1 week differentiation through growth factor-withdrawal. Bold indicates the transcription factor that is reduced in all G0-skip mutants, FOXJ1. H Gene expression of GFAP (glial marker) and DRC1 (cilia component) before (NPC) and after 1 week of mitogenic GF withdrawal (diff) for control cells and after KO of G0-skip genes. These genes are significantly induced during differentiation and reduced in the G0-skip mutants to varying degrees.