Enhancer RNAs predict enhancer–gene regulatory links and are critical for enhancer function in neuronal systems

Abstract Genomic enhancer elements regulate gene expression programs important for neuronal fate and function and are implicated in brain disease states. Enhancers undergo bidirectional transcription to generate non-coding enhancer RNAs (eRNAs). However, eRNA function remains controversial. Here, we combined Assay for Transposase-Accessible Chromatin using Sequencing (ATAC-Seq) and RNA-Seq datasets from three distinct neuronal culture systems in two activity states, enabling genome-wide enhancer identification and prediction of putative enhancer–gene pairs based on correlation of transcriptional output. Notably, stimulus-dependent enhancer transcription preceded mRNA induction, and CRISPR-based activation of eRNA synthesis increased mRNA at paired genes, functionally validating enhancer–gene predictions. Focusing on enhancers surrounding the Fos gene, we report that targeted eRNA manipulation bidirectionally modulates Fos mRNA, and that Fos eRNAs directly interact with the histone acetyltransferase domain of the enhancer-linked transcriptional co-activator CREB-binding protein (CBP). Together, these results highlight the unique role of eRNAs in neuronal gene regulation and demonstrate that eRNAs can be used to identify putative target genes.


INTRODUCTION
To orchestrate the precise gene expression patterns that give rise to the phenotypic and functional diversity of complex biological systems, mammalian genomes utilize millions of regulatory elements known as enhancers. Enhancers, often located many kilobases from genes that they regulate, direct transcriptional dynamics at linked genes by activation of proximal gene promoters (1)(2)(3). Enhancer-promoter interactions help to ensure cell-and tissue-specific gene expression profiles in the brain, defining which genes can be turned on during neuronal specification and which genes remain accessible in adult neurons (4)(5)(6)(7)(8)(9). In addition to regulating neuronal development, enhancer regions direct activity-and experience-dependent gene expression programs required for neuronal plasticity, memory formation and behavioral adaptation to environmental stimuli (10)(11)(12)(13)(14)(15)(16)(17)(18). Moreover, the majority of DNA sequence variants that possess a causal relationship to neuropsychiatric disease and intellectual disability fall in non-coding regions of DNA (19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29), and the association between these polymorphisms and altered enhancer function is becoming increasingly clear. Thus, understanding how genomic enhancers regulate individual genes in neuronal systems is critical for unraveling transcriptional contributions to brain health and disease.
Recent advances in DNA sequencing have revealed that the transcriptional landscape of all mammalian organisms is far more complex than previously appreciated. In contrast to earlier predictions, a significant fraction of mammalian genomes is transcribed into non-coding RNAs, which include long non-coding RNAs (lncRNAs; generally defined as non-coding RNAs longer than 200 nucleotides) (30)(31)(32). Much of this lncRNA landscape is composed of enhancer regions which undergo bidirectional, RNA polymerase II (RNAP2)-dependent transcription to yield enhancer RNAs (eRNAs) that are generally not spliced or polyadenylated (5,13,33,34). Critically, RNA synthesis from enhancers that regulate cellular differentiation and stimulus-dependent genes precedes mRNA synthesis from these genes (33). eRNA synthesis also precedes impor-Primary rat neuronal cultures were generated from embryonic day 18 rat cortical, hippocampal, or striatal tissue as described previously (43,45). Briefly, cell culture wells were coated overnight at 37 • C with poly-L-lysine (0.05 mg/ml for culture wells supplemented with up to 0.05 mg/ml Laminin) and rinsed with diH2O. Dissected tissues were incubated with papain for 25 min at 37 • C. After rinsing in Hank's Balanced Salt Solution (HBSS), a single cell suspension of the tissue was re-suspended in Neurobasal media (Invitrogen) by trituration through a series of large to small firepolished Pasteur pipets. Primary neuronal cells were passed through a 100 M cell strainer, spun and re-suspended in fresh media. Cells were then counted and plated to a density of 125 000 cells per well on 24-well culture plate and 250 000 cells per well on 12-well culture plate with or without glass coverslips (60,000 cells/cm). Cells were grown in Neurobasal media plus B-27 and L-glutamine supplement (complete Neurobasal media) for 11 DIV in a humidified CO2 (5%) incubator at 37 • C. [4][5][6][7][8][9][10][11], neuronal cultures were treated as described. For KCl stimulation experiments, KCl (Sigma) was added to complete Neurobasal media, and KCl solution or vehicle (complete Neurobasal media alone) was added for the indicated final concentrations. Cells were incubated with KCl for described time points prior to RNA extraction. For TTX inactivation experiments, cells were treated with 1 M TTX (Tocris Bioscience) in Neurobasal media for the described time points prior to RNA extraction. S-AMPA, NMDA, and FSK (Sigma) were resuspended in sterile water, diluted in Neurobasal media, and added to cultures for 1 h at equal volumes (final concentrations of 1, 10 or 100 M). The same volume of Neurobasal media was added as a vehicle control. For time course experiments, samples were flash frozen on dry ice at the described time points and stored at −80 • C until further processing. For experiments involving RNAP inhibitors, cultures were treated for 4 or 4 h followed by a 1 h, 25 mM KCl stimulation. The RNAP2-dependent transcriptional inhibitor DRB (Sigma) was dissolved to a 20 mM stock solution in 100% cell culture grade DMSO (Sigma) and diluted in Neurobasal media to described experimental concentrations. For CREB inhibitor experiments CREBi (666-15, Torcis) also called 3-(3-Aminopropoxy)-N- [2-[[3-[[(4-chloro-2-hydroxyphenyl)amino]carbonyl]-2-naphthalenyl]oxy]ethyl]-2-naphthalenecarboxamine hydrochloride was dissolved to a 10 mM stock solution in 100% cell culture grade DMSO (Invitrogen) and diluted in Neurobasal media to for a final treatment concentration of 1M. In experiments involving DRB and 666-15, vehicletreated cells received equal concentrations of DMSO in Neurobasal media.
For viral transduction, cells were transduced with lentiviruses on DIV 4 or 5 (only viruses with a minimum titer of 1 × 10 9 GC/ml were used for target multiplicity of infection (MOIs) of at least 1000). After an 8-16 h incubation period, virus-containing media was replaced with conditioned media to minimize toxicity. A regular half-media change followed on DIV 8. On DIV 11, transduced cells were imaged and virus expression was verified prior to KCl-treatment and/or RNA extraction. All samples were frozen immediately and stored at −80 • C until further processing. Immunocytochemistry for FLAG was performed as described previously (43) with an anti-FLAG antibody (MA1-91878, Thermo Fisher Scientific, RRID AB 1957945). EGFP and mCherry expression was also used to visualize successful transduction using a Nikon TiS inverted epifluorescence microscope.
Nuclei from rat embryonic cortical, hippocampal, or striatal neurons (50 000/region) were used for ATAC-seq library preparation, following a modified protocol (46)(47)(48). Briefly, 25 l of 0.4 M KCl (Sigma) or vehicle (Neurobasal media) was added to cell culture wells to achieve a final concentration of 10 mM KCl and incubated at 37 • C for 1 h. Following treatment, media containing KCl or Vehicle was aspirated, and cells were washed with 1× cold PBS. Then, 600 l of lysis buffer (10 mM Tris-HCl, 10 mM NaCl, 3 mM MgCl 2 , 0.1% IGEPAL CA-630, Molecular Grade H 2 O) was added to cell culture wells and incubated for 5 min on ice. Lysed cells were then transferred to a 1.5 ml Eppendorf tube and centrifuged at 500 g for 10 min in a swinging bucket centrifuge. Following removal of supernatant containing lysis buffer and cell debris, 500 l of 1× cold PBS was added to each tube, and nuclei were counted using the Countess II (Life Technologies). The approximate volume needed to achieve 50 000 nuclei was then placed in new Eppendorf tube, and centrifuged at 500 g for 10 min in a swinging bucket centrifuge. The nuclei pellet was then resuspended in 22.5 l of tagmentation reaction mix containing Molecular Grade H2O, 2× Tagment DNA Buffer (Illumina), 0.011% Digitonin, 0.1% Tween-20, followed by the addition of 2.5 l TDE1 (Illumina) and incubated at 37 • C for 1 h. Following tagmentation, libraries were PCR amplified and purified using the Qiagen Minelute PCR purification kit (Qiagen). To remove primer dimers, libraries were bead purified using one round of 1.0× Ampure XP beads with 80% ethanol, followed by a second round of 0.8× Ampure beads with 80% ethanol. Libraries were sequenced (75-bp paired end reads) on the NextSeq500 platform (Illumina).

ATAC-Seq mapping and peak calling
Paired-end FASTQ files were aligned and trimmed using a custom bioinformatic pipeline initialized in snakemake v5.3.0 (49). Briefly, low-quality bases (Phred < 20) and Nextera adapters (5 -CTGTCTCTTATA-3 ) were identified and trimmed using FastQC and TrimGalore (v0.4.5), respectively. FASTQ files containing trimmed sequences were then aligned to Rn6 Ensembl genome assembly (v95) to generate binary alignment map (BAM) files with Bowtie2 (v.2.3.4.2) with custom options:`-local -very-sensitive-local`, and`-X 3000` (50). Sequences were then ordered by genomic position using Samtools (v.1.9) (51). Next, BAM files for each brain region were merged to generate three metasamples that were used for downstream data analysis. Peaks for each region were called using MACS2 (52) callpeak with the options --qvalue 0.00001 --gsize 2729862089 --format BAMPE. These options utilize the default behavior of the MACS2 algorithm to ignore duplicates. Within each brain region, peaks closer than 1000 bp were merged with BEDtools (53) and peaks less than 146 bp, the specific length of DNA wrapped around a single nucleosome, were removed in R. Finally, peaks from each brain region were merged with BEDtools to create an additive peak set containing 192 830 peaks.

Bulk RNA-sequencing
Bulk RNA-sequencing (RNA-seq) was carried out at the Heflin Center for Genomic Science Genomics Core Laboratories at the University of Alabama at Birmingham. RNA was extracted, purified (RNeasy, Qiagen), and DNasetreated for three to four biological replicates per brain region and experimental condition. 1 g of total RNA underwent quality control (Bioanalyzer) and was prepared for directional RNA sequencing using NEBNext reagents (New England Biolabs) according to manufacturer's recommendations. Specifically, the NEBnext rRNA depletion kit was used to remove ribosomal RNA to include both polyadenylated and non-polyadenylated RNA. RNA-seq libraries underwent sequencing (75 bp paired-end directional reads; ∼12.4-31.1M reads/sample) on an Illumina sequencing platform (NextSeq500).

Bulk RNA-seq data analysis
Paired-end FASTQ files were uploaded to the University of Alabama at Birmingham's High-Performance Computer cluster for custom bioinformatics analysis using a pipeline built with snakemake (49). Read quality was assessed with FastQC and low-quality bases (Phred < 20) and Illumina adapters were trimmed with TrimGalore (v0.4.5). Spliceaware alignment to the Rn6 Ensembl genome assembly (v95) was performed with STAR (v2.6.0; (54)). Binary alignment map (BAM) files were merged and indexed with Samtools (v1.9).

HOMER motif analysis
Region-specific TAPE genomic positions were compiled into three separate BED files. The findMotifsGenome.pl function within the HOMER (v4.11.1) package was used to identify enriched motifs and their corresponding transcription factors within these genomic positions with optionssize 1000 -len 8,10,12 -mask -preparse -dumpfasta. As we were interested in the motifs that were specific to each region, a BED file containing peaks from the two other regions was used as background. For example, when identifying enriched motifs for striatal-specific TAPEs, a background file containing genomic positions for cortical-and Nucleic Acids Research, 2020, Vol. 48, No. 17 9553 hippocampal-specific TAPEs. Furthermore, only those denovo motifs not marked as possible false positives are reported.

ROC & TAPE identification
General enhancer identification was performed on merged ATAC-seq and total RNA-seq BAM files from all replicates for Veh and KCl treated samples from primary cultures generated from cortical, hippocampal and striatal tissue. A complete data analysis and replicate handling pipeline, with all custom code, is provided in Supplementary Methods. MACS2 -defined ATAC-seq peaks from main autosomal, X or Y chromosomes were used to generate a list of 191 857 regions of open chromatin (ROCs) spanning ±500 bp of identified ATAC-seq peaks. Due to the difficulty in separating intronic enhancers from potential promoters or other elements, we used Seqmonk software (Babraham Institute) to identify ROCs that fall >1kb outside of genes curated by Refseq, UCSC and Ensemble. To account for unannotated or misannotated genes, we also filtered out ROCs that overlapped contiguously transcribed regions (>100 bp in length, merging regions closer than 1 kb) or known noncoding RNAs such as miRNA, rRNA, snoRNA, snRNA and tRNA. These filtering steps provided a list of 100 767 intergenic ROCs (iROCs). Bidirectional transcription at iROCs was calculated in Seqmonk by quantification of strand-specific RNA counts from merged total RNA-seq files from all replicates and culture systems, defined as the number of forward reads as a percentage of total reads. iROCs containing between 5% and 95% of total RNA-seq reads originating from the forward strand were considered bidirectional. This filtering produced a list of 28,492 transcriptionally active putative enhancers (TAPEs) with potential bidirectional transcription. Possible TAPE-gene pairs were identified by mapping all gene promoters within 1 Mb upstream or downstream from the center of the TAPE. CPKM values for each TAPE and associated gene were correlated using a Pearson's correlation in R. Pearson's correlations were calculated using all samples (to yield a global correlation across culture systems), as well as using only samples from each region-specific culture system (to generate region-specific correlation values). TAPE-gene pairs with global correlations of NA were removed as these values were generated when the corresponding gene contained count values of 0 for every sample. Removing these pairs left 388 605 potential TAPE-gene pairs. TAPE-gene pairs with correlations >0.5 were deemed high-confidence pairs. These high-confidence pairs were then used to investigate distance and gene position distributions. Region-specific TAPEs were identified in Seqmonk using DESeq2 to compare TAPE RNA counts in cortical, hippocampal, and striatal primary neuron cultures (FDR < 0.05). DESeq2 identified 2467 potential region-selective TAPEs. However, this list includes TAPEs which may be highly expressed in two of the three culture systems. To ensure region selectivity, we conducted hierarchical clustering on TAPE RNA values from these 2467 TAPEs (z-scored across region) in Seqmonk to identify nodes of TAPEs specific to each region. This analysis identified 776, 390 and 898 region-selective TAPEs for hippocampus, cortex and striatum, respectively.
To identify whether identified TAPE loci correspond to other marks of active enhancers, we used publicly available ENCODE project datasets of histone modifications in the mouse forebrain at postnatal day 0. ENCODE datasets were obtained from the UCSC Table Browser (H3K4me1,  ENCFF172LKQ; H3K4me3, ENCFF022PPU; H3K27ac,  ENCFF348JQQ; H3K27me3, ENCFF915VHI). Mouse genome coordinates (mm10) were transformed to rat genome coordinates (Rn6) using LiftOver. CTCF binding motifs were identified across the Rn6 genome build by searching the CTCF consensus motif (CCGCGNGGNG GCAG) in MEME suite. The relationship between these datasets and identified TAPEs was queried using Seqmonk to align identified peak locations or motifs to the center of all TAPEs, and graphed as counts per million.

CRISPR-dCas9 and shRNA construct design
To achieve transcriptional activation, lentivirus-compatible plasmids were engineered to express dCas9 fused to VP64 or VPR constructs (Addgene plasmid # 114196 (55)). dCAS9-VP64 GFP was a gift from Feng Zhang (Addgene plasmid # 61422 (56)). The pcDNA-dCas9-p300 Core construct was a gift from Charles Gersbach (Addgene plasmid # 61357 (57)). VP64-and VPR-expressing constructs were co-transduced with sgRNA-containing constructs. A guide RNA scaffold (a gift from Charles Gersbach, Addgene #47108 (58)) was inserted into a lentivirus-compatible backbone, and EF1␣-mCherry was inserted for live-cell visualization. A BbsI cut site within the mCherry construct was mutated with a site-directed mutagenesis kit (NEB). Gene-specific gRNAs were designed using an online sgRNA tool, provided by the Zhang Lab at MIT (crispr.mit.edu). To ensure specificity, all CRISPR crRNA sequences were analyzed with the National Center for Biotechnology Information's (NCBI) Basic Local Alignment Search Tool (BLAST). sgRNAs were designed to target Fos, Fosb and Nr4a1 enhancers, respectively, as well as the promoter and control regions (a list of the target sequences is provided in Supplementary Data Table S4). cr-RNA sequences were annealed and ligated into the sgRNA scaffold using the BbsI or BsmBI cut site. For CRISPR-Display, lentiCRISPR v2 from Feng Zhang (Addgene plasmid # 52961 (59)) was modified and engineered to express dCas9 (instead of Cas9) and GFP under an hSYN promoter, as well as additional restriction sites (Esp3I) for subsequent acRNA insertion using PacI and Xbal. sgRNA, and acRNA sequences of eRNA1, eRNA3 and control RNA were inserted via restriction enzyme cloning using gBlocks for acRNA insertion (cut with Esp3I). As another control, a plasmid lacking the eRNA sequence was targeted to the same genomic sites. To achieve RNA knockdown, shRNA sequences targeting the gene of interest were designed using the Broad TRC shRNA design tool (http://portals.broadinstitute.org/gpp/public/) according to the Addgene pLKO.1 protocol (https://www.addgene. org/tools/protocols/plko/) and inserted into a lentiviruscompatible shRNA construct. Briefly, a pLKO.1-TRC vector (a gift from David Root; Addgene plasmid #10879; RRID:Addgene 10879) (60) was cloned into an expression vector with an mCherry reporter (Addgene plasmid #114199) (55) to generate distinct U6-shRNA and EF1a-mCherry expression cassettes. Oligonucleotides containing the shRNA sequence were cloned into this backbone using AgeI and EcoRI restriction digest.
To ensure specificity all shRNA sequences were analyzed with BLAST. All targeting and scrambled control shRNA sequences were annealed and ligated into the shRNA cassette-containing construct using the AgeI and EcoRI cut sites. Plasmids were sequence-verified with Sanger sequencing; final crRNA insertion was verified using PCR.

Lentivirus production
Viruses were produced in a sterile environment subject to BSL-2 safety by transfecting HEK-293T cells with specified CRISPR-dCas9 plasmids, the psPAX2 packaging plasmid, and the pCMV-VSV-G envelope plasmid (Addgene 12260 & 8454) with FuGene HD (Promega) for 40-48 h as previously described (61). Viruses were purified using filter (0.45 m) and ultracentrifugation (25 000 rpm, 1 h 45 min) steps. Viral titer was determined using a qPCR Lentivirus Titration Kit (Lenti-X, qRT-PCR Titration Kit, Takara). For smaller scale virus preparation, each sgRNA plasmid was transfected in a 12-well culture plate as described above. After 40-48 h, lentiviruses were concentrated with Lenti-X concentrator (Takara), resuspended in sterile PBS, and used immediately. Viruses were stored in sterile PBS at −80 • C in single-use aliquots.

Antisense oligonucleotide (ASO) design and treatment
To manipulate Fos mRNA or eRNA levels, we designed 20 bp ASOs that targeted distinct transcripts from the Fos gene locus (see Supplementary Data Table S4 for target sequences). ASOs targeting exon 3 of Fos mRNA or Fos eRNA1 were synthesized with two chemical modifications: an all phosphorothioate backbone and five 2 Omethyl RNA bases on each end of the oligonucleotide (Integrated DNA Technologies). Primary neuronal cultures were treated with scrambled or Fos targeted ASOs (15 M in buffer EB, for a final concentration of 1.5 M) and incubated for 24 h (basal experiments) or 23 h followed by 1 hr neuronal depolarization with 25 mM KCl (or vehicle control). Following ASO treatment, RNA was extracted (Qiagen RNeasy kit) and Fos mRNA and eRNA levels were determined using RT-qPCR with custom primers.

RNA single molecule fluorescent in situ hybridization (sm-FISH)
smFISH probe design. We designed and ordered Stel-laris® FISH probe sets for Gapdh mRNA, Fos eRNA1, Fos eRNA3 and Fos mRNA carrying a fluorophore (Quasar ® 570 for eRNA1 and Gapdh mRNA probes, Quasar ® 670 for both Fos and Gapdh mRNA probes). We preferred probes of 20-mer oligonucleotides. Multiple probes per set targeting the same RNA molecule were designed for an adequate signal to background ratio and to optimize signal strength. Target sequences of each probe set are provided in Supplementary Data Table S4).
Sample preparation and hybridization. Day 1: Primary neuronal cultures (∼250 000 neurons per coverslip/well) were KCl-or vehicle-treated for 1 h on DIV 11. After treatment cells were cross-linked with 3.7% formaldehyde (paraformaldehyde in 1× PBS) for 10 min at room temperature (21 • C) on a rocking platform. Wells were washed twice with PBS and permeabilized in 70% ethanol for at least 3 h at 4 • C. Wells were washed in Stellaris® Wash Buffer A with for 5 min at room temperature. Coverslips were transferred to a humidifying chamber and incubated with hybridization buffer (0.5 nM mRNA probe, 0.5 nM eRNA probe) for 14 h at 37 • C.
Day 2: Coverslips were washed three times in Stellaris® Wash Buffer A for 30 min at 37 • C. After a 5 min wash in Stellaris ® Wash Buffer B at room temperature, coverslips were mounted using ProLong™ antifade with DAPI for imaging.

Quantification of expression.
A number of freely available programs have been developed to quantify smRNA FISH results. We used StarSearch (http://rajlab.seas.upenn.edu/ StarSearch/launch.html), which was developed by Marshall J. Levesque and Arjun Raj at the University of Pennsylvania to automatically count individual RNAs. mRNA and eRNA detection involved two major steps. First, images for each probe set as well as a DAPI image are merged and cells were outlined. Punctae detection was carried out and additional adjustment of thresholds was performed. The same threshold range was used for all images, and this analysis was performed blind to treatment group. As a negative control, we quantified processed samples without FISH probes to determine non-specific background signals. Background signal for the Quasar ® 570 channel (which was used to Nucleic Acids Research, 2020, Vol. 48,No. 17 9555 image Fos eRNA1 and Gapdh mRNA) was close to zero (0.3761 ± 0.07906 spots/cell). We did not detect any background spots in the Quasar ® 670 channel, which was used to image Fos mRNA and Gapdh mRNA.

RNA electrophoretic mobility shift assay (REMSA)
Mobility shift assays were conducted with synthetic Fluorescein-labeled ∼150-base RNA oligonucleotides and specified concentrations of recombinant CBP HAT domain (Sigma, 1319-1710) and CBP bromodomain (Abcam, ab198130). Oligonucleotides were synthesized using acRNA-containing Display plasmids as template DNA and the HiScribe T7 High Yield RNA synthesis kit (NEB) with Fluorescein-12-UTP (Sigma). RNA oligonucleotides (eRNA1 46.11 ng, eRNA3 47.19 ng, ctrl RNA 46.51 ng for 20 nM final concentrations; up to 5 M for competition assay) were heated to 95 • C for 5min and refolded at room temperature prior to incubation with CBP protein (0-0.6 M) in REMSA buffer (20 mM HEPES, 40 mM KCl, 1 mM EDTA, 0.2 mM DTT, 0.1 mg ml -1 BSA, 0.1% Tween-20 and 20% Glycerol, and 0.1 mg/ml yeast tRNA (Sigma)) for 1 h at 37 • C. REMSA was performed using 1% agarose gel electrophoresis in 0.5× TBE. Electrophoretic mobility of fluorescein-labeled RNA was assayed using fluorescence imaging on the Azure c600 Imaging System (Azure biosystems). RNA-CBP complex formation was quantified as the Fluorescein signal intensity appearing at the higher band (corresponding to CBP-bound RNA with lower electrophoretic mobility) divided by the total signal intensity (bound RNA plus free probe).

Statistical analysis
Required sample sizes were calculated using a freely available calculator (Lenth, R. V. . Java Applets for Power and Sample Size [Computer software]. Available at http://www.stat.uiowa.edu/~rlenth/Power). Transcriptional differences from PCR experiments were compared with one-way ANOVA with post hoc tests where appropriate, two-way ANOVA with post hoc tests where appropriate, or Student's t-tests/Mann-Whitney tests. Significance of sm-FISH data was assessed with Mann-Whitney test or Pearson correlation test. Statistical significance was designated at ␣ = 0.05 for all analyses. Statistical and graphical analyses were performed with Graphpad software (Prism). For all statistical tests, assumptions (e.g. normality and homogeneity for parametric tests) were formally tested. In cases where these assumptions were violated, non-parametric tests were used to evaluate statistical significance.

Chromatin accessibility and transcription predict enhancer location, activity state and target genes
To map enhancers genome-wide for three different brain regions (primary neuronal cultures generated from cortex, hippocampus and striatum), we took advantage of the fact that active enhancers are associated with an open chromatin structure and are bidirectionally transcribed. We first identified 191 857 regions of open chromatin (ROCs) by generating Assay for Transposase-Accessible Chromatin sequencing (ATAC-Seq) libraries from each primary neuronal culture system. Consistent with common cutoffs used to dissociate enhancers from more proximal promoters, we filtered for intergenic ROCs (iROCs) that fell outside of canonical gene boundaries and >1 kb away from annotated transcription start sites. Next, we capitalized on the characteristic bidirectionality of enhancer transcription and incorporated direction-specific total RNA-seq datasets from the same neuronal culture systems to identify bidirectionally transcribed regions of open chromatin, identifying 28,492 transcribed putative enhancers (TAPEs) ( Figure  1A, top; Supplementary Data Table S1). Confirming our pipeline, these bidirectionally transcribed regions exhibited high chromatin accessibility and dual peaks of RNA expression from each strand marking TAPE centers ( Figure  1A, bottom). Identified TAPEs exhibited characteristic patterns of transcription, chromatin landscape and high transcriptional correlations to potential target genes as shown in the representative examples at the Klf4 and Sik1 loci ( Figure 1B). To further validate our enhancer identification pipeline, we capitalized on publicly available ENCODE datasets from postnatal day zero (P0) mouse forebrain for the major histone modifications commonly used for enhancer identification, including H3K4me1, H3K4me3 and H3K27ac (marks for active or poised enhancers; (2), Figure  1C). TAPE centers were enriched for the histone modifications H3K4me1, H3K4me3 and H3K27ac (but not the repressive H3K27me3 modification), as well as the enhancerlinked chromatin looping factor CTCF (CCCTC-binding factor) motifs. Additionally, TAPE regions exhibited enhanced sequence conservation (PhastConsElements20way) compared to surrounding regions. While many pipelines have been developed to enable putative enhancer identification, the ability to predict which gene or genes are controlled by specific enhancers has remained challenging and often requires many overlapping ChIP-seq, ATAC-seq and HiC-seq datasets for accurate prediction within a cell type. Here, we leveraged the inherent variability in eRNA levels across identified TAPEs to construct pairwise correlation matrices with mRNA estimates from all annotated protein-coding genes falling within 1 Mb of TAPE boundaries. This strategy makes the simple assumption that enhancers and linked genes will correlate in their transcriptional output across different cell classes and/or activity states. By filtering positively correlated, high-confidence TAPE-gene pairs from 433 416 TAPE-gene matches, we defined predicted pairs for further investigation and functional validation (Fig. 1D). Globally, identified TAPEs correlated more strongly with proximal genes than with distal genes, although average correlations were relatively weak (r < 0.05; Figure 1E, F). Likewise, in contrast to traditional enhancer-gene pair prediction pipelines that simply annotate the nearest gene to an identified enhancer, our pipeline demonstrated that only 20.9% of TAPEs exhibited the strongest correlation with the closest gene, whereas 79.1% of all maximal enhancer-gene correlations occurred between TAPEs and more distal genes (ranked gene position; Figure 1F). Importantly, the maximal correlation value of top TAPE-gene pairs did not decrease at more distal target genes relative to their respective  TAPE ( Figure 1G), suggesting the presence of long-range enhancer-gene interactions with many intervening genes.

Selective enhancer subsets are linked to gene expression programs that underlie region-specific development and function
Harnessing region-specific variation in transcription, we next quantified count data at individual TAPEs and used DESeq2 to identify TAPEs that exhibited region-selective expression patterns. This analysis yielded 390 cortical, 776 hippocampal and 898 striatal putative enhancers ( Figure  2A Table S2). Transcription factor binding motif enrichment analysis at TAPES and gene ontology term analysis at predicted TAPE target genes indicated that region-selective TAPEs play integral roles in several region-specific processes ( Figure 2C; Supplementary Data Table S3). Hippocampus-selective TAPEs were correlated with genes linked to cellular and nervous system development and cell-cell signaling, and they display an enrichment of NEUROG2 (Neurogenin2) binding mo-

D E
Prox1 mRNA P56 mouse brain

Allen Brain Atlas
Mn1 mRNA P56 mouse brain Allen Brain Atlas GO Term analysis on region-selective gene targets Enriched motifs regulation of cell development anatomical structure morphogenesis epithelial cell differentiation regulation of nervous system development cell-cell signaling regulation of protein serine/threonine kinase activity forebrain development central nervous system neuron development negative regulation of Ras protein signal transduction central nervous system neuron axonogenesis heart trabecula morphogenesis actin polymerization or depolymerization mitochondrial respiratory chain complex I assembly N-terminal protein amino acid modification nuclear polyadenylation-dependent rRNA catabolic process dTMP biosynthetic process protein maturation by iron-sulfur cluster transfer positive regulation of oxidative stress-induced cell death  Table  S3). (D-F) Genome browser tracks showing ATAC-seq and total RNA-seq (reads mapped to + and -strand) signal at example loci separated by the three brain regions of interest (top, square brackets indicate y-axis max CPM value for all ATAC-seq and RNA-seq tracks at the respective loci). Region-selective TAPEs are represented by Kcnf1 for cortex, Prox1 for hippocampus, and Mn1 for striatum. Distance between TAPE and predicted target gene and global Pearson's correlation is indicated above each representative TAPE-gene pair, Rn6 chromosome coordinates are indicated below. Region-selective expression of the respective target genes is also evident in in situ hybridization images from the adult mouse brain (bottom, Image credit: Allen Brain Atlas, Allen Institute). tifs, a transcription factor crucial for dentate gyrus development (62,63). Similarly, cortex-selective TAPEs correlated with genes implicated in forebrain and CNS development, neuron axonogenesis, and actin polymerization, and exhibit enrichment in binding motifs for SOX5, a transcription factor that regulates neuronal migration and differentiation in the neocortex (64). Likewise, striatum-selective TAPEs correlated with genes important for several amino acid metabolism and modification pathways as well as mitochondrial functions. These TAPEs are enriched in motifs for the transcription factor ISL1, which is required for differentiation of striatonigral pathway projection neurons (65). Not surprisingly, genes corresponding to these TAPEs also demonstrated region-selective expression, as seen in the examples Kcnf1, Prox1 and Mn1 for cortex, hippocampus, and striatum, respectively ( Figure 2D-F). These observations were confirmed by in situ hybridization images obtained from the Allen Brain Mouse Atlas.

Activity-dependent enhancer RNA precedes and predicts mRNA induction
To determine whether eRNAs are correlated with activitydependent alterations in protein-coding genes, we exam-ined RNA transcription from TAPEs following neuronal depolarization with 10 mM potassium chloride (KCl) for 1 h ( Figure 3A). Globally, we identified 96 activity-regulated TAPEs that were significantly altered by KCl treatment in at least one cell type, with 20 selective for cortex, 22 for hippocampus, and 66 for striatum ( Figure 3B; Supplementary Data Table S5). Activity-regulated TAPEs demonstrated either increased or decreased transcription in response to neuronal depolarization, here termed up-or downregulated TAPEs. As expected, we found that mRNA expression of predicted target genes correlated with TAPE transcription ( Figure 3C). However, surprisingly we found that chromatin accessibility (ATAC-seq signal) at either promoter or TAPE regions was not predictive of gene expression levels ( Figure 3C), but often shifted in the opposite direction of TAPE/gene RNA estimates. While the decreased TAPE ATAC signal could be specific to this time point, it is worth noting that enhancer transcription in this case provides better predictions of basal and stimulus-dependent gene expression than chromatin accessibility. Figure 3D-F shows ATAC-seq and RNA-seq results from three representative IEGs (Fos, Fosb and Nr4a1) that are significantly induced by KCl depolarization. Each of these genes displayed distal activity-regulated TAPEs, including at least three distinct enhancers near the Fos gene. The locations of these enhancers are consistent with locations of enhancer elements in other species relative to the Fos gene (12,13) and map to DNA sequences that are enriched for histone modifications associated with active enhancers (H3K4me1 and H3K27ac). Further, each of these elements undergoes bidirectional transcription to yield strand-specific eRNAs. Using ChIP-PCR and RT-qPCR, we also demonstrate that RNAP2 is recruited to the most distal Fos enhancer (here termed E1) after neuronal depolarization, and that transcription from this enhancer requires RNAP2 (Supplementary Figure S1).
To further explore this TAPE-gene relationship, we performed a KCl stimulation time-course experiment in which cultured neurons were depolarized with 25mM KCl, and RNA was isolated from neurons at multiple time points (0, 3.75, 5, 7.5, 15, 30 and 60 min) after treatment. Here, we focused on enhancer RNAs transcribed from the two most distal and most conserved Fos enhancers (upstream enhancer-1 and downstream enhancer-3) as eRNAs transcribed from enhancer-2 showed the weakest correlation with Fos mRNA in RNA-seq datasets ( Figure 3D). RT-qPCR using transcript-specific primers at the Fos locus revealed that following KCl depolarization, Fos eRNA1 and eRNA3 are significantly upregulated within 7.5 min, whereas Fos mRNA is not significantly upregulated until 15 min after stimulation ( Figure 3G). We observed similar patterns at Fosb and Nr4a1 loci, indicating that for many IEGs, eRNA induction precedes mRNA induction in response to neuronal depolarization ( Figure 3H, I). It is important to note that we used intron-spanning RT-qPCR primers for mRNA (but not eRNA; see Supplementary Data Table S4). These primers require a spliced transcript for efficient cDNA amplification, and therefore may produce a delay in mRNA response curves. However, given that mRNA splicing is co-transcriptional in neurons (66), we predict that this would only have minor effects. While we included ear-lier time points than previous studies to capture the window in which eRNA and mRNA are first induced, our data corroborate the described dynamics of eRNA transcription as well as recently identified activity-induced rapid enhancerpromoter interactions (15,33,39). Our results extend these findings and show that eRNA transcription in response to activity is rapid and follows distinct temporal profiles as compared to mRNA from linked genes ( Figure 3J).
To determine whether eRNAs are sensitive to other forms of neuronal and synaptic activation or inactivation, we treated cortical neurons with a variety of pharmacological compounds, including KCl, specific glutamate receptor agonists (AMPA and NMDA), the adenylyl cyclase activator Forskolin (FSK), or the sodium channel blocker tetrodotoxin (TTX) at 11 days in vitro (DIV). We found increased transcription of Fos eRNA1 and eRNA3 in response to KCl, AMPA, and NMDA in a dose-dependent fashion ( Figure 3K). Interestingly, while FSK treatment resulted in strong Fos mRNA induction, it did not affect Fos eRNA levels. Likewise, only mRNA levels were reduced by TTX. Together, these results suggest that Fos eRNA levels are modulated by neuronal activity states in a similar but distinct fashion compared to mRNA levels.
To gain insight into the spatial distribution of eRNAs and their response to stimulation, we performed single molecule fluorescent in situ hybridization (smFISH), a technique that allows visualization of individual eRNA and mRNA transcripts on a single-cell level. Our data confirmed a KClmediated induction of Fos eRNA1 and revealed a correlation of Fos eRNA1 and Fos mRNA transcript numbers (Supplementary Figure S2). These results suggest that eR-NAs contribute to transcriptional regulation of their target genes not only on a cell population level, but also on a single-cell level, and that enhancer transcription in single cells may explain at least part of the variability in expression from linked genes.

Functional validation of enhancer-gene pairs using CRISPR activation tools
To verify predicted enhancer-gene pairs, we sought to determine whether transcriptional activation at selected candidate enhancers was sufficient to induce mRNA at linked genes. To test this, we employed a CRISPR-dCas9 activation (CRISPRa) system in which dCas9 is fused to a strong transcriptional activator (such as VPR or VP64), enabling selective activation of targeted genomic sites ( Figure 4A-B, Supplementary Figure S3; (55,67)). We designed CRISPR single guide RNAs (sgRNAs) to target the CRISPRa system to transcriptionally active enhancer loci near the Fos, Fosb and Nr4a1 genes, as well as sgRNAs targeting proximal promoters to drive mRNA transcription directly (Figure 4C, D, Supplementary Figure S3 and S4). We chose Fos, Fosb, and Nr4a1 based on their dynamic and activitydependent nature, which makes them promising targets to study the mechanistic interactions between eRNAs and enhancer function in neurons. As a non-targeting negative control, we employed a sgRNA for lacZ, a bacterial gene that is not present in eukaryotes.
At DIV 4-5, cortical cultured neurons were transduced with separate lentiviruses expressing dCas9-VPR and      sgRNA constructs. On DIV 11, we confirmed transgene expression (indicated by mCherry reporter for sgRNA constructs and FLAG immunocytochemistry for the VPR construct; Figure 4B) and extracted RNA for RT-qPCR. At all four candidate eRNA-mRNA pairs, CRISPRa-mediated transcriptional activation of enhancers not only increased eRNA expression but also significantly induced corresponding mRNA levels ( Figure 4C). In contrast, dCas9-VPR targeting to gene promoters specifically increased target mRNA at all candidate genes but did not alter eRNA levels at three out of four candidate pairs ( Figure 4D). For example, activation of distinct enhancers either upstream (enhancer-1) or downstream (enhancer-3) of the Fos gene produced local eRNA (eRNA1 and eRNA3) induction, but also significantly increased Fos mRNA expression. Notably, transcriptional activation of the upstream enhancer did not activate the downstream enhancer, and vice versa (Supplementary Figure S3 and S4). This strongly indicates that enhancers and eRNAs can be induced by transcriptional activators, that this activation can drive mRNA expression, and that there is little crosstalk between transcriptional activation states at enhancers. Interestingly, dual Nucleic Acids Research, 2020, Vol. 48, No. 17 9561 activation of both enhancers with multiplexed sgRNAs targeting Fos enhancer-1 and enhancer-3 had additive effects and stronger mRNA induction compared to individual enhancer activation (Supplementary Figure S4). Given that enhancers can interact with promoters in enhancerpromoter loops, it is possible that transcriptional activators are close enough to act simultaneously on enhancers and promoters. However, we observed little or no effect on eRNA expression when we targeted gene promoters to drive mRNA expression ( Figure 4D, Supplementary Figure  S4), suggesting that enhancer regulation of linked mRNA is a unidirectional phenomenon. Moreover, we did not observe any effects of enhancer activation on non-targeted eR-NAs or mRNAs, supporting the site-specificity of observed CRISPRa effects (Supplementary Figure S4).
To determine whether these results translate to nonneuronal cell types that were not used to generate enhancergene pair predictions, we repeated selected experiments in C6 cells, a rat glioma dividing cell line ( Supplementary Figure S3). As in neurons, we found that recruitment of transcriptional activators (VPR or VP64) to selected Fos enhancers not only induced transcription at enhancers but also upregulated Fos mRNA. In contrast, Fos promoter targeting increased mRNA levels without altering eRNA levels. Together, these findings imply that enhancers can be activated in a site-specific manner and that observed increases in mRNA are due to enhancer activation and potentially increased eRNA levels. Further these results validate enhancer-gene pair predictions based on eRNA transcription abundance.

Enhancer RNAs are necessary and sufficient for induction of mRNA
To further interrogate the functional role of eRNAs, we explored the effect of eRNA localization on the expression of linked genes. To do so, we employed CRISPR-Display (68), a novel CRISPR approach that allowed us to tether a specific accessory RNA (acRNA) sequence to chosen target sites in the genome and investigate local effects, as compared to global over-expression approaches. Given eR-NAs from Fos enhancers showed high sequence conservation and robust effects on mRNA expression in neurons as well as in C6 cells, we designed Display acRNA sequences based on conserved regions within transcribed enhancer elements of this gene. We packaged dCas9 along with either sgRNA-eRNA (eRNA-tethering Display construct) or sgRNA-alone (no-acRNA control construct) cassettes into a single plasmid expression vector ( Figure 5A). Constructs containing either Fos enhancer-1 sgRNA or a nontargeting lacZ control were nucleofected into C6 cells, followed by RT-qPCR after a 16 h incubation period. Anchoring of an acRNA sequence based on Fos eRNA1 in close proximity to its parent enhancer (Fos enhancer-1) resulted in increased Fos mRNA levels compared to the dCasonly control ( Figure 5B). Importantly, overexpression of the eRNA1 sequence without enhancer targeting did not affect Fos mRNA expression, indicating that the effects of this eRNA are location-dependent ( Figure 5B, left). To determine whether the length of the acRNA contributes to the observed effects, we constructed eRNA-tethering CRISPR-Display plasmids with increasing acRNA lengths of 150, 300 and 450 nucleotides (nt; Figure 5A, bottom). Intriguingly, RNA length did not further increase the effect of CRISPR-Display targeting on mRNA expression (Figure 5B, right), suggesting that eRNA-mediated increases in mRNA expression are not directly proportional to the size of eRNA delivered. Since the effects of eRNA targeting appeared to be location-dependent, we next sought to determine their location specificity. We targeted enhancer-1, -3, and a non-regulatory control region between enhancer-1 and the promoter with Display constructs tethering 150 nt long sequences of eRNA1, 3, or a control RNA based on the targeted control region ( Figure 5A, bottom). These experiments revealed that only Fos eRNA1 tethered to its own enhancer induced mRNA ( Figure 5C, middle). Importantly, no significant effects were observed when other RNAs were tethered to enhancer-1 ( Figure 5C, middle), nor when eRNA1 was tethered to either enhancer-3 or the control locus ( Figure 5C, left and right). These results suggest that the observed effects are specific to eRNA1 when localized to its origin enhancer. Interestingly, eRNA3 did not produce the same effects on mRNA expression as eRNA1 when targeted to its corresponding enhancer-3 ( Figure 5C, right). This could either be due to technical differences between the acRNA designs or hint towards different functional roles between the two eRNAs. While none of the acRNAs contain specific binding motifs that the authors are aware of, it is possible that the chosen 150 nt of eRNA3 does not carry the required properties to be functional in this assay and that a different region of the eRNA3 sequence would show similar results to the eRNA1 acRNA. This discrepancy could also reflect differences in chromatin looping at baseline between the two enhancers. A recent study by Beagan et al. assessed activity-dependent chromatin restructuring and enhancerpromoter loops and found activity-responsive looping between the Fos promoter and enhancer-1 but not enhancer-3 (39). Existing enhancer-promoter interaction could potentially be required for acRNA mediated effects. More importantly, these experiments provide novel evidence that Fos eRNA1 acts locally and is sufficient to induce the Fos gene.
Based on the CRISPR-Display results, we next sought to address the functional requirement of eRNAs in activitydependent gene transcription in cortical neurons. We employed an anti-sense oligonucleotide (ASO) strategy to directly target eRNA1 while leaving mRNA and other enhancer functions unperturbed. Rat primary cortical cultures were treated with sequence-specific eRNA1 ASOs for 24 h prior to RNA harvesting followed by RT-qPCR (Figure 6A). ASOs targeted to Fos eRNA1 induced a robust decrease in eRNA1 expression but did not alter expression of eRNAs from other Fos gene enhancers (reinforcing the concept of functional independence of Fos eRNAs). Notably, Fos eRNA1 ASOs also produced a significant decrease in Fos mRNA levels, both at baseline and following neuronal depolarization with KCl ( Figure 6A, C). These results suggest that Fos eRNA1 is not only required for normal expression from the Fos gene, but also for neuronal activitydependent expression of this IEG. In contrast, we found that knockdown of Fos mRNA with an ASO targeted to the mRNA had no effect on eRNA synthesis from any en- hancer, further supporting a unidirectional model of eRNA function ( Figure 6B). Overall, these findings demonstrate that altering the levels of a single eRNA is sufficient to modulate gene expression.

Enhancer RNAs interact with but do not depend on chromatin remodelers
The mechanisms by which eRNAs can regulate proximal mRNA transcription remain poorly understood, and while there is evidence for specific eRNA-protein interactions in the literature, no general mechanism has been identified. eRNAs have been demonstrated to be involved in the release of transcriptional repressors, enhancer looping, and epigenetic modifications (15,(69)(70)(71)(72). One possible role of eRNAs lies in the regulation of dynamic chromatin reorganization. Enhancer activation is often accompanied by the recruitment of CBP, CREB, MEF2, NPAS4 and FOS proteins to enhancers near activity-regulated genes (e.g. Fos, Rgs2 and Nr4a2; (13)). Furthermore, enhancer manipulations have been shown to induce changes in H3K27ac both directly and indirectly. Chen et al. demonstrated increased H3K27ac at Fos enhancers after dCas9-p300 recruitment to these enhancers (10)). Similarly, other studies have reported decreased histone acetylation at enhancers not only when using histone deacetylase tethering constructs, but also other repressors such as KRAB and histone demethylases (76). To interrogate the relationship between eRNAs and these enhancer-binding epigenetic modifiers, we fo-cused CBP, which is recruited to enhancers upon activation and has been shown to bind eRNAs (69). First, to determine whether Fos eRNA and mRNA expression is regulated by CBP levels, we designed a custom shRNA expression vector targeting Crebbp mRNA (which codes for CBP protein). Specific shRNA sequences targeting Crebbp mRNA or a scrambled control sequence were cloned into a lentivirus-compatible expression vector and used for lentiviral packaging. At DIV 4-5, cortical cultured neurons were transduced with either control or Crebbp shRNA. On DIV 11, we confirmed transgene expression (indicated by mCherry reporter) and extracted RNA for RT-qPCR. Intriguingly, Crebbp knockdown reduced Fos mRNA by ∼40%. However, we observed no reduction in eRNA1 or eRNA3 levels ( Figure 7A). Even though CBP has been shown to interact with enhancers, our data indicate that Crebbp is crucial for Fos mRNA induction but is not required for eRNA transcription.
We next tested whether CREB-CBP interactions were critical for Fos eRNA and mRNA expression using a selective small molecule compound (666-15) that blocks CREBmediated gene transcription. DIV 11 cortical neurons were treated with 666-15 and KCl for 1hr followed by RNA harvesting and RT-qPCR. CREB inhibition significantly blunted the Fos mRNA response to KCl ( Figure 7B). However, CREB inhibition did not alter Fos eRNA1 or eRNA3 expression either at baseline or in response to KCl, demonstrating that eRNA expression does not require CREB-CBP interactions.
Nucleic Acids Research, 2020, Vol. 48, No. 17 9563 Fos mRNA ASO Since CBP is not required for Fos eRNA transcription, we next investigated the role of CBP in enhancer function. First, we tested the effects of CBP's core histone acetyltransferase (HAT) domain on enhancer function. We paired a sgRNA specific to the Fos enhancer-1 locus with a dCas9 fusion containing the core HAT domain from p300, which is nearly identical to the HAT domain from CBP (57). HAT targeting to either Fos enhancer in C6 cells elevated Fos mRNA to a similar extent as CRISPRa-based activation and CRISPR-based eRNA1 targeting ( Figures 5B, C, 7C). While eRNA1 and eRNA3 levels in the same experiment did appear to be elevated following dCas9-HAT targeting to their respective enhancers, these changes were more variable and did not meet statistical significance as compared to lacZ controls (Supplementary Figure S5). Overall, these data suggest that eRNAs may link enhancer transcription to downstream chromatin remodeling via histone acetylation.
Recent work suggested that eRNAs interact with CBP and modulate the HAT activity of this protein (69). However, eRNAs have also been shown to interact with proteins containing a bromodomain (which binds acetylation marks on histones) (72). Intriguingly, CBP (and the dCas-HAT fusion used in Figure 7C) contains both a bromodomain and a HAT domain ( Figure 7D), implying either function is possible. Furthermore, enhancer acetylation has recently been shown to be a crucial regulatory mechanism of activity-induced dynamic gene expression and transcriptional bursting (10). To investigate whether eRNA binding via either of these mechanisms contributes to eRNA function at Fos enhancers, we performed RNA electrophoretic mobility shift assays (REMSAs) using ∼150nt RNA probe sequence based on Fos eRNA1, eRNA3 and a control RNA (based off the same control locus as CRISPR-Display control RNA, Figure 5A, C). REMSA probes were transcribed in vitro with fluorescein-labeled bases, and subsequently incubated with recombinant CBP HAT or bromodomains ( Figure 7E-I). To ensure that eRNA-protein interactions were specific and not due to general RNA binding of CBP domains, all REMSA experiments used buffer containing 0.1mg/ml of yeast tRNA. Strikingly, we observed almost complete binding of Fos eRNA1 and the CBP HAT domain at the highest protein concentration (0.6 M, Figure 7E). While both tested eRNAs bound the CBP HAT domain, it is noteworthy that all tested RNAs, including the control RNA, interacted with the HAT domain. In line with our CRISPR-Display results, these findings suggest that the observed interactions between eRNAs or control RNA and the CBP HAT domain are not entirely sequence-specific, but likely regulated by their short half-life and local transcription. As eRNA1 showed almost complete binding and the strongest effects on mRNA expression in our Display experiments, we used this sequence for control experiments. These experiments demonstrate that the eRNA1-HAT interaction was competitively inhibited by unlabeled eRNA1, suggesting that binding was not mediated by the fluorescent label. Additionally, there was no appreciable binding of eRNA1 to the CBP bromodomain at identical concentra-G CBP -Bromo   tions ( Figure 7H-I). This result thus supports recent findings by Bose et al. that various eRNAs can interact with CBP and potentially increase its HAT activity (69). Together with our eRNA-tethering results, this data extends the previous findings, suggesting that functional specificity of eRNAs is driven by the location of eRNA transcription rather than their sequence.
Taken together, these results suggest that while eRNAs are not dependent on CREB and CBP, they can facilitate transcriptional induction through direct interaction with the CBP HAT domain. Our findings suggest a model for eRNA function in which eRNAs participate in enhancerpromoter communication where they can interact with epigenetic modifiers such as CBP ( Figure 7J).

DISCUSSION
Distal enhancer elements in DNA enable higher-order chromatin interactions that facilitate gene expression programs and thus contribute to cellular phenotype and function (1-3). In the developing brain, the majority of enhancer elements exhibit temporally specific emergence during precise developmental windows, with only ∼15% of enhancers being utilized continually from late embryonic development into adulthood (5,7). These developmentally regulated enhancers contribute to cell-and tissue-specific gene expression patterns that establish communication within and between brain structures (7,73,74). Not surprisingly, enhancers utilized in early embryonic brain development possess the highest degree of sequence conservation across species, suggesting that robust evolutionary pressures drive enhancer function (7). In the postnatal and mature brain, enhancers continue to play a widespread role in the activity-dependent transcriptional programs that regulate key aspects of neuronal plasticity and function (5,12,13,(16)(17)(18)27,41,75). Repression or deletion of enhancer elements has profound effects on the genes that they control, including complete inactivation (12,16,41,76). Likewise, targeted enhancer activation induces robust upregulation of linked genes, suggesting that enhancers serve as bidirectional regulators of gene activity (57,73).
While genome-wide enhancer identification traditionally relies on extensive ChIP-seq, ATAC-seq, and/or HiC-seq datasets to asses location and activity states, our findings are in line with recent work by Wang et al. and emphasize the advantages of transcriptional information as an indicator of local chromatin states (77). In this study, the authors developed a machine learning tool to predict chromatin landscapes with nucleosome resolution based exclusively on nascent transcription. Along with this line of work, others have underlined the tight relationship between enhancer transcription and transcription factor activity (78), as well as enhancer and promoter function (79). Our results demonstrate that eRNA levels can be a useful measure of enhancer activity state, and often represent a distinct readout from that provided by ATAC-seq signal at the same enhancer. For example, whereas we detected many eR-NAs induced by neuronal depolarization, we observed that ATAC-seq signals at these activated enhancers (or the promoters of linked genes) decreased following depolarization. While the molecular origins of this effect are unclear, this finding may suggest that at the timescales examined here, ATAC-seq quantification cannot discriminate between active and poised enhancer states for constitutively poised activity-regulated enhancers.
Our results also extend recent efforts to identify and validate enhancer-gene pairs on a large scale using CRISPR interference (CRISPRi). One study perturbed enhancer function using CRISPRi sgRNA libraries followed by singlecell sequencing to assess enhancer-gene links. The resulting dataset demonstrated that at least 33% of the investigated 770 enhancer-gene pairs skipped at least one proximal TSS (80). Similarly, another recent study by Fulco. et al. developed CRISPR-FlowFISH, an approach that uses CRISPRi to perturb hundreds of regulatory elements in tandem with single-cell analysis of mRNA abundance (36). The results from this study suggest that genes can be regulated by multiple distal enhancers and that some enhancers regulate more distal genes while skipping over their more proximal genes. Furthermore, this work demonstrates poor accuracy of several commonly used enhancer-gene pair prediction methods (e.g. proximity, Hi-C data and chromatin marks) when used individually. To solve this problem the authors developed a new activity-by-contact (ABC) model in which enhancer-gene predictions are based on the activity state of the enhancer (measured by DNase-seq and H3K27ac ChIPseq), the frequency of 3D enhancer-promoter interactions (measured by Hi-C), and the gene expression response to CRISPRi targeting. In contrast to this study our method capitalizes on inherent transcriptional variability between samples and therefore does not require large-scale enhancer perturbations to predict enhancer-gene pairs. Furthermore, our results demonstrate that it is possible to resolve activity states of regulatory elements and their corresponding target genes with only a few datasets and suggests that eRNA induction levels can predict mRNA response to stimulation more reliably than changes in chromatin accessibility. Recent work by Beagan et al. examined three-dimensional chromatin structure changes in response to neuronal activity, and demonstrated the formation of enhancer-promoter loops within 20 min of stimulation at IEGs including the Fos gene (39). Consistent with our findings, this report also showed increased eRNA and enhancer-promoter interaction changes prior to induction of mRNA from these loci. This study highlights the need to correlate activity states at promoters and enhancers to generate higher confidence enhancer-promoter link predictions, and shows that loop strength translates directly to transcriptional output. Building on these findings, our data emphasize the importance of incorporating transcriptional information in enhancer identification, as eRNA levels may change more rapidly than overlapping chromatin status and potentially even three-dimensional structures. Using measures of eRNA abundance may therefore provide higher temporal resolution when investigating stimulus-dependent changes in enhancer landscape.
Although it is well accepted that genomic enhancers play critical roles in tuning the spatiotemporal nature of transcription from linked genes, techniques typically used to examine enhancer function (e.g. enhancer deletion (81), Cas9based mutation (59,82), or activation/inactivation with dCas9 fusion proteins (12,57,83-86) interfere with both the genomic locus and eRNAs transcribed from that locus. Therefore, these approaches cannot dissociate the effects of enhancer function and eRNA function. To address this problem, we first implemented genome-wide transcriptional profiling in our enhancer identification pipeline to gain a better understanding of the relationship between enhancer transcription, enhancer activity, and downstream transcriptional regulation. Moreover, we took two different approaches that directly target eRNAs in order to examine their function separately from enhancer function. First, we used a novel CRISPR-Display approach to target Fos eRNAs to their own enhancer. These results demonstrate that Fos eRNA1 is sufficient to induce Fos mRNA and provide novel evidence for a location-dependent functional role of eRNA ( Figure 5). Furthermore, our data shows that while Fos eRNA expression is independent of CREB and CBP function, eRNAs can interact with CBP through direct binding to the HAT domain. Likewise, CRISPR-dCas9 mediated recruitment of a HAT domain recapitulated the effects of eRNA-tethering at the same Fos enhancer loci (Figure 7). Secondly, we employed stable, cell-penetrating ASOs to target eRNA for degradation. These results suggest that eRNA is necessary for normal expression of Fos mRNA, both under basal conditions and after neuronal depolarization. (Figure 6).
Overall, these results agree with a previous report demonstrating that eRNAs transcribed from activity-dependent enhancers are necessary for induction of mRNA from linked genes (15). This report utilized lentiviral shRNA knockdown approaches to directly target activity-induced eRNAs near Arc and Gadd45b genes and followed this knockdown with KCl depolarization to induce mRNAs. Targeted shRNA knockdown of eRNA specifically blocked mRNA induction at these genes but not other IEGs induced by neuronal activation (Fos, Egr1). Our results extend these important findings in two ways. First, given that the Fos gene exhibits multiple enhancers and activity-dependent eRNAs, we were able to address the functional relationship between eRNAs near the same gene. Our results suggest that while eRNAs do regulate mRNA induction at linked genes, eRNAs are functionally independent of each other. Thus, ASO-mediated knockdown of eRNAs transcribed from the most distal Fos enhancer did not downregulate eRNAs transcribed from other enhancers ( Figure 6). Secondly, in parallel experiments we were able to target Fos mRNA for knockdown using an identical approach. These results demonstrate that the relationship between eRNA and mRNA levels at the same gene is unidirectional, i.e. that mRNA knockdown does not also reduce eRNA levels. This is a critical control at autoregulating IEGs like Fos, given that the protein product of this gene is a transcription factor that localizes to enhancers in an AP1 complex with Jun family members (14).
Biological roles of lncRNAs are generally linked to their ability to bind functionally active proteins to operate as molecular guides, decoy molecules, scaffolding, or even allosteric modulators (32,87). In agreement with this concept, a large number of chromatin-associated proteins bind RNA in addition to DNA (43,88,89), and several well-characterized transcriptional regulators have recently been shown to possess functional interactions with eR-NAs (2,15,(69)(70)(71)90,91). For example, eRNAs have been shown to bind the ubiquitous transcription factor Yin-Yang 1 (YY1) to 'trap' YY1 at the enhancer, thus facilitating its action at local YY1 motifs in DNA (91). In this study, a similar CRISPR-dCas9 system was used to tether Arid1a RNA in close proximity to an enhancer YY1-binding motif. Our CRISPR-Display experiments build on this work by showing direct changes in target gene expression that are dependent on the target location but not sequence length. Our data confirms and highlights the importance and sufficiency of eRNAs as transcriptional organizers.
Similarly, eRNAs can act as decoy molecules for negative elongation factor (NELF) complexes, which are important regulators of RNAP2 pausing and transcriptional bursting (15). In line with previous findings that eRNAs interact with CBP and stimulate its activity as a HAT at enhancer loci (69), we found that Fos eRNAs can interact with CBP through direct binding to the HAT domain and that eRNA or HAT recruitment to the enhancer increases mRNA expression. Additionally, our data suggest that the specificity of this interaction is based on the location of eRNA transcription, likely in combination with their short half-life rather than sequence specificity.
While our results do not rule out other eRNA-protein interactions, our findings are consistent with recent observations that histone acetylation plays a key role at enhancers by influencing transcriptional properties of corresponding genes (10). This study demonstrated that histone acetylation at enhancers increases transcriptional bursting at linked genes, and that these effects are mediated by BRD4, a bromodomain-containing protein important for RNAP2 phosphorylation and transcription. Our results provide the first evidence that eRNA function is dependent on eRNA location and partially dependent on sequence, but not sequence length. Future studies will be required to explore how different factors, such as distance from their target gene, influence enhancer and eRNA function and characterize other eRNA/protein interactions in more detail. This will help to determine whether common regulatory mechanisms dictate expression of different eRNAs targeting the same gene, as we observed some characteristic and functional differences between Fos eRNA1 and eRNA3. It will further help to unravel whether a group of eRNAs that regulate the same gene have distinct functional mechanisms. Finally, it will be crucial to understand the interplay of different enhancers and eRNAs and how they orchestrate gene expression programs and potentially fine-tune responses to specific stimuli.
The vast majority of gene variants linked to human health and disease by genome-wide association studies are located in non-coding regions of the genome (25,27,92), with putative enhancers containing more disease-linked single-nucleotide polymorphisms than all other genetic loci combined (93). Disease-linked genetic variants could affect enhancer activity either via direct modification of enhancer DNA sequence (e.g. disruption of a transcription factor motif) or by alterations in long-range chromatin interactions between enhancers and gene promoters (reviewed in (94)). Indeed, numerous diseases have already been linked to sequence variations in enhancer regions (92,(95)(96)(97), including complex polygenic conditions such as depression (23,98), obesity (28,98), schizophrenia (22,99), bipolar disorder (22), Alzheimer's disease (38,100), and autism spectrum disorders (24,29). This growing link between enhancer activity and brain function strongly highlights the need to better understand the mechanistic interactions that regulate enhancer function at the molecular level, and also suggests that enhancers could be attractive targets for a new generation of disease therapeutics.

DATA AVAILABILITY
All relevant data that support the findings of this study are available on day-lab.org/resources or by request from the corresponding author (J.J.D.). The CRISPR-Display construct will be made available in the Addgene plasmid repository (Addgene #159086). Bash scripts and R code for TAPE analysis will be made available via GitLab (https://gitlab.rc. uab.edu/day-lab).
Sequencing data that support the findings of this study are available in Gene Expression Omnibus (GSE150499 and GSE150589).