3D Epigenome Mapping of Glial Cell Types in the Developing Human Cortex

Tissue dissociation, sample fixation and storage

Cells were isolated from the developing human cortex between gestational weeks (GW) 15 and 24 using a protocol adapted from a previously published method.³ The dissociated cells were washed twice with PBS and fixed in 2% paraformaldehyde (PFA) for 10 min at room temperature. Fixation was stopped by adding glycine to a final concentration of 200 mM and incubating the samples for an additional 5 min at room temperature. All subsequent steps were performed on ice or at 4 °C. Fixed cells were centrifuged at 1,000g, washed with PBS, passed through a 70-μm nylon mesh and washed again with PBS. The cells were then pelleted, snap-frozen and stored at −80 °C until further analysis.

Fluorescence-activated cell sorting

All fluorescence-activated cell sorting (FACS) procedures were performed on ice or at 4 °C. Approximately 1.5 × 10⁸ fixed cells were thawed and permeabilized in 1 ml of PBS containing 0.1% Triton X-100 for 15 min. Bovine serum albumin (BSA) was added to a final concentration of 1%, and cells were centrifuged at 1,000g for 8 min. The pellet was washed once with staining buffer (PBS containing 1% BSA) and resuspended in 100 µl of the same buffer. Cells were blocked with FcR Blocking Reagent (Miltenyi Biotech, 1:20) for 10 min and incubated with primary antibodies for 30 min.

FACS antibodies included PerCP-Cy5.5 anti-SOX2 (BD Biosciences, 561506, for radial glia (RG)), PE-Cy7 anti-EOMES (Invitrogen, 25-4877-42, for RG), unconjugated anti-HOPX (Proteintech, 11419-1-AP, for RG), Alexa Fluor 647 anti-OLIG2 (Abcam, ab225100, for oligodendrocyte precursor cells (OPCs) and microglia (MG)) and PE anti-PU.1 (Cell Signaling Technology, 81886, for OPCs and MG). Anti-HOPX was used at a 1:250 dilution; all other primary antibodies were used at 1:20. After staining, cells were washed twice with staining buffer, resuspended in 300 µl of buffer and incubated with Alexa Fluor 647 donkey anti-rabbit secondary antibody (Invitrogen, A-31573, for RG FACS only) at a 1:300 dilution.

Cells were sorted on BD FACSAria II instruments into collection buffer consisting of PBS with 5% BSA (Supplementary Figs. 1 and 2). Sorted cells were centrifuged at 1,000g for 10 min, snap-frozen on dry ice and stored at −80 °C before downstream processing. For FACS samples used in RNA sequencing, 1% RiboLock RNase Inhibitor (Thermo Scientific, EO0384) was added to all buffers.

RNA-sequencing library preparation and data analysis

Total RNA was extracted from sorted cell populations using the RNA FFPE kit (Qiagen, 73504), starting with 3 × 10⁵ to 1.8 × 10⁶ cells. RNA quality was assessed with an Agilent 2100 Bioanalyzer by calculating the percentage of RNA fragments longer than 200 bp (DV200). Samples with DV200 ≥ 30% were used for library preparation. Ribosomal RNA was depleted with the KAPA RNA HyperPrep Kit with RiboErase (HMR KK8560), followed by first- and second-strand cDNA synthesis, dA-tailing and sequencing-adapter ligation. Final adapters were added through PCR amplification. Libraries were paired-end sequenced on a NovaSeq S4 system using 100-bp paired-end reads.

Raw RNA-seq reads were trimmed to 100 bp with fastp⁷⁵ (v.0.22.0) and aligned to the hg38 human reference genome with STAR (v.2.7.10a) using standard ENCODE parameters. Strand-specific gene quantification was performed with RSEM (v.1.2.28) and GENCODE 38 annotations. Library quality was evaluated using the median TIN⁷⁶ score, which exceeded 60 for all samples. TMM-normalized reads per kilobase per million mapped reads (RPKM) were calculated with edgeR⁷⁷ (v.3.32.1). Mean expression values across biological replicates were used for downstream analyses.

Differentially expressed genes (DEGs) were identified with DESeq2⁷⁸ using a multifactor design that included genotype or individual to compare ventricular radial glia (vRG) and outer radial glia (oRG). Clustering was performed using regularized log-transformed data and hierarchical clustering based on sample distances. Batch effects associated with genotype or individual were removed with the removeBatchEffect function in limma⁷⁹ (v.3.46.0).

Characterizing bulk RNA sequencing with single-cell RNA sequencing

Cellular composition in bulk RNA-seq samples was estimated with CIBERSORTx⁸⁰ using a matched single-cell RNA-seq (scRNA-seq) reference dataset.⁹ Embryonic neurons (eN), inhibitory neurons (iN) and intermediate progenitor cell (IPC) subclusters were combined into a single cluster before the reference matrix was generated. The final reference included eN, iN, IPC, vRG, oRG, truncated RG, OPC and MG populations.

Raw counts from bulk vRG, oRG, OPC and MG libraries were imported from RSEM output using tximport(). CIBERSORTx was then used to estimate cell-type fractions with B-mode batch correction, relative run mode and 100 permutations.

ATAC-sequencing library preparation and analysis

ATAC sequencing was performed using a protocol adapted from a previously published method.³ Briefly, 50,000–100,000 formaldehyde-fixed and sorted cells were resuspended in nuclei extraction buffer containing 10 mM Tris-HCl pH 7.5, 10 mM NaCl, 3 mM MgCl₂, 0.1% Igepal CA630 and 1× protease inhibitor. Samples were incubated at 4 °C for 5 min. Cells were then resuspended in 50 μl of 1× TD buffer from the Nextera DNA Library Prep Kit (Illumina, FC-121–1030) and treated with 2.5 μl of TDE1 enzyme for 45 min at 37 °C with shaking at 4,500 rpm.

To reverse crosslinks, 150 μl of reverse-crosslinking solution was added. This solution contained 50 μl of 1 M Tris pH 8.0, 100 μl of 10% SDS, 2 μl of 0.5 M EDTA, 10 μl of 5 M NaCl, 800 μl of water and 2.5 μl of 20 mg ml⁻¹ proteinase K. Samples were incubated overnight at 65 °C. DNA was purified with the Qiagen MinElute kit (28004), PCR amplified and size-selected for fragments between 300 and 1,000 bp using AMPure beads (Beckman Coulter, wsr-450437). Libraries were sequenced on a NovaSeq S4 instrument using 100-bp paired-end reads.

ATAC-seq reads were trimmed to 100 bp with fastp⁷⁵ (v.0.22.0), aligned to hg38 and processed with the ENCODE ATAC-seq pipeline (https://github.com/kundajelab/atac_dnase_pipelines) using default settings. All replicates showed transcription start site enrichment scores above 7 (18.54–27.69) and fractions of reads in peaks above 0.3 (0.32–0.53). Optimal overlapping peaks for each cell type were used in downstream analyses.

Differentially accessible regions (DARs) between vRG and oRG were identified with DiffBind⁸¹ (v.3.4.11), using DESeq2 and the multifactor design described for DEG analysis. A consensus peak set was generated, followed by normalization and testing with a false discovery rate below 0.01. ATAC-seq clustering was based on Spearman correlations of reads within a merged peak set across all four cell types, calculated with multiBamSummary from deepTools⁸² (v.3.5.1).

Whole-genome bisulfite sequencing library preparation and analysis

DNA was isolated from sorted cell populations with the MagMAX FFPE DNA/RNA Ultra kit (Applied Biosystems, A31879), starting with 2.9 × 10⁵ to 1.5 × 10⁶ cells. Genomic DNA was sonicated to 300–600 bp using a Covaris M220, and fragment size was confirmed with an Agilent Bioanalyzer. Unmethylated lambda DNA was added at 0.5% of the total sample to monitor bisulfite conversion efficiency. Bisulfite conversion was performed using the EZ DNA Methylation Direct kit (Zymo, D5020). Whole-genome bisulfite sequencing (WGBS) libraries were prepared with the Accel-NGS Methyl-seq Combinatorial Dual Indexing kit (Swift Biosciences, 38096). Library size and adapter removal were verified using an Agilent Bioanalyzer. Libraries were sequenced on a NovaSeq instrument with 100-bp or 150-bp paired-end reads.

Paired-end WGBS FASTQ files were adapter-trimmed with Trim Galore v.0.6.6 (https://github.com/FelixKrueger/TrimGalore). An additional 15 bases were removed from the 5′ end of read 2 and the 3′ end of read 1. FastQC was used to evaluate raw and trimmed reads and to confirm removal of methylation bias introduced during end repair. Reads were aligned to hg38 with Bismark⁸³ v.0.16.2 and Bowtie2⁸⁴ v.2.4.1. Duplicate reads were removed with deduplicate_bismark before biological replicates from the same cell type were merged.

Lambda DNA reads were aligned separately to confirm bisulfite conversion efficiencies above 99% for all samples. Following replicate merging, methylation calls for each cytosine context were generated with bismark_methylation_extractor. Only cytosines covered by at least ten reads were retained for downstream analysis. Pairwise differentially methylated regions (DMRs) were identified with methylKit⁸⁵ and defined as regions showing at least a 25% methylation difference. MethylSeekR⁸⁶ was used to identify unmethylated regions, low-methylated regions (LMRs) and partially methylated domains.

PLAC-sequencing library preparation and analysis

PLAC sequencing was performed with the Arima-HiC+ Kit. Each library was prepared from 2–4 million cells fixed with 2% formaldehyde (F79-500). Following chromatin digestion and ligation, samples were sonicated with a Covaris S220 at 4 °C using a peak power of 105 W, duty factor of 5%, 200 cycles per burst and a treatment time of 300 s. Chromatin immunoprecipitation was performed with 2.5 μl of an H3K4me3 antibody (Millipore, 04-745). Sequencing adapters were added with the Swift Biosciences Accel-NGS 2S Plus DNA Library Kit, and libraries were amplified with KAPA HiFi HotStart ReadyMix. Libraries were paired-end sequenced on NovaSeq S4 instruments using 100-bp reads. Raw reads were trimmed to 100 bp with fastp⁷⁵ (v.0.22.0).

The MAPS⁸⁷ pipeline was used to identify significant H3K4me3-associated chromatin interactions at 2-kb resolution across interaction distances of up to 2 Mb. Raw reads were aligned to hg38 with BWA-MEM. Unmapped reads and read pairs with low mapping quality were removed, and the remaining pairs were processed as previously described. Read pairs were classified as AND, XOR or NOT interactions according to whether both, one or neither read overlapped a universal anchor.

To generate the universal anchor set, H3K4me3 peaks were called for each cell type with MACS2 using ‘-g hs –broad –nolambda –broad-cutoff 0.01’. Approximately 30 million read pairs with interaction distances shorter than 1 kb were used for each cell type. This produced 19,826 (37,241), 18,198 (35,783), 20,107 (38,480) and 18,730 (35,361) peaks (2-kb bins) in MG, OPC, oRG and vRG, respectively. Bedtools merge was then used to create a universal H3K4me3 anchor set containing 24,167 peaks (47,412 2-kb bins) (Extended Data Fig. 2b).

HPRep⁸⁸ confirmed the reproducibility of biological replicates. Approximately 10 million usable reads were used per sample to control for sequencing-depth differences, and all biological replicates showed Pearson correlations above 0.9. For final analyses, merged cell-type datasets were downsampled to approximately 60 million usable AND and XOR reads. Significant chromatin interactions were identified using a Poisson regression-based approach.

Defining candidate cis-regulatory elements

Candidate cis-regulatory elements (cCREs) were classified as general cCREs, accessible-region cCREs (cCREs^AR) or low-methylated-region cCREs (cCREs^LMR) using bedtools. Any base-pair overlap between cCREs^AR and cCREs^LMR was classified as a cCRE. In contrast, cCREs^AR and cCREs^LMR were defined as regions exclusively overlapping accessible regions or LMRs, respectively, using the ‘-v’ option to exclude the corresponding feature. Chromosomes X and Y were excluded. Intervene⁸⁹ was used to generate upset plots showing overlaps among cCREs and other genomic feature sets.

Transcription-factor motif enrichment analysis

Motif enrichment in cell-type-specific cCREs and DARs was assessed with HOMER³⁴ using findMotifsGenome.pl. Default settings were used, except that ‘-size given’ was specified. Significant motif-associated transcription factors were retained only when their expression exceeded 10 RPKM in the corresponding cell type.

Heatmap analysis of cell-type-specific cCREs

Cell-type-specific XOR interactions were first identified from datasets containing a cCRE within the distal interaction bin. For each bin pair, the normalized contact score was calculated as the observed contact count divided by the MAPS-generated expected count. Average ATAC-seq signal within each distal bin was calculated with bigWigAverageOverBed using counts-per-million (CPM)-normalized signal.

The heatmap displays the percentage contributed by each cell type relative to the summed signal across all cell types. When multiple peaks occurred within one distal bin, their average signals were summed. ATAC-seq counts were further corrected for sequencing depth by quantile normalization. Transcriptome values were calculated similarly using summed RPKM gene expression within the H3K4me3 anchor bin. Methylation values for distal bins were calculated as the mean CpG methylation of ATAC-seq peaks located within each bin. Finally, cell-type-specific interactions were retained only when they overlapped at least 50% of cCREs from the corresponding cell type.

Mouse enhancer transgenic assay

Candidate elements for VISTA mouse transgenic assays were selected from ATAC-seq peaks identified in RG, IPC, eN and iN populations.³ ATAC-seq reads were quantified within a merged peak set and quantile normalized. Peaks overlapping transcription start sites identified by cap analysis gene expression sequencing (CAGE-seq) were removed. Remaining regions were required to participate in a three-dimensional chromatin interaction and overlap the top 15,000 accessible orthologous regions in the mouse embryonic brain at E11.5 (refs. ^90,91).

Candidate peaks were further filtered for strong ATAC-seq signal (>0.8 CPM) and cell-type-specific enrichment (>0.5 CPM difference). These criteria yielded 61 candidate regions, from which 29 elements were selected manually. Most overlapped cCRE annotations (20 of 29) or chromatin-interacting regions (20 of 29; 16 were cCREs) identified in the vRG, oRG, OPC and MG datasets generated in this study (Supplementary Table 4).

Transgenic mouse embryos were generated as previously described,⁹² except that embryos were collected at embryonic day 12.5 (E12.5). Assays were performed in Mus musculus FVB mice. Animals were maintained on a 12-h light/12-h dark cycle, with lights on from 06:00 to 18:00, at 20.6–23.9 °C (69–75 °F) and 30–70% humidity.

VISTA element enrichment analysis

To identify cCREs associated with functional VISTA elements,^24,25 we analyzed neural VISTA elements (n = 1,233) and negative controls (n = 1,868). Neural elements were defined by activity in the neural tube, forebrain, midbrain, hindbrain, dorsal root ganglion, cranial nerve or trigeminal nerve. Negative controls showed no detectable signal. These elements were overlapped with cCREs involved in three-dimensional chromatin interactions and compared with distance-matched control elements generated as previously described.⁸⁷ The proportions of neural and negative elements overlapping cCREs or control regions were compared using Fisher’s exact test.

Potential target genes for VISTA elements were identified by linking positive neural VISTA elements that were both accessible and classified as LMRs to genes with bedtools. Candidate targets were retained when expression exceeded 1 RPKM in the matching cell type.

TOBIAS footprinting and transcription-factor network analysis

Footprinting analysis was performed using merged BAM files from multiple sequencing runs for each cell type. Transcription-factor motifs were obtained from the HOCOMOCO v.11 database.⁹³ TOBIAS³⁰ was run with its standard workflow and parameters to identify transcription-factor footprints in each cell type. ENCODE blacklist regions were supplied through the TOBIAS ATACorrect tool to correct for Tn5 insertion bias.

Differential binding was evaluated pairwise. Only transcription factors expressed at more than 10 RPKM were included in volcano plots and Pearson correlation analyses with gene expression. Motif-binding predictions were considered cell-type-specific when a motif was predicted to be bound in the corresponding cell type and the absolute log₂[fold change] between cell types exceeded one. Target genes were identified as described above. For Cytoscape network visualization, only differentially expressed genes were displayed.

Isolation and in vitro culture of outer radial glia

The ventricular zone, inner subventricular zone and outer subventricular zone were dissected from a primary human cortical tissue sample collected at GW20 and dissociated with the Papain Dissociation System (Worthington Biochemical). Cells were transduced with lentiviruses expressing GFP and either scrambled control shRNA (shCTRL) or shRNAs targeting LHX2 (shLHX2_1 and shLHX2_2) (Supplementary Table 7).

After 72 h, cells were blocked with FcR Blocking Reagent (Miltenyi Biotech, 1:20) for 10 min and incubated with antibodies for 30 min. FACS antibodies included APC-conjugated LIFR (leukaemia inhibitory factor receptor; R&D Systems, FAB249A) and PE-Cy7 anti-ITGA2 (BioLegend, 359314). GFP⁺ITGA2⁺LIFR⁺ oRG cells were collected. For each condition, 50,000 cells were seeded into four wells of a 24-well plate and cultured for 7 days in RG differentiation medium containing DMEM/F12, 2 mM GlutaMAX, 2% B27 without vitamin A, 1% N2 and 1× penicillin–streptomycin.

Samples were processed with the 10× GEM-X Universal 3′ 4-plex on-chip multiplexing assay, targeting 1,300–5,000 cells per replicate. Individual libraries were pooled and sequenced on an Illumina NovaSeq X Plus instrument.

Single-cell RNA-sequencing analysis of outer radial glia

The Cell Ranger (v.9.0.1) multi pipeline was used for cell-barcode calling, read alignment and quality assessment. A custom human reference genome based on GRCh38 and GENCODE v.32/Ensembl98 was used, with enhanced GFP (eGFP) added according to 10× Genomics protocols. Ambient RNA was removed from pooled gel bead-in-emulsion libraries with CellBender⁹⁴ (v.0.3.2) using the remove-background command.

High-quality cells were retained when they contained more than 1,000 detected genes, fewer than 5% mitochondrial reads and a scDblFinder⁹⁵ (v.1.23.4) doublet score below 0.3. Data-quality metrics are summarized in Supplementary Table 7. Log normalization with a size factor of 10,000, data scaling and regression of G2/M and S-phase cell-cycle markers were performed in Seurat⁹⁶ (v.5.4.0).

A nearest-neighbour graph was constructed using the first 30 principal components, and clusters were identified with the Louvain algorithm. Clusters with low unique molecular identifier counts, which were likely to represent low-quality cells, were removed before repeating the clustering analysis. Cell types were assigned according to established marker genes (Extended Data Fig. 8c).

scDist³⁵ (v.1.1.5) was used to determine which cell type showed the greatest transcriptional difference between LHX2 knockdown and control conditions. SCTransform⁹⁷ (v.0.4.3) was used for normalization and scaling in this analysis. Changes in cell-type composition were assessed with scCODA⁹⁸ (v.0.1.9) using default settings.

Linkage disequilibrium score regression

Linkage disequilibrium score regression (LDSC) was performed for each complex neuropsychiatric disorder using joint models. Models incorporated either cCREs involved in H3K4me3-mediated interactions or distal 2,000-bp bins targeting anchor bins that overlapped H3K4me3 signal across all cell types, together with a baseline model.⁹⁹ These approaches were used for Fig. 4a and Extended Data Fig. 8c, respectively. For Extended Data Fig. 8a,b, joint models included the baseline model, data generated in this study and ATAC-seq peaks from RG, IPC, eN and iN datasets, classified according to whether they participated in three-dimensional chromatin interactions.

Training the GKM accessibility prediction model

Following previously published approaches,^69,70 a GKM-SVM classifier was trained to predict chromatin accessibility. For each cell type, the top 75,000 MACS2 peaks ranked by peak score were selected. Peaks containing ambiguous N bases were excluded, and all remaining peaks were standardized to 1,000 bp by extending 500 bp from each summit.

Negative training sequences were generated with genNullSeqs to match repeat content and GC composition. For model evaluation, all chromosomes except chromosome 2 were used for training, while chromosome 2 served as the independent test set. The gkmtrain function from the LS-GKM package was run with default parameters and the wgkm kernel (t = 4). Performance was evaluated on the test set using the PRROC R package. After successful validation, the classifier was retrained using the complete dataset. For deltaSVM analysis, all 11-mers were generated with nrkmers.py and scored with gkmpredict.

In silico analysis of genetic variants

Reference and alternative alleles were evaluated in silico using methods adapted from previous work.⁶⁹ Alzheimer’s disease,⁴⁷ rare non-coding⁵⁰ and schizophrenia (SCZ)⁵³ variants were filtered with bedtools to retain variants located in accessible regions in at least one cell type. Each variant was extended by 100 bp upstream and downstream to generate a 200-bp sequence.

Three methods—deltaSVM,⁴⁵ in silico mutagenesis (ISM) and GkmExplain⁴⁶—were used to identify candidate regulatory SNPs. For GkmExplain, only the difference in the central 50-bp region was considered. The methods showed highly consistent results, with Pearson correlation coefficients above 0.99, although some outliers were observed (Extended Data Figs. 8a and 9a).

SNPs were considered significant when GkmExplain, ISM and deltaSVM scores fell outside the 95% confidence interval of null t-distributions. Prominence and magnitude scores derived from seqlets matching potential transcription-factor motifs were also used to support variant prioritization, as previously described.⁶⁹ Scores were calculated for both positive and negative contributions because the model may learn motifs associated with transcriptional repression and reduced chromatin accessibility. These values are provided in Supplementary Tables 9 and 10.

Human accelerated region enrichment testing

Enrichment of cCRE groups within 3,168 annotated human accelerated regions (HARs)⁵⁸ was evaluated by comparing the observed overlap with an empirical null distribution. For example, 4,515 unique vRG cCREs overlapped 30 HARs. A total of 1,000 random samples of 4,515 cCREs were drawn from the 121,317 cCREs obtained by merging vRG, oRG, OPC and MG cCREs, and the number of overlapping HARs was recorded for each sample.

The resulting empirical null distribution was used to calculate a z score by comparing the observed overlap of 30 HARs with the randomized overlaps. For Extended Data Fig. 9, 241,128 cCREs representing the union of vRG and oRG cCREs and IPC, iN and eN ATAC-seq peaks were sampled.

Identification of variants in human accelerated regions

HAR variants were obtained using a method adapted from previous work.⁶⁰ Alignments for HARs accessible in oRG cells were retrieved with mafsInRegion. The msa_view command was used to convert alignments into multiple-sequence alignment format, retaining only hg38 and pantro4 sequences. The resulting alignment was converted to VCF format with snp-sites. Each variant was extended by 100 bp upstream and downstream to create a 200-bp sequence, which was evaluated with deltaSVM, ISM and GkmExplain.

Prediction of transcription-factor binding changes

Potential transcription-factor binding changes between human and chimpanzee sequences were predicted with motifbreakR.⁶¹ Motifs were obtained from HOCOMOCO v.11 core datasets A, B and C. Default parameters were used, including a threshold of 1 × 10⁻⁴ and the ‘ic’ method. Only strong predicted effects were retained, and motifs corresponding to transcription factors expressed above 5 RPKM were included in downstream analyses.

CRISPR interference knockdown of HARsv2_1313

The CROP-seq-opti-eGFP vector was developed from the CROP-seq-opti backbone (Addgene, 106280) to co-express puromycin resistance (PuroR), eGFP and dual guide RNAs (gRNAs). eGFP was inserted downstream of the PuroR coding sequence and connected by a P2A self-cleaving peptide. Two gRNA pairs targeting HARsv2_1313 and one pair targeting the ROCK2 promoter were designed with CHOPCHOP¹⁰⁰ (Supplementary Table 10). Dual gRNAs were cloned into the CROP-seq-opti-eGFP vector using a previously described protocol.¹⁰¹

For lentivirus production, 7.5 μg of each gRNA plasmid per T-75 flask was cotransfected with pMD2.G (1.5 μg; Addgene, 12259) and psPAX2 (4.5 μg; Addgene, 12260) into 293T-LentiX cells (Takara Bio, 632180) using PolyJet transfection reagent (SignaGen, SL100688). Culture medium was replaced 18 h after transfection, and viral supernatants were collected once daily for 3 consecutive days. Lentivirus was passed through a 0.45-μm syringe filter and concentrated using Amicon Ultra-15 centrifugal filters (Millipore, UFC901024).

Differentiation of induced pluripotent stem cells into neural progenitor cells

Human WTC11 induced pluripotent stem (iPS) cells stably expressing dCas9-KRAB⁴⁷ were maintained in mTeSR medium (STEMCELL Technologies, 100-0274) and routinely tested for mycoplasma. Cells were dissociated with Accutase and seeded on Matrigel-coated plates at 2.5 × 10⁵ cells per cm² in basal medium containing DMEM/F12, 1× N2, 1× B27 without vitamin A, 100 μM non-essential amino acids, 0.5 mg ml⁻¹ BSA, 1× penicillin–streptomycin and 100 μM 2-mercaptoethanol. The medium was supplemented with 20 ng ml⁻¹ FGF2 and 10 μM Y-27632.

At approximately 90% confluency, designated day 0, the medium was replaced with neural induction medium consisting of basal medium supplemented with 10 μM SB431542, 100 nM LDN193189 and 1 μM cyclopamine. Fresh neural induction medium was added daily from day 1 through day 11. On day 12, neural progenitor cells (NPCs) were dissociated with Accutase and replated on Matrigel-coated plates at 4.5 × 10⁵ cells per cm² in neural stem-cell medium (NSCM; Thermo Fisher, A10509-01) supplemented with 20 ng ml⁻¹ FGF2, 20 ng ml⁻¹ epidermal growth factor, 1× GlutaMAX, 30 μg ml⁻¹ heparin, 0.2 mM ascorbic acid and 1× penicillin–streptomycin. NSCM also contained 10 μM Y-27632 during the first passage.

From day 13 onward, NSCM was replaced every other day. Cultures were passaged approximately every 5 days after reaching superconfluency. NPCs were maintained in NSCM and used for downstream experiments between days 26 and 30.

Reverse-transcription quantitative PCR

NPCs were transduced with lentivirus on day 21. Three days later, infected cells were enriched by puromycin selection for 3 days and allowed to recover for a further 3 days. Total RNA was extracted with the AllPrep DNA/RNA Mini Kit (Qiagen, 80204). For cDNA synthesis, 200 ng of RNA was reverse transcribed using the iScript cDNA Synthesis Kit (Bio-Rad, 1708891). Quantitative PCR was performed with NEBNext Ultra II Q5 Master Mix (NEB, M0544) containing 1× SYBR Green. ROCK2 expression was normalized to GAPDH.

Immunocytochemistry and image analysis

NPCs were transduced with lentivirus on day 21, replated in 96-well plates on day 27 and fixed with 4% paraformaldehyde on day 29. Cells were blocked in PBS containing 0.1% Triton X-100 and 5% horse serum. Primary antibodies were diluted in blocking solution and applied overnight at 4 °C, followed by secondary-antibody incubation for 1 h at room temperature. Antibodies included rat anti-Ki67-Alexa Fluor 647 (BioLegend, 151206; 1:200), rabbit anti-EMX2 (GeneTex, GTX17164; 1:400) and donkey anti-rabbit-Alexa Fluor 546 (Invitrogen, A10040; 1:500).

Images were acquired with an Opera Phenix Plus high-content imaging system (Revvity) using a ×20 water-immersion objective and analyzed with Harmony software. Nuclei were identified using a basal Ki67 signal-intensity threshold of 0.30. Condensed Ki67 puncta were identified as high-intensity spots with relative spot intensity above 0.105. A cell was classified as Ki67-positive (Ki67⁺) when at least one condensed Ki67 punctum was detected within its nucleus. Mean GFP intensity was measured for every cell.

All downstream analyses were performed in R. GFP intensity was log₁₀-transformed for thresholding. Cells with log₁₀[GFP intensity] above 2.4 were considered GFP-positive (GFP⁺) and included in subsequent analyses. This threshold was selected from the GFP-intensity distribution of non-transduced controls to distinguish background from true expression.

GFP⁺ cells from all treatment groups were pooled and ranked by raw GFP intensity. Cells were divided into four quartiles: 0–25%, 25–50%, 50–75% and 75–100%. Quartile assignments were then returned to the individual treatment groups. For each biological replicate and condition, the percentage of Ki67⁺ cells was calculated within each GFP quartile. Ki67⁺ percentages in the two highest quartiles, where the CRISPRi knockdown effect reached saturation, were compared between conditions. Bar plots show mean ± s.e.m., and P < 0.05 was considered statistically significant.

To test whether Ki67 positivity changed across GFP quartiles, simple linear regression was performed within each treatment condition, with quartile number coded from 1 to 4 as the independent variable and the percentage of Ki67⁺ cells as the dependent variable. The slope and associated P value were used to assess trend direction and significance. Regression lines are shown with 95% confidence intervals. Differences in trends between control and treatment groups were tested using models containing an interaction between quartile and treatment. The interaction-term significance indicated whether the slope differed between conditions.

Luciferase reporter assays

The Dual-Luciferase Reporter Assay System (Promega, E1910) was used to compare the activity of chimpanzee and human HAR variants. HAR sequences were amplified from human WTC11 iPS-cell genomic DNA or chimpanzee C3649 genomic DNA (ref. ¹⁰²) using NEBNext Ultra II Q5 polymerase (NEB, M0544L) and identical primers for both genomes (Supplementary Table 11). The pGL4.13 vector (Promega, E6681) was digested with Xho I and Nco I. HAR elements and a synthesized minimal promoter were inserted by Gibson assembly (NEB, E2621L), and constructs were validated by Sanger sequencing.

The ventricular zone, inner subventricular zone and outer subventricular zone were dissected from primary human cortical samples collected between GW17 and GW20 and dissociated with the Papain Dissociation System (Worthington Biochemical). Cells were plated in poly-d-lysine-coated 24-well plates at 1.5 × 10⁶ cells per well. Culture medium consisted of DMEM/F12 with GlutaMAX, 1× B27 without vitamin A, 1× N2, 0.1 mM 2-mercaptoethanol, 1× non-essential amino acids, 20 ng ml⁻¹ FGF2, 20 ng ml⁻¹ brain-derived neurotrophic factor, 20 ng ml⁻¹ pleiotrophin, 20 ng ml⁻¹ platelet-derived growth factor DD and 1× penicillin–streptomycin.

After 24 h, cells were cotransfected in triplicate with approximately 495 ng of HAR reporter vector and pRL-CMV-Renilla luciferase vector (Promega, E2261) at a 10:1 ratio using Lipofectamine 3000 (L3000001). After 48 h, cells were lysed with passive lysis buffer. Each well was washed with 500 μl PBS, treated with 100 μl passive lysis buffer and rotated at room temperature for 15 min. Luciferase activity was measured with a GloMax plate reader. Background activity from non-transfected controls was subtracted, and firefly luciferase activity was normalized to the mean minimal-promoter signal for each sample.

Ethics statement

Human prenatal tissue samples were collected as de-identified specimens following informed consent and in accordance with all applicable legal, institutional and ethical requirements. Sample acquisition, collection and use were approved by the Human Gamete, Embryo and Stem Cell Research Committee and Institutional Review Board at the University of California, San Francisco, USA. The Human Gamete, Embryo and Stem Cell Research Committee protocol number was 10-05113. All experiments followed the approved protocol.

Animal procedures were reviewed and approved by the Lawrence Berkeley National Laboratory Animal Welfare and Research Committee. Animals were inspected weekly by the committee chair and the head of the animal facility in consultation with veterinary staff. The LBNL Animal Care Facility is accredited by the American Association for the Accreditation of Laboratory Animal Care International.

Reporting summary

Additional information about the research design is available in the Nature Portfolio Reporting Summary linked to this article.