Main

Human radial glia (RG) evolved an increased proliferative capacity, altered cell fate potential and protracted period of maturation, supporting the increased number and complexity of daughter cells2,15,16. Gene regulatory networks governing RG self-renewal, differentiation and maturation have long been implicated in cortical expansion17,18. RG sequentially produce distinct subtypes of neuron followed by glial cell types19,20. Although cortical inhibitory neurons (INs) are generated mainly in the ganglionic eminences and migrate to the cortex21,22,23, recent lineage tracing studies1,24,25,26,27 and developmental cell atlases9,12 have demonstrated that human RG also produce cortical-like INs at late stages of neurogenesis that share transcriptional signatures with caudal ganglionic eminence (CGE) and lateral ganglionic eminence derived INs. However, the role of TFs in regulating RG differentiation decisions, lineage plasticity to produce INs and maturation remains largely unexplored.

Genetic perturbations combined with measurements of single-cell gene expression provide a powerful approach, termed Perturb-seq, for high-throughput dissection of gene function13,14. Cas9-based CRISPR loss-of-function perturbation approaches have been applied recently to study gene regulatory networks influencing cortical development using induced pluripotent stem (iPS)-cell-derived organoid28,29,30 and mouse31,32 models. However, Cas9-induced double-stranded DNA breaks can cause cytotoxicity, influencing proliferation and differentiation decisions33,34, whereas acquired genetic and epigenetic variation in source iPS cells35,36, and patterning biases37 and cell stress38 during differentiation, can influence cell-type fidelity and fate specification in organoid models. CRISPR interference (CRISPRi) screens using dCas9-KRAB enable efficient and uniform repression of target genes with limited cellular toxicity39,40, whereas primary cortical culture41,42 captures physiological aspects of human cortical neurogenesis, supporting sensitive detection of cell fate choices.

Here we established a primary cell model system that recapitulates in vivo differentiation dynamics and performed Perturb-seq to measure the transcriptional and cell fate consequences of repressing 44 TFs that are expressed robustly in the human cortical RG lineage. Extending human primary cell culture approaches1,42,43, we targeted TFs in a homogeneous RG population and then removed growth factors to permit cortical neurogenesis, cell fate choice and early subtype specification. Our screening revealed the role of ZNF219—not previously described in cortical development—in repressing neuronal differentiation, opposing roles for NR2E1 and ARX in regulating the balance of human excitatory versus inhibitory neurogenesis, and the role of ARX in safeguarding IN subtype specification through transcriptional repression of downstream transcription cofactor LMO1. Intersecting dysregulated genes under different perturbations revealed candidate hub effector genes downstream of several TFs enriched for roles in neurodevelopmental and neuropsychiatric disorders. Coupling CRISPRi screening with barcoded lineage tracing demonstrated the potential to engineer lineage plasticity and developmental tempo of individual RG through TF perturbation. Together with parallel screening in rhesus macaque, our data illuminate conserved mechanisms governing cortical RG lineage progression across primates.

Systematic TF repression during human corticogenesis

We designed a primary cell model of neurogenesis and lineage progression to evaluate the impacts of TF repression on cell fate choice during human cortical development. To direct gene targeting to RG at the start of differentiation, we first enriched for RG isolated from primary human tissue samples by adding epidermal growth factor (EGF) and fibroblast growth factor 2 for 5 days before infecting with an all-in-one CRISPRi lentivirus (Fig. 1a), and we further expanded RG for 1 week allowing for target gene knockdown (KD) before differentiation39 (Extended Data Fig. 1a–c). We then replaced growth factors with brain-derived neurotrophic factor (BDNF) to support spontaneous differentiation of perturbed RG. Consistent with recent studies1,42,43, this model recapitulates in vivo RG lineage progression and generation of excitatory neurons (ENs) and INs (Fig. 1a and Extended Data Fig. 1a,b), enabling detection of perturbations that affect differentiation dynamics and cell fate choice. Repression of PAX6 and EOMES recapitulated the effects of these TFs described in other model systems in promoting excitatory neurogenesis, highlighting the potential of our system to uncover additional regulators44,45,46,47 (Extended Data Fig. 1d).

Fig. 1: Single-cell Perturb-seq on 44 TFs in primary human neuronal differentiation.

a, Experimental design for high-throughput perturbation of 44 TFs in primary human cortical progenitors from four individuals, with immunocytochemical labelling of EOMES and NEUROD2 in in vitro 2D culture of human cortical RG before (day 0) and after (day 3, day 7) induced differentiation. Representative images from four individuals are shown. FACS, fluorescence-activated cell sorting. b, Uniform manifold approximation and projections (UMAPs) of cells collected on day 0 (21,151 cells, n = 2 individuals) and day 7 (116,166 cells, n = 4 individuals), coloured by cell class, individual and sex. c, UMAPs highlighting cells from different timepoints. d, UMAP coloured by supervised cell type, with stacked barplots (bottom) showing cell-type distributions across individual at each timepoint. e, UMAPs showing expression of HOPX, MKI67, NEUROD2, DLX2 and BCAS1. f, Left, dotplot of marker gene expression across annotated clusters; right, barplot of cell numbers. Dot size denotes the expressing-cell fraction and colour denotes mean expression. g, Heatmap of normalized expression of 44 TFs targeted in this study along pseudotime with branches for IN and EN (left) lineages, with a tree plot showing cell-type distribution along pseudotime (right). Pseudotime was calculated using NT cells on day 7. h, Stacked barplot showing sgRNA distributions per TF target on day 7, coloured by individual sgRNA. i, Dotplot showing target gene KD efficiency compared with NT cells on day 7, filled by individual sgRNA with borders highlighting significance. log2FC values were calculated with DEseq2 in each cell class in n = 4 individuals and the lowest value for each sgRNA was used for visualization and filtering. sgRNAs with less than 25% reduction were removed from downstream analyses; remaining active sgRNAs showed a median 72% KD (28% residual expression). Astro, astrocytes; DIV, dividing; ImN, immature neuron; MG, microglia; OPC, oligodendrocyte progenitor cell; Vasc, vascular. Scale bars, 200 μm. Panel a created in BioRender. Nowakowski, T. (https://BioRender.com/7cr6o12).

To systematically identify regulators of human cortical neurogenesis, we first prioritized TFs using single-cell RNA sequencing (scRNA-seq) and single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq) multiome data10,12. We selected 44 TFs based on robust expression, motif accessibility, gene regulatory network size, target gene expression and predicted transcriptional consequences in the RG lineage (Supplementary Table 1; Methods). We then adapted an all-in-one CRISPRi vector co-expressing green fluorescent protein (GFP)48 by inserting a capture sequence in the single guide RNA (sgRNA) scaffold, enabling direct capture of sgRNAs during scRNA-seq using a 10x Genomics platform49, thereby supporting screening in primary cell culture systems50. We synthesized and cloned a library containing 164 sgRNAs, targeting the 44 TFs, and included 20 non-targeting control (NT) sgRNAs. For each TF, we targeted active promoters with accessible chromatin during human cortical neurogenesis9 using three sgRNAs per promoter51,52 (Supplementary Table 2). We derived primary human cultures from cryopreserved cortical tissue of four individuals at stages of peak neurogenesis from gestational weeks (GW) 16–18 (ref. 41). We delivered the CRISPRi library by lentivirus into primary RG, targeting less than 30% infection rate to generate libraries with mostly singly infected cells (Extended Data Fig. 1a and Supplementary Table 3), removed growth factors, and confirmed efficient gene repression before and throughout differentiation (Extended Data Fig. 1c,e).

scRNA-seq confirmed a 95% population of cycling RG on day 0 of differentiation, marked by co-expression of RG marker HOPX and proliferation marker MKI67 (Fig. 1b–f). By day 7, three principal cortical cell class trajectories emerged—EN, IN and oligodendrocyte lineages marked by NEUROD2, DLX2 and BCAS1, respectively (Fig. 1c–e)—in addition to a continuum from RG to astrocytes marked by AQP4 (Extended Data Fig. 1f). Integrating both experimental timepoints highlighted homogeneity of day 0 populations with minimal spontaneous differentiation (Fig. 1c and Extended Data Fig. 1g), as well as comparable representation of cells derived from independent technical and biological replicates across cell types with minimal batch effects (Fig. 1b,d and Extended Data Fig. 1h).

Pearson correlation and reference mapping to a developmental cortical cell atlas12 confirmed the recapitulation of in vivo-like gene expression, cell types, states, differentiation dynamics and cell fate choice in the primary culture differentiation system (Extended Data Fig. 1i,j). Most neurons exhibited immature states (Extended Data Fig. 1j), consistent with recent generation from RG after differentiation. These differentiation dynamics were also recovered by RNA velocity and pseudotime analysis (Fig. 1g and Extended Data Fig. 1l,m). Notably, comparison with iPS-cell-based models revealed a near twofold increase in the Pearson correlation coefficients of marker gene expression compared with in vivo cell types in the primary two-dimensional (2D) system, supporting improved fidelity to normal development (Extended Data Fig. 1i,k). In addition, primary culture in 2D and organotypic slice both showed reduced transcriptional signatures of cellular stress across principle dimensions, including glycolysis, endoplasmic reticulum stress, oxidative stress and apoptosis53,54,55 (Extended Data Fig. 1n), supporting the physiological relevance of the primary cell models.

We assigned sgRNAs to 118,456 cells (Methods) across timepoints, with a mean of 200 singly infected cells per sgRNA on differentiation day 7 (Fig. 1h and Extended Data Fig. 1o). KD efficiency was calculated in each cell class and 18 sgRNAs with less than 25% KD (log2 fold change (FC) > −0.4) were removed from downstream analyses (Fig. 1i and Supplementary Table 4; Methods). Remaining active sgRNAs exhibited a median KD efficiency of 72%, with comparable transcriptional responses to independent sgRNAs targeting the same promoter (Fig. 1i and Extended Data Fig. 2). For further analysis, we collapsed all active sgRNAs sharing the same target TF, which yielded a median KD efficiency of 80% and a mean of 600 cells per gene (Extended Data Fig. 1o,p).

Transcriptional regulators of human corticogenesis

We next examined the consequences of TF repression on gene expression and cell-type composition by differentiation day 7 (Fig. 2a and Supplementary Table 5; Methods). Gene expression changes were correlated between individual sgRNAs and genes (average Pearson R = 0.88) and cell-type abundance was preserved upon downsampling of cell number (average Pearson R = 0.94), supporting the power of the screen to detect changes in both modalities (Extended Data Fig. 3a,b).

Fig. 2: Convergent transcriptional changes reveal neuropsychiatric-disorder-enriched hub effector genes.

a, Schematics showing analytical workflow applied to the day 7 Perturb-seq dataset obtained from Fig. 1. b, Scatterplot for summed |estimated coefficient| in cell abundance (x axis) and number of DEGs (y axis) across seven main cell types in the EN and IN lineage for each TF perturbation. Dot size reflects perturbations with the strongest phenotypes calculated by (summed |estimated coefficient| × total number of DEGs), with embedded pie charts indicating fractions of up- and downregulated DEGs. c, Volcano plot of DEGs shared across different numbers of perturbations. DEGs with |log2FC | > 0.2 are shown. DEGs that overlap with seven sets of neuropsychiatric disorder-associated gene sets are coloured red and SCZ-associated genes are labelled. d, Dotplot showing gene ontology biological processes enrichment among convergent DEGs that were detected in at least two TF perturbations, using non-convergent DEGs as background. Dot size denotes Z score; colour denotes −log10[adjusted P]. e, Dotplot showing enrichment of seven disease-associated gene sets (Methods) among convergent DEGs versus non-convergent DEGs. One-sided Fisher’s exact test; dot size denotes odds ratio; colour denotes −log10[P]. f, Line plot showing odds ratio of disease-associated gene enrichment in DEGs shared across the minimum number of perturbations indicated. g, Regulatory network plot showing eight prioritized TFs (NR2E1,* ARX*, SOX2, ZNF219, NEUROD2, PHF21A,* SOX9* and CTCF) and connected convergent DEGs. Edge colour indicates direction of regulation (pink, positive regulation (where DEGs were downregulated upon TF perturbation)). Pink nodes or labels denote disease-associated genes or TFs. AD, Alzheimer’s disease; ADHD, attention-deficit/hyperactivity disorder; BP, bipolar disorder; NDD, neurodevelopmental disorder. Panel** a** created in BioRender. Nowakowski, T. (https://BioRender.com/7cr6o12).

We observed a positive correlation between the effects of individual TFs on cell-type composition and on gene expression across cell types (Fig. 2b) and within cell classes (Extended Data Fig. 3c). Comparing the extent of changes in both dimensions prioritized TFs whose depletion caused the strongest phenotypes: NR2E1, ARX, ZNF219, SOX2, SOX9, CTCF, NEUROD2 and PHF21A (Fig. 2b). These TFs also showed the strongest impact on other molecular phenotypes, including Euclidean distance and energy distance56 to NT, which compare the average expression between groups and the distance between and within groups, respectively, as well as maximum composition change among all cell types in the RG lineage (Extended Data Fig. 3d). Notably, the TFs with strong cellular and transcriptional phenotypes have been implicated previously in neurological disorders: ARX in X-linked lissencephaly57, epilepsy58,59 and intellectual disability60,61, and autism spectrum disorder (ASD)58; NR2E1 in schizophrenia (SCZ)62; SOX2 in intellectual disability and epilepsy63,64; CTCF in intellectual disability with microcephaly65,66; NEUROD2 in intellectual disability67 and early infantile epileptic encephalopathy68; PHF21A in intellectual disability with epilepsy and ASD69 and ZNF219 in a case of low IQ ASD70. Although we cannot rule out that additional TFs may impact differentiation or maturation phenotypes not captured in this assay due to incomplete KD, redundancy or stage-specific roles, these results highlight the functional importance and disease relevance of TFs prioritized by screening.

At the gene expression level, the predominance of upregulated differentially expressed genes (DEGs) following NR2E1 (88%) and ARX (85%) KD was consistent with their role as transcriptional repressors71,72 (Fig. 2b). Approximately 25% of DEGs were affected by the perturbation of more than one TF (Fig. 2c). These convergent DEGs showed enrichment for cell adhesion- and synaptic development-related terms, suggesting that effector genes regulating neurogenesis, migration and maturation are highly interconnected across TF regulatory networks (Fig. 2d). Moreover, intersecting these convergent DEGs with seven sets of neurodevelopmental and neuropsychiatric disorder-related genes revealed significant overlaps with SCZ and major depressive disorder (MDD)-associated genes, in comparison with non-convergent DEGs (one-sided Fisher’s exact test, P < 0.05) (Fig. 2e,f). This overlap includes PTPRD, modulated by ARX, SOX2 and NEUROD2; and IL1RAPL1, modulated by SOX2, ZNF219, CTCF and TBR1 (Fig. 2c,g). Both PTPRD and IL1RAPL1 are associated with SCZ and depletion of both genes have been shown to induce aberrant neurogenesis, maturation and behaviours in mouse models73,74,75. This finding highlights strong connections among disease-related genes in TF-driven regulatory networks, indicating their potential roles as hub effector genes downstream of developmental TFs.

Distinct fate outcomes of RG upon TF perturbations

We further investigated the impact of TF perturbations on cell-type composition using cluster-free approaches76 to map cell fate changes with higher resolution along developmental trajectories (Fig. 3a). These approaches yielded consistent results with cluster-aware methods at the gene and individual sgRNA level and were robust to downsampling (Extended Data Figs. 3e,f and 4).

Fig. 3: TF perturbations alter EN/IN output of human cortical RG.

a, Neighbourhood graphs based on the UMAP in Fig. 1d showing results of differential abundance testing using Milo under perturbation of eight prioritized TFs at the gene level. NT cells were downsampled to match perturbed cell numbers. Node sizes denote cell numbers in each neighbourhood; edges depict the number of cells shared between neighbourhoods; node colour denotes log2FC (false discovery rate (FDR) = 0.1) when perturbed. b, Top, schematics illustrating RG of distinct fate biases; bottom, experimental design for combined TF perturbation and lineage tracing in cortical RG of nine human individuals. c, UMAP of lineage-resolved 129,003 human cells collected on differentiation day 7, coloured by supervised cell type. A median of 15,469 cells per TF perturbation were recovered after filtering. d, Upset plot of cell class compositions in clones containing cells from more than one class (multi-class clones). e, Heatmap for lineage coupling score matrix in main cell types. f, UMAPs showing cell fate distributions in clonal clusters identified using scLiTr89: RG-, EN-, EN/IN-, IN- and oligodendrocyte-biased clusters. g, Box plots showing distributions of clonal cluster fractions in eight human individuals in NR2E1, ARX and ZNF219 KD versus NT. Two-sided paired Wilcoxon tests; P = 0.0078, 0.64, 0.016 (NR2E1); 0.0078, 0.0078, 1 (ARX) and 0.38, 0.023, 0.0078 (ZNF219). Centre line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers. *P < 0.05; **P < 0.01; ***P < 0.001. Panel** b** created in BioRender. Nowakowski, T. (https://BioRender.com/7cr6o12).

Different perturbations elicited a range of composition consequences in differentiated cell types (Fig. 3a and Extended Data Fig. 3e,f). These phenotypes included examples consistent with mechanistic studies. For example, repression of SOX2 resulted in accumulation of RG77, whereas repression of the proneural factor NEUROD2 enriched for progenitors and immature ENs at the expense of mature ENs (Fig. 3a and Extended Data Fig. 3e). However, screening also revealed additional roles for previously characterized genes implicated in disease. For example, NR2E1 knockout was previously shown to drive premature neural differentiation in mouse models78,79, but repression in human RG specifically enriched for post-mitotic INs, in addition to depleting RG (Fig. 3a and Extended Data Fig. 3e). Similarly, ARX knockout was previously shown to affect IN differentiation and migration in mouse models80,81,82, but repression in human RG specifically enriched for the EN lineage at the expense of INs.

Our screen also uncovered disease-linked genes with composition phenotypes, including ZNF219 and PHF21A, that to our knowledge have yet to be examined functionally during cortical neurogenesis. Repression of ZNF219 mirrored the effects of repressing SOX9—a known ZNF219 interaction partner in chondrocytes83—by increasing the proportion of ENs, but showed different effects from SOX9 in dividing RG and INs (Fig. 3a and Extended Data Fig. 3e). Repression of PHF21A—a member of the BRAF/HDAC complex implicated in synapse formation84—led to an overall trend towards EN over IN differentiation (Fig. 3a and Extended Data Fig. 3e,f), matching the effects of targeting other ASD causal genes ARX and CTCF85,86.

Flow cytometry analysis of the top TF candidates, ARX,* NR2E1* and ZNF219, with finer temporal resolution validated the cellular phenotypes observed in our Perturb-seq (Extended Data Fig. 5 and Supplementary Table 6), and further revealed early manifestations of compositional changes in KI67+ and EOMES+ progenitors on differentiation day 4, confirming that these perturbations affect self-renewal and cell fate decisions at the stage of neurogenesis instead of neuronal maturation (Extended Data Fig. 5d).

TF perturbations alter EN to IN output of individual RG

The influence of multiple TFs on the relative abundance of ENs and INs nominated RG lineage plasticity as a candidate developmental mechanism for observed composition differences. Composition differences could arise from perturbations of RG that preferentially affect EN-biased clones, mixed clones containing both ENs and INs, or IN-biased clones (Fig. 3b). To distinguish between these possibilities, we combined Perturb-seq with lineage tracing using STICR—a GFP-expressing lineage tracing library with static barcodes that can be measured by scRNA-seq1. We constructed an mCherry-expressing sub-library targeting NR2E1, ARX and ZNF219, and also included SOX2 and NEUROD2 as controls, the former known to impact both lineages and the latter known to promote the EN lineage (Fig. 3b).

As RG fate plasticity changes with maturation, we extended the experiments using this targeted dual library to primary human RG culture isolated from nine human individuals from GW16 to GW22, spanning the peak of excitatory and inhibitory neurogenesis12,27 to early stages of gliogenesis. Cells co-expressing mCherry and GFP were isolated on day 7 for scRNA-seq, and transcriptomes were reference mapped to the initial Perturb-seq data for cell-type annotation. This analysis yielded 129,003 cells with assigned sgRNAs and lineage barcodes (Fig. 3c and Extended Data Fig. 6a–j; Methods) and recapitulated the gene expression and cell-type composition consequences observed in the initial screen (Extended Data Fig. 6k,l).

We next investigated the landscape of clonal lineage relationships. We considered multicellular clones containing at least three cells with consistent sgRNA assignment, recovering 8,937 clones (Extended Data Fig. 7a; Methods). One person with low overall cellular coverage was dropped from this analysis. Remaining clones contained a mean of five cells in NT conditions (Extended Data Fig. 7b). The most abundant multi-class clones contained RG and ENs (35%) and RG and INs (20%) but, consistent with recent studies, we also observed around 15% clones containing both ENs and INs (Fig. 3d). Lineage coupling, defined as the normalized barcode covariance between cell types87,88, further supported the presence of mixed clones with a comparable linkage of dividing RG between both excitatory and inhibitory intermediate progenitor cells (IPCs) and immature neurons (Fig. 3e). Unsupervised clustering of clones based on cell-type composition using scLiTr[89](https://www.nature.com/articles/s41586-025-09997-7#ref-CR89 “Erickson, A. G. et al. Unbiased profiling of multipotency landscapes rev