Gene switches far from the gene control brain support cells linked to Alzheimer's disease

To identify functional enhancers in human primary astrocytes, we used normal human astrocytes (NHAs) from fetal brain (Methods). The transcriptomes of NHAs with or without dCas9–KRAB integration (Methods) were highly correlated (r = 0.95, ρ = 0.9; Extended Data Fig. 1), and clustered with fetal astrocytes (ref. ¹⁶; Fig.1b). We further confirmed that the developmental stage of NHAs was preserved at transcriptome and chromatin state level—99.7% of dCas9–KRAB–NHA cells were annotated as fetal astrocytes based on single-cell RNA-seq (scRNA-seq) with a human brain reference dataset¹⁷ (Fig.1c), and ATAC–seq data from dCas9–KRAB–NHAs clustered with snATAC–seq from fetal astrocytes¹⁷ (Extended Data Fig. 1). For the rest of the article, we use ‘NHAs’ to refer to dCas9–KRAB–NHAs.

Candidate regions were selected from intergenic PsychENCODE¹³ candidate enhancers (>2 kb from any gene) and required to be in an open chromatin state in NHAs and in at least two other astrocytes or NeuN⁻ datasets^18–21. In addition, we required candidate enhancers to reside in a topologically associated domain (TAD)²² containing at least one gene highly expressed in NHAs (Methods). This yielded 979 intergenic candidate enhancers (Fig. 1a, Supplementary Table 1 and Extended Data Fig. 2), none of which overlapped with copy number variations identified in NHAs (Supplementary Table 1).

The candidate regions were enriched for distal enhancers²³ (odds ratio (OR) = 43.8, P = 4 × 10⁻²⁵⁹, Fisher’s exact test), depleted of promoters (OR = 0.02, P = 2 × 10⁻²⁷, Fisher’s exact test) and proximal enhancers (OR = 0.05, P = 4 × 10⁻⁴⁴, Fisher’s exact test; Fig.1d). Among the four ATAC–seq datasets used for candidate selection^18–21, the chromatin state of candidate regions best resembled that of brain organoid astrocytes (93% overlap; Fig.1e) consistent with NHA’s fetal origin. Most candidate regions (61%, n = 602) were accessible in both fetal-derived and adult-derived astrocytes (Fig.1e).

To further characterize the regulatory landscape captured by candidate enhancers, we constructed a co-accessibility network using snATAC–seq from the human brain across developmental stages¹⁷ (Methods). Candidate enhancers clustered into six modules, and 86% of enhancers belonged to modules enriched in astrocytes or more broadly glial cells (Fig.1f,g and Supplementary Table 1). Five modules showed developmental stage-specific chromatin state (analysis of variance P < 0.05; Fig.1f,g), and one captured developmental variation within astrocytes (P = 0.03; Fig.1f,g). Bulk brain-derived ATAC–seq and H3K27ac from ref. ²⁴ revealed cell-type-specific profiles consistent with those observed in snATAC–seq data, with module M3 showing the strongest astrocyte-specific signal, while H3K4me3 showed low coverage in astrocytes, as expected (Extended Data Fig. 3).

Overall, our results demonstrate that NHAs retain a fetal astrocyte identity, and that the candidate enhancer set captures astrocyte-specific regions, as well as regions shared across brain cell types and developmental stages. Although our selection prioritized broadly relevant candidate enhancers, that is, in open state across fetal and adult astrocytes, the candidate set also includes regions with developmentally regulated chromatin states.

We constructed an sgRNA library using the CROP-seq-opti-GFP backbone containing 2,478 sgRNAs targeting the 979 candidate enhancers, 20 positive and 50 negative control sgRNAs (Supplementary Table 2), as well as separate positive and negative control libraries (n = 250 sgRNAs each; Methods). Candidate enhancers varied in length between 87 and 1,136 bp with a median of 298 bp (Extended Data Fig. 3), and were tiled with 1 sgRNA per 200 bp and ≥2 sgRNAs per enhancer (Extended Data Fig. 3).

sgRNAs were transduced into NHAs at a low multiplicity of infection¹ (MOI) of 5, to preserve viability. scRNA-seq (Fig.1a,c) was performed 7 days post-transduction (12.5B reads, 47,577 high-quality cells, median depth of 52,472 unique molecular identifiers (UMIs) per cell; Methods). Each enhancer was targeted in an average of 327 cells (range = 65–1,301), and a median of 7.4 guides was detected per cell (Extended Data Fig. 4).

For each enhancer, differential expression (DE) analyses compared the targeted and nontargeted cells, using SCEPTRE²⁵, assessing all expressed genes within 500 kb from the enhancer midpoint (Methods; Extended Data Fig. 4). Significant effects were defined using an empirical false discovery rate (FDR) of <0.1 as previously described¹ (Methods). Robust dCas9–KRAB silencing activity was confirmed by significant downregulation of 122 of 125 positive control genes (Supplementary Table 3), while a well-calibrated statistical approach without P-value inflation was confirmed based on the negative control library (Extended Data Fig. 4).

To identify well-powered nonfunctional candidate regions (that is, regions showing no effect in the CRISPRi screen despite sufficient statistical power), we used a simulation-based method¹⁴, and identified 605 regions with >80% power to detect a 15% expression change for at least one gene (Methods; Supplementary Table 3). This set was used for downstream comparisons of functional and nonfunctional candidate regions.

Functional enhancers were enriched for prior CRISPRi/MPRA hits in K562↗ (ref. ⁶⁷) (P = 2.2 × 10⁻³, OR = 2.71, Fisher’s exact test), and for predicted NHA superenhancers²⁶ (P = 1.0 × 10⁻⁶, OR = 4.0, Fisher’s exact test; Fig.2c). Functional enhancer–gene pairs (EGPs) were depleted of interactions crossing TAD boundaries (OR = 0.28, P = 4.0 × 10⁻⁶, Fisher’s exact test; Fig.2c). Furthermore, enhancer-regulated genes significantly overlapped genes previously predicted to contribute to the core regulatory circuitry in astrocytes²⁷ (P = 2.7 × 10⁻⁶, OR = 4.40, Fisher’s exact test).

To validate the results with an independent method, we individually silenced 19 functional enhancers (corresponding to 20 EGPs), prioritizing those that controlled genes dysregulated in AD (Fig. 5b) and spanned a range of effect sizes. Eighteen of the twenty EGPs (90%) showed significant effects (FDR < 0.05) on the target gene expression measured with NanoString nCounter (Methods; Supplementary Table 3). All 20 EGPs showed effects in the same direction as in the single-cell screen, with high correlation between effect sizes measured with the two methods (r = 0.86; Fig.2d). Seven EGPs were further tested through qRT–PCR in an independent NHA line of a different sex (Methods), confirming six of seven (Extended Data Fig. 4). To assess the effect of enhancer activity on protein expression, we performed a western blot for FTH1 after individually silencing Enh210, which regulates FTH1 (Methods). This led to a 24% decrease in FTH1 protein level (Extended Data Fig. 4).

These findings demonstrated that the identified enhancer–gene interactions aligned with known enhancer properties and were well-validated across methods and cell lines.

Although enhancers were significantly more likely to regulate the nearest gene than expected by chance (P < 2.2 × 10⁻¹⁶, OR = 11.7, Fisher’s exact test), only 38.6% EGPs involved the nearest gene. An additional 15.8% involved in a more distal gene with no intervening gene, 29.1% skipped a highly expressed gene to regulate a more distal gene and 16.4% skipped only lowly or nonexpressed genes (Fig.3b). When not targeting their nearest gene, enhancers bypassed up to 34 genes (mean = 7.4).

To assess enhancer RNA transcription in astrocytes, we applied transient-transcriptome sequencing (TT-seq²⁹) to NHA cells (Extended Data Fig. 5 and Supplementary Table 1; Methods). Functional enhancers were strongly associated with transcription (P = 1.3 × 10⁻¹¹, OR = 3.7, Fisher’s exact test), with 73.1% being transcribed (Fig. 3c,d). Among transcribed enhancers, functional enhancers were more likely to be bidirectional (P = 0.01, OR = 1.8, Fisher’s exact test; Fig. 3d), and had higher transcription levels (P = 0.045, Wilcoxon test; Fig. 3e).

Using snATAC–seq data across cell types and developmental stages (Fig.1f,g), we found that functional enhancers showed higher module membership (kME values) than inactive candidates for the astrocyte-enriched module M2 (P = 0.04, Wilcoxon test). Moreover, 48% of functional enhancers (n = 70), considered individually, showed significant variation in chromatin state across cell types, 4% (n = 7) across developmental stages and 34% (n = 50) across both cell type and developmental stages (Fig. 3f,g; analysis of variance adjusted P < 0.05).

When comparing H3K27ac coverage in astrocytes²⁴ between functional and nonfunctional candidate regions, there was no significant difference (P = 0.21, t test), indicating that, among open chromatin regions, higher H3K27ac coverage values are not predictive of enhancer function. Furthermore, when considering functional enhancers only, there was no significant difference in H3K27ac between astrocytes and other cell types (P = 0.66, t test). However, H3K27ac at functional enhancers belonging to the astrocyte-specific co-accessibility module M3 was significantly enriched in astrocytes compared to other cell types (P = 3.6 × 10⁻⁷, t test). The promoter mark H3K4me3 (ref. ²⁴) at functional enhancers was significantly depleted in astrocytes compared to other cell types, as expected (P = 0.005, t test).

Together, these data add several critical layers to the AstroREG resource, including eRNA expression data and characterization of chromatin state across cell types and developmental stages, and underscore the limitation of noncoding variant annotation to the nearest gene.

To investigate astrocyte-specific regulatory circuits, we performed TF footprinting using TOBIAS³⁷ on the deeply sequenced NHA ATAC–seq data (>800 million reads). Footprinting identifies TF motifs overlapping ATAC–seq footprints, regions of depleted signal within open chromatin indicative of protein binding. Motifs are classified as ‘bound’ or ‘unbound’ based on the footprint overlap, although overlapping TF motifs cannot be resolved to a specific TF. Despite this limitation, ATAC–seq footprinting enables genome-wide detection of candidate TF binding events. All ‘bound’ TF calls for screened regions are listed in Supplementary Table 5.

Functional enhancers were significantly more likely to be bound by TFs expressed in astrocytes and showed a stronger footprinting signal (Fig. 4c). A total of 83 TFs showed enriched binding to functional enhancers compared to nonfunctional enhancers (FDR < 0.05, Fisher’s exact test; Extended Data Fig. 7 and Supplementary Table 5), including 13 with particularly strong enrichment (OR > 3), of which FOXK1 is involved in astrocyte activation³⁴ while RUNX2 mediates astrocyte maturation³⁸.

We then constructed a distal regulatory network by integrating TF footprinting, CRISPRi screening and human brain snRNA-seq and snATAC–seq¹⁷. We defined astrocyte-specific enhancers using snATAC–seq, and astrocyte-specific genes and TFs using snRNA-seq from the same samples¹⁷ (Methods). Functional EGPs were more likely to be coregulated across cell types (P = 0.01, OR = 1.64, Fisher’s exact test), and called astrocyte-specific (P = 4.8 × 10⁻⁹, OR = 3.87, Fisher’s exact test). The regulatory network consisted of (1) 18 astrocyte-specific genes regulated by 19 astrocyte-specific enhancers (Fig.4d,e) and (2) 36 astrocyte-specific TFs belonging to 20 motif clusters from the JASPAR database³⁹, which bind at least one of the 19 astrocyte-specific enhancers (Methods; Fig.4e, Extended Data Fig. 7 and Supplementary Table 5). Our data allowed to construct network edges that are experimentally supported by CRISPRi screen data (enhancer–gene edges), and based on ATAC–seq footprinting for TF–enhancer edges (Fig.4e).

To determine TF–enhancer interactions supported by chromatin immunoprecipitation followed by sequencing (ChIP–seq) data from glial cells, we integrated ChIP–seq data from the following two complementary resources: ref. ⁴⁰ and the ChIP–Atlas⁴¹ (Methods). The data in ref. ⁴⁰ include TF–ChIP–seq from bulk brain tissue and sorted nuclei (neurons (NeuN⁺), oligodendrocytes (Olig2⁺) and astrocytes/microglia (NeuN⁻/Olig2⁻)) for the following six TFs in our network: FOS, ZBTB7B, MEIS2, NFIB, REST and SOX9. The ChIP–Atlas provides additional glia-relevant data, including TP53 in iPSC-derived astrocytes, TEAD1 in glioblastoma (GBM) cells and JUN in RPE-1 cells (a glial-like retinal cell type of neuro-ectodermal origin).

Across 33 predicted TF–enhancer pairs for which ChIP–seq data were available, 24 (72%) were supported by binding of the relevant TF to the predicted enhancer in at least one dataset (Fig. 4e, Extended Data Fig. 7 and Supplementary Table 5). Notably, all TF–enhancer overlaps identified in ref. ⁴⁰ were detected in the astrocyte-enriched fraction (NeuN⁻/Olig2⁻) or bulk tissue, with none observed exclusively in neurons or oligodendrocytes—providing further support for the astrocyte-specific nature of this regulatory network.

The AP1 complex components FOS and JUN are known primarily as immediate early genes, upregulated in response to neuronal activity and have a role in cellular proliferation in astrocytes⁴². Our data demonstrate that, in a resting state, astrocytes show the highest expression of JUN and FOS among brain cell types, and these TFs bind to enhancers that regulate astrocyte-specific genes (Fig.4e, Extended Data Fig. 7 and Supplementary Table 5). Notably, 11 of 12 FOS/JUN–enhancer interactions detected by ATAC–seq footprinting were confirmed by TF–ChIP, with 7 in both datasets (ref. ⁴⁰ and ChIP–Atlas).

Overall, these results uncover a regulatory network involved in astrocyte-specific gene expression, highlighting key transcription factors (TFs) and distal enhancers that contribute to the regulatory landscape in astrocytes.

Several enhancer-regulated genes were differentially expressed in AD cases versus controls, and had expression levels associated with neuropathological features and cognitive function (Fig. 5b). Among these, LGALS3 is an astrocyte-enriched gene (Fig. 4b,d) encoding Galectin-3, which modulates microglial activation during neurodegeneration⁴⁸. LGALS3 also mediates the astrocyte-microglia crosstalk involved in the neuropathological effects of APOE4, the strongest genetic risk factor for AD⁴⁹, and was downregulated in the brains of female AD donors carrying the APOE4 allele⁴⁹. We identified three enhancers controlling the expression of LGALS3, two showing astrocyte-specific open chromatin state, thereby serving as potential targets for cell-type-specific manipulation of LGALS3 expression in the brain. In addition to LGALS3, we uncovered functional enhancers that show cell-type-specific chromatin state for 11 other genes dysregulated in AD astrocytes (Fig. 5b).

Individual enhancer silencing experiments with NanoString readout (summarized in Fig. 2d) confirmed the regulatory effects of 16 enhancers on genes dysregulated in AD, including all three enhancers controlling LGALS3 expression (Fig. 5c).

Collectively, these data reveal that enhancer–gene regulatory interactions identified through CRISPRi screening are relevant to diseases, particularly in the context of AD.

Using all those significant eQTLs reported by each dataset, only 27 of the 158 regulatory interactions (17%) were captured by any of the eQTL datasets (Fig.6b and Extended Data Fig. 8). In addition, 43% had no significant eQTL reported, and, in 40% of cases, enhancers were linked to a different gene by eQTL than by CRISPRi screening in astrocytes. Using the more stringent fine-mapped eQTLs, 20 of the 158 regulatory interactions were captured (13%), 57% had no eQTL assignment and 30% were assigned to a different gene. Notably, eQTLs linking the enhancer to the same gene as identified in the CRISPRi screen data were more likely to be replicated across datasets (OR = 5.1, P = 2 × 10⁻⁴, Fisher’s exact test; Fig.6c). Moreover, astrocyte-specific eQTL calls obtained by deconvolution of bulk RNA-seq⁵² or scRNA-seq⁵³ were markedly underpowered, capturing none of the 158 regulatory interactions (Extended Data Fig. 8).

For example, we identified the following two enhancers regulating ID3 (inhibitor of DNA binding, also known as inhibitor of differentiation): Enh53 (chr1: 23,554,524–23,555,279, 4.6-kb downstream) and Enh54 (chr1: 23,562,311–23,562,640, 3.0-kb upstream; Extended Data Fig. 8). ‘ID3’ encodes a transcriptional regulator highly expressed in neural progenitors, where it controls proliferation⁵⁴, while in differentiated cells it is most highly expressed in astrocytes (Fig.4d). Enh53 and Enh54 overlap superenhancers⁵⁵, were previously shown to control ‘ID3’ expression in human primary macrophages⁵⁵ and we confirmed their regulatory effects in astrocytes by NanoString and qRT–PCR in an independent cell line (Fig. 5c and Extended Data Fig. 4). However, no brain eQTL dataset captured these interactions, possibly due to a combination of low ID3 expression, low astrocyte abundance and/or the lack of common variants within enhancers.

Together, these data demonstrate the power of CRISPRi screening of enhancers to provide a new functional annotation for noncoding regions in brain cell types.

To progress from region-level to nucleotide-level annotation, we performed in silico mutagenesis using a deep-learning framework (Beluga) in ref. ⁵⁶, predicting the effect of sequence changes on chromatin state in astrocytes (Methods). We also obtained disease-impact scores (DISs)⁵⁶ for all SNPs located within 1 kb of the tested enhancers, which estimates the probability of a variant being pathogenic⁵⁶, and found that SNPs within 1 kb from functional enhancers were more likely to have a significant DIS than those in nonfunctional candidates (P = 0.009, OR = 1.16, Fisher’s exact test).

Using WGS data from the AD sequencing project (ADSP)⁵⁷, we found that 57% of the functional enhancers contained AD-associated variants in at least one ADSP analysis (Methods). One such enhancer, Enh427 (chr17: 34,252,507–34,253,084), regulates CCL2 (Extended Data Fig. 4), a chemokine dysregulated in AD astrocytes and associated with AD neuropathology (Fig. 5b), which recruits microglia and monocytes to amyloid-β plaques, thereby amplifying the inflammatory response and contributing to neuronal damage^58,59. The CCL2 enhancer contains three nominally AD-associated variants (Fig. 6d), and the variant with the most significant association with AD exhibited the most significant DIS (Fig. 6d), demonstrating the added value of variant-level prioritization by in silico mutagenesis. Supplementary Table 8 provides DIS data for all functional enhancers, and the in silico mutagenesis data are available at https://voineagulabunsw.github.io/astrocyte_crispri_resource/↗.

To expand the enhancer–gene regulatory map, we evaluated four predictive models (ABC^4,60, TAP-seq RF¹⁵, ENCODE’s rE2G-extended and rE2G-DNase) that outperform a broad range of deep-learning and standard machine learning models¹⁴. Because these models were developed using data from K562 cells^4,14,15,61, we assessed their performance on both astrocyte and K562 CRISPRi screen datasets as ground truth.

We then developed an RF model of enhancer–gene regulation in astrocytes (EGrf), using the following seven variables: three enhancer chromatin state features (ChIP H3K27ac, ChIP H3K4me3, and ATAC–seq), gene stability index³⁶, enhancer–gene distance, the number of transcription start sites (TSSs) between the gene and enhancer, and a binary indicator of whether the gene in the enhancer–gene pair is the closest gene to the enhancer. Of these, the gene stability index and the closest gene variables had not been included in the other predictive models we evaluated.

EGrf achieved an AUPRC of 0.57 and 35% precision at 70% recall, a higher predictive accuracy compared to the rE2G-DNase model in astrocytes (P = 0.04; Fig. 7a,b). While several other variables were predictive on their own (enhancer transcription, Hi-C contacts, maximum DeepSEA-based variant scores), their inclusion in the EGrf model did not improve the performance alongside the existing seven variables (Extended Data Fig. 9). Notably, replacing ATAC–seq with DNase-seq significantly reduced the model’s performance (Extended Data Fig. 9). In Supplementary Note 2, we describe the contribution of each variable to EGrf, and introduce an RF model to predict EG pairs in K562 cells, which performs as well as rE2G-extended, while using fewer variables.

We then applied EGrf to intergenic ATAC–seq peaks not tested by CRISPRi to assess their effects on expressed genes within a 500 kb region. At a stringent threshold of 80% precision (prediction score = 0.535), EGrf-predicted 1,348 enhancer–gene regulatory interactions (Supplementary Table), involving 1,223 unique enhancers and 805 unique genes, substantially expanding the catalog of predicted functional regulatory interactions in astrocytes. 9

When assessing the properties of the predicted EGPs, we found that they were enriched for enhancers previously identified by CRISPRi and MPRA⁵⁸, and for regions carrying superenhancer marks (Fig. 7c). The median enhancer–gene distance was 26.69 kb. The predicted enhancer-regulated genes were depleted of housekeeping genes, while the predicted functional enhancers were much more likely to be transcribed as measured by TT-seq (P = 8.1 × 10⁻²¹⁷, OR = 6.6), supporting the validity of the predicted EGPs.

We then carried out individual CRISPRi experiments for five enhancers predicted by EGrf to control genes implicated in AD^58,62–64 and neurodegeneration⁶⁵, including an enhancer controlling CLU, a major risk gene for late onset AD. Each enhancer was targeted by two sgRNAs and each sgRNA was tested in duplicate, with the expression of the predicted target gene assessed by qRT–PCR. These data validated four of the five predicted regulatory interactions (Fig. 7d), consistent with EGrf’s 0.8 precision.

Together, our data establish a framework for functional characterization of enhancer–gene regulatory networks in a given cell type, through CRISPRi screening followed by genome-wide prediction using an RF approach.

Additional details on experimental methods are provided in: characterization of NHA and NHA–dCas9–KRAB cells (genotyping, ATAC-seq, RNA-seq and TT-seq library preparation and sequencing) and validation experiments. Supplementary Methods

NHA cells were transduced with Lenti–dCas9–KRAB–blast (a gift from G. Hon; Addgene, 89567) and selected with 10 µg ml⁻¹ blasticidine S hydrochloride (Sigma-Aldrich, 15205-25MG) after 48 h to create a stable cell line expressing dCas9–KRAB.

CROP-seq-opti-GFP was cloned by inserting the sequence for green fluorescence protein (GFP) and a T2A linker connected to the puromycin sequence in CROP-seq-opti-GFP (a gift from J. Shendure; Addgene, 106280). A synthetic gBlock (IDT) was generated with the sequence for GFP and a T2A linker, with flanking regions to recapitulate the CROP-seq-opti lost due to restriction enzyme digestion (Supplementary Table 10). The gBlock and CROP-seq-opti-GFP were digested with BsiWI (New England Biolabs, R3553S). The digested plasmid was dephosphorylated with Shrimp Alkaline Phosphatase (New England Biolabs, M0371S), and ligated with the digested gBlock using T4 DNA ligase (New England Biolabs, M0202S). The ligated plasmid was transformed into NEB 5-α competent Escherichiacoli (high efficiency) (New England Biolabs, C2987H) and grown on LB agar plates (UNSW; Lowy Solutions, UCS-BIO-0038) overnight at 30°C with 100 µg ml⁻¹ ampicillin (Sigma-Aldrich, A5354-10ML). Colonies were screened by culturing picked colonies in LB broth (UNSW; Lowy Solutions, UCS-BIO-0046) with 100 µg ml⁻¹ ampicillin overnight at 30°C and isolating plasmids using the QIAprep Spin Miniprep Kit (Qiagen, 27106) before Sanger sequencing across the cloned region (Ramaciotti Centre for Genomics, UNSW; Supplementary Table 10). CROP-seq-opti-GFP is available upon request.

Additional details on data analysis methods are provided in Supplementary Methods:Bulk ATAC–seq, RNA-seq and TT-seq data processingSNP genotyping and copy number variation callingNanoString data analysisProcessing and analysis of data from ref. ¹⁷Functional annotation of genesFunctional annotation of enhancersAnnotation of SNPs overlapping candidate enhancerseQTL data processingTF footprinting analysisRegulatory network constructionProcessing and analysis of TF–ChIP–seq dataDeep-learning-based in silico mutagenesis and DISsProcessing of WGS data from the ADSPPredictive models of enhancer–gene regulatory interactionComparison of chromatin state at intergenic regions across ENCODE astrocyte linesLarge language models

Candidate NHA regulatory regions were selected by filtering NHA ATAC–seq peaks on five criteria (Extended Data Fig. 2):ATAC–seq read count (that is, MACS2 Pileup) >2 counts per million (CPM).Intergenic, that is, at least 2,000 bp from gene bodies (GENCODE, v32).Overlapping with the PsychENCODE consortium’s list of brain-active regulatory regions¹³.Overlapping with a TAD in cultured induced glial cells²², where the TAD contained at least one gene expressed above a Seurat-normalized scRNA-seq expression of 0.3 (log(((counts/total count) × 10,000) + 1)). This expression threshold was selected based on a small-scale pilot experiment assessing the minimal gene expression required for detecting silencing induced with CRISPRi (Extended Data Fig. 2).Overlapping with peak calls from at least two of the four ATAC–seq datasets generated from the following: cultured hiPSC-derived astrocytes²¹; GFAP⁺ or HEPACAM⁺ cells isolated from cortico-spheroids²⁰; NeuN⁻ nuclei isolated from the prefrontal cortex¹⁹; and NeuN⁻ nuclei isolated from the dorsolateral prefrontal cortex, ventrolateral prefrontal cortex and orbitofrontal cortex¹⁸. All datasets were filtered to include peaks with pileup counts per >2 million.

All intersections were performed using intersectBed from the Bedtools suite⁷⁷, requiring at least 1-bp overlap with no window. Where necessary, peak coordinates were lifted over to hg38 using rTracklayer⁷⁸ with the UCSC hg19-to-hg38 chain (https://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz↗).

FASTQ demultiplexing, mapping and quantification of gene expression were performed using CellRanger (v6.0.2; ref. ⁷⁹). Default parameters were followed, with mapping to the CellRanger 2020-A transcriptome reference (hg38 and GENCODE32).

Cells were called using the EmptyDrops package⁸⁰, with a threshold of an FDR of <0.01. Data were processed using Seurat (v4.3.0.1; ref. ⁸¹). To remove potential doublets, we excluded cells whose total UMIs were above the 99.8th percentile within their library or the 99th percentile of all cells with only one assigned gRNA. Low-quality cells were identified and removed when mitochondrial genes, ribosomal genes or MALAT1 accounted for >20% of total UMIs, or dropletQC⁸² annotated them as damaged (Supplementary Fig. 1).

Associations between cell barcodes and guide RNAs were determined using the get_barcodes Python script (https://github.com/shendurelab/single-cell-ko-screens↗) from ref. ⁸³. This script was applied to BAM files from both the 10x libraries and hemi-nested PCR-enriched sequencing runs. Default parameters were used, but with the following flags: --no_swalign --barcode_length 16, and --seach_seq GTGGAAAGGACGAAACACC. Cell–gRNA assignments required a minimum of three UMIs to be called, while also accounting for >2% of total gRNA UMIs within that cell.

scRNA-seq data generated in the CRISPRi screen were annotated to validate the cellular and developmental origin of NHAs. Cell-type annotation was based on normalized cell-level expression data, which were generated using the Seurat pipeline—log(((counts/total count) × 10,000) + 1)⁸¹. The SingleR algorithm (v1.8.1; ref. ⁸⁴), which provides cell-type annotation based on external reference data using an unsupervised correlation-based algorithm with iterative tuning, was applied with default parameters. The reference dataset was constructed using an snRNA-seq time course of human brain maturation¹⁷. Cells were pseudobulked across batches according to their broad cell type (excitatory neurons, inhibitory neurons, astrocytes, oligodendrocytes, oligodendrocyte precursor cells and microglia), and developmental stage (fetal, neonatal, infancy, childhood, adolescence and adult). Pseudobulked data were converted to CPM and log transformed (offset of +0.5).

DE was performed using the improved negative binomial regression model implemented in SCEPTRE²⁵ (v0.1.0), an algorithm designed for high MOI, single-cell-level CRISPR screens. DE analyses tested the effect of silencing each enhancer on all expressed genes with a TSS within ±500 kb of the enhancer midpoint, as in previous studies^6,85. Genes were considered expressed provided that the mean Seurat-normalized expression was greater than 0.015625 (2⁻⁶), equivalent to ~1.64 CPM and the gene being detected in at least ~5.6% of cells (Supplementary Fig. 2). This threshold was chosen based on the distribution of mean gene expression in the scRNA-seq library (Supplementary Fig. 2). Because some TADs were longer than 500 kb, 22 enhancers had no expressed gene within 500 kb, and thus 957 enhancers were tested in the DE analyses.

We thus tested the effect of 957 enhancers on 3,125 genes, for a total of 7,759 EGPs tested. DE analyses were performed at the enhancer level, comparing the targeted group (that is, cells expressing any guide targeting the enhancer) with the nontargeting group (that is, all other cells in the library with at least one other guide) as previously described⁶. Guide-level DE analyses were carried out for negative control sgRNAs, comparing each to all unique genes tested in the enhancer-level analyses. Positive control analyses were conducted by pooling both guides for the given TSS and testing against only the TSS’s gene, noting that one TSS was excluded due to a lack of cells expressing its guides (Pos_EIF3A, with only one cell expressing either guide).

The improved negative binomial test implemented in SCEPTRE was applied with default parameters, as a two-sided test and with the number of resamplings set to 1,000. The following cell-level covariates were adjusted for: log total UMIs; log MOI (where MOI is defined as the number of distinct gRNAs detected in each cell); sequencing batch; percent of UMIs aligning to mitochondrial genes; percent of UMIs aligning to ribosomal genes; nuclear fraction, as calculated using DropletQC⁸²; and the first principal component of cell-cycling gene expression, determined using Tricycle⁸⁶.

To avoid P-value inflation in the EGP tests, we considered the following two approaches: SCEPTRE’s conditional randomization test²⁵ and the empirical P-value correction proposed in ref. ¹—P_empirical = ((number of negative control guide–gene pairs with a smaller P value) + 1)/((the total number of negative control guide–gene pairs + 1)). Multiple testing correction was then applied using the Benjamini–Hochberg FDR approach⁸⁷. Significant EGPs were defined at FDR < 0.1 as previously described^1,25. We found the SCEPTRE approach to have lower sensitivity, identifying 90 significant EGPs versus 158 identified using the empirical P-value approach, of which 88 overlapped (Supplementary Table 3). To determine if the EGPs identified using the approach mentioned in ref. ¹, not the SCEPTRE conditional randomization test, were valid, we included two such enhancers in our validation experiments, and found both to be experimentally validated (Enh219 and Enh220 regulating NEAT1; Fig. 5c). Thus, we selected the empirical P value with the FDR correction approach for our final result set (Supplementary Table 3). We do acknowledge that this more permissive set may, overall, possess a higher false positive rate and thus also specify the more stringent set in Supplementary Table 3 by including the SCEPTRE results.

All statistical analyses were performed in R (v4.1.2).

No statistical methods were used to predetermine sample sizes. However, a power simulation was used post hoc to remove low-powered enhancers, genes or EGPs from statistical background distributions in bioinformatic analyses.

Enhancers were excluded from the CROP-seq scRNA-seq analyses if no expressed genes were detected within 500 kb (n = 22). In qRT–PCR, where technical replicates had a s.d. of >0.4, the replicate with a Ct with the largest difference from other replicates was discarded as an outlier if its removal reduced the s.d.

Data from the CROP-seq scRNA-seq were assumed to follow a negative binomial distribution. Data from qPCR and log₂-transformed NanoString experiments were assumed to be normally distributed. For all Fisher’s exact tests, both variables were categorical, independent and mutually exclusive.

Randomization is inherent to the CROP-seq experimental design, as single cells are randomly transduced with sgRNAs from the pooled sgRNA library. For all other experiments involving NHAs, randomization was not performed as all cells were derived from the same parent population.

Major statistical findings from the CROP-seq experiment assessed for robustness by using two algorithms (SCEPTRE and the empirical FDR), with validation of key results sought using multiple independent orthogonal approaches (for example, NanoString and qPCR). Where possible, multiple independent resources were consulted from the literature (for example, four resources of brain eQTLs and three lists of genes dysregulated in astrocytes in AD).

Data collection and analysis were not performed blindly to the experimental conditions.

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary