A map of enhancer regions in primary human neural progenitor cells using capture STARR-seq

Abstract

Genome-wide association studies (GWASs) and expression analyses implicate noncoding regulatory regions as harboring risk factors for psychiatric disease, but functional characterization of these regions remains limited. Here, we perform capture STARR-sequencing of over 70,000 candidate regions to identify active enhancers in primary human neural progenitor cells (phNPCs). We select candidate regions by integrating data from NPCs, prefrontal cortex, developmental timepoints, and GWASs. Over 8000 regions demonstrate enhancer activity in the phNPCs, and we link these regions to over 2200 predicted target genes. These genes are involved in neuronal and psychiatric disease-associated pathways, including neuronal system, nervous system development, and developmental delay. We functionally validate a subset of these enhancers using mutation STARR-sequencing and CRISPR deletions, demonstrating the effects of genetic variation on enhancer activity and enhancer deletion on gene expression. Overall, we identify thousands of highly active enhancers and functionally validated a subset of these enhancers, improving our understanding of regulatory networks underlying brain function and disease.

Figure 1.

Panel design and quality control results. (A,B) The process for candidate enhancer selection for Panel 1 (A) and Panel 2 (B). See “Results” for a complete description of candidate selection. (C) The experimental workflow. Sheared genomic DNA was hybridized to probes specific for the candidate enhancer regions. These regions were then cloned into the STARR-seq plasmid and transfected into phNPCs. (D,E) The fold change correlation between the two technical replicates for Panel 1 (D) and Panel 2 (E). Pearson's r² values are included on the graphs. (F,G) Volcano plots representing the tested enhancer regions. Regions that had significant peaks as determined by STARRPeaker are in dark blue and nonsignificant regions are in light blue. (PFC) prefrontal cortex, (PEC) PsychENCODE Consortium, (BP) bipolar disorder, (SZ) schizophrenia, (GWAS) genome-wide association study, (phNPC) primary human neural progenitor cell, (PCW) postconception week, (FDR) false discovery rate, (FC) fold change.

Figure 2.

Expression of transcription factors (TFs) with enriched binding site motifs. Single-cell RNA sequencing expression score is displayed on the y-axis and measured in counts per million (CPM). A higher expression score indicates a more highly expressed gene. For each plot, the solid line represents the median expression score whereas the “X” represents the mean expression score. Individual TF expression scores are shown as data points on the box plots. Expression scores were compared using a two-sample, two-tailed t-test.

Figure 3.

Metascape pathway analysis of the 2288 predicted target genes. (A) Pathway and process enrichment results fielded from KEGG Pathway, GO Biological Processes, Reactome Gene Sets, Canonical Pathways, CORUM, WikiPathways, and PANTHER Pathway. P-values were calculated using the cumulative hypergeometric distribution (Zar 1999). (B) Associated diseases as identified through DisGeNET (Piñero et al. 2017). Enriched cell and tissue types are depicted in C (Subramanian et al. 2005) and D (Pan et al. 2013).

Figure 4.

MutSTARR-seq results comparing enhancer activity (log₂fc) between the reference and alternate alleles. Box plots represent the distribution of activity across four technical replicates. Enhancer activity is defined by log₂ fold change, which represents the normalized output/input ratio in log₂ space. Variants that had a significant effect on enhancer activity (P < 0.05) are boxed in green. P-values were calculated using a t-test comparing the log₂fc(output/input) of the reference and alternate alleles across the replicates. (log₂fc) log₂ fold change, (SNV) single nucleotide variant.

Figure 5.

CRISPR-Cas9 enhancer knockout (KO). (A–D) Left: DNA agarose gel image of the genotyping PCR results after KO of candidate enhancers in phNPCs through ribonucleoprotein (RNP)-mediated CRISPR-Cas9 genome editing. Control cells undergoing the same electroporation without any RNPs showed a clear strong wild-type (WT) band. For enhancer KO samples, besides the higher WT band, there is an additional clearly visible lower band in both BR1 and BR2 samples. The sizes of these lower bands are the same as the expected size of genome edited bands after enhancer KO. (A–D) Right: TaqMan qPCR showed diminished expression level of the target gene after enhancer KO. C_T values from triplicates were used to calculate the expression of the target gene relative to control cells using the Pfaffl method. Averages of the BR1 and BR2 and standard deviations are shown as error bars. (A) Left: EH37E1198822 KO genotyping PCR result (WT band: 3375 bp, genome edited band: 850 bp). Right: Relative expression level change of the target gene NGEF after enhancer KO. (B) Left: EH37E1000386 genotyping PCR result (WT band: 1938 bp, genome edited band: 809 bp). Right: Relative expression level change of the target gene RORB after enhancer KO. (C) Left: EH37E0114246 genotyping PCR result (WT band: 1849 bp, genome edited band: 1205 bp). Right: Relative expression level change of the target gene PLEKHO1 after enhancer KO. (D) Left: EH37E0426064 genotyping PCR result (WT band: 3069 bp, genome edited band: 2203 bp). Right: Relative expression level change of the target gene TOM1L2 after enhancer KO. (Cells, +e) cells without RNPs that underwent the same electroporation and served as control, (NTC) PCR no-template control, (BR) biological replicate.