Altered enhancer transcription underlies Huntington’s disease striatal transcriptional signature

Epigenetic and transcriptional alterations are both implicated in Huntington’s disease (HD), a progressive neurodegenerative disease resulting in degeneration of striatal neurons in the brain. However, how impaired epigenetic regulation leads to transcriptional dysregulation in HD is unclear. Here, we investigated enhancer RNAs (eRNAs), a class of long non-coding RNAs transcribed from active enhancers. We found that eRNAs are expressed from many enhancers of mouse striatum and showed that a subset of those eRNAs are deregulated in HD vs control mouse striatum. Enhancer regions producing eRNAs decreased in HD mouse striatum were associated with genes involved in striatal neuron identity. Consistently, they were enriched in striatal super-enhancers. Moreover, decreased eRNA expression in HD mouse striatum correlated with down-regulation of associated genes. Additionally, a significant number of RNA Polymerase II (RNAPII) binding sites were lost within enhancers associated with decreased eRNAs in HD vs control mouse striatum. Together, this indicates that loss of RNAPII at HD mouse enhancers contributes to reduced transcription of eRNAs, resulting in down-regulation of target genes. Thus, our data support the view that eRNA dysregulation in HD striatum is a key mechanism leading to altered transcription of striatal neuron identity genes, through reduced recruitment of RNAPII at super-enhancers.


Results
Differential expression of eRNAs in the striatum of HD R6/1 mice. To explore the hypothesis that alteration of eRNA transcription might be a component of HD pathogenesis, we assessed eRNAs at genome-wide scale, analyzing non-coding RNAs from RNA sequencing (RNAseq) data previously generated in the striatum of control and HD R6/1 transgenic mice 2 . A strand-specific total RNA sequencing protocol was used to generate sequencing reads, thereby allowing analysis of long non-coding RNAs (see Methods). To identify eRNAs, i.e. non-coding RNAs synthesized from enhancers, we first excluded signals within genic regions, defined as the interval starting 3 kb upstream of the transcription start site and ending 10 kb downstream of the transcription termination site, since they might result from polymerase read-through of genic transcripts (ref. 19 and Methods). Second, we filtered RNA signals resulting from enhancer regions using H3K27ac ChIP-seq data, generated from the striatum of control (WT) and HD R6/1 mice 2 . As a result, a total of 6068 eRNAs were selected based on H3K27ac occupancy (Fig. 1A). Analysis of differentially expressed eRNAs between R6/1 and WT striata was performed, as well as eRNA annotation, which provided gene-eRNA associations (see Methods). 677 and 335 eRNAs were found decreased and increased, respectively (Fig. 1B and Table S1). Noticeably, down-regulated eRNAs were globally expressed at higher levels than up-regulated eRNAs in WT mice (Fig. 1C). Decreased expression of selective eRNAs, including eRNAs in neighborhood of Rgs4, Rgs9, Slc24a4, Chn1, Gpr6, Ajap1, Bcr and Asphd2 genes, was confirmed by q-RT-PCR (Fig. 1D,E and S1A). Expression of Hps1-associated eRNA, which was unchanged between R6/1 and WT, was used as a negative control (Fig. 1D,E and S1A). mRNAs transcribed from Rgs4, Rgs9, Slc24a4, Chn1, Gpr6, Ajap1, Bcr and Asphd2 were also decreased, in contrast to Hps1 mRNA (Fig. S1B), suggesting a link between eRNA and mRNA deregulation in R6/1 mouse striatum. To evaluate the degree of conservation of the mechanism, we analyzed another HD mouse model widely used in the field, the Q140 knockin model, expressing full-length mutant Htt. These mice display progressive transcriptional dysregulation, particularly in the striatum 10 . Rgs4, Rgs9, Slc24a4, Chn1, Gpr6, Ajap1, Bcr and Asphd2 mRNAs were also decreased in the striatum of 12 month-old Q140 mice (ref. 10 and S1B). eRNAs associated with these genes were significantly decreased or showed a tendency to the decrease, except Slc24a4-associated eRNA (Fig. S1C), suggesting that eRNA dysregulation in HD striatum is a general mechanism.
Target genes of decreased eRNAs in R6/1 striatum are enriched in neuronal function genes.
We investigated whether down-regulated eRNAs in R6/1 striatum were enriched in genes implicated in specific functions. Gene ontology analysis (GO) using GREAT 22 showed that enhancer regions involved in decreased eRNAs in R6/1 striatum were strongly associated with genes enriched in biological processes linked to neuronal activity, including neuronal transmission, synaptic plasticity and learning and memory ( Fig. 2A). In contrast, regions involved in increased eRNAs were close to genes enriched in biological processes related to stem cell proliferation ( Fig. 2A). Thus, down-and up-regulated eRNAs in R6/1 striatum associate with genes that display neuronal and developmental signatures, respectively.

Target genes of decreased eRNAs in R6/1 striatum are enriched in down-regulated genes.
Down-regulated genes in R6/1 striatum also present a strong neuronal signature 2 . Since eRNAs positively regulate their target genes, this suggests that decreased eRNAs might modulate expression of genes down-regulated in R6/1 striatum. Integrated analysis showed that target genes of enhancers associated with decreased eRNAs in R6/1 striatum were enriched in down-regulated genes (Fig. 2B). Moreover, levels of eRNAs in the neighborhood of decreased mRNA in R6/1 striatum were globally reduced (Fig. 2C). As expected, the subset of down-regulated genes associated with decreased eRNAs in R6/1 striatum displayed a clear neuronal signature (Fig. 2D). Interestingly, they were enriched in genes controlling neuronal excitability, including genes coding for voltage-gated potassium channels such as Knca4, Kcnab1 and Kcnj4 (Fig. 2E). In contrast, target genes of increased eRNAs were not enriched in up-regulated genes in R6/1 striatum (Fig. 2B). However, eRNAs associated with up-regulated genes in R6/1 vs WT striatum were globally increased (Fig. 2C), and these genes included developmental genes such as Onecut1 and Onecut2, expressed in neural stem cells ( Fig. 2E and ref. 23). Together, these results suggest that eRNA down-regulation has a broader influence on gene expression than eRNA up-regulation in R6/1 striatum, with biological impact of decreased eRNAs in R6/1 striatum affecting neuronal activity, including neuronal excitability, and that of increased eRNAs influencing neuronal fate. These two effects of eRNA dysregulation might synergistically contribute to loss of neuronal differentiated state of HD striatal neurons.    Altered transcribed enhancers in R6/1 striatum are enriched in super-enhancers. Target genes of super-enhancers, a category of broad enhancers regulating cell type-specific identity genes, are preferentially down-regulated in R6/1 striatum 2 . We therefore hypothesized that decreased eRNAs in R6/1 striatum were transcribed from super-enhancers. Integrated analysis of eRNA and H3K27ac ChIP-seq data on R6/1 and WT striatum supported this hypothesis, since H3K27ac-enriched regions associated with decreased eRNAs were broader in comparison with those associated with increased eRNAs (Fig. 3A). Consistently, enhancers associated with down-regulated eRNAs in R6/1 striatum were enriched in super-enhancers, in contrast to up-regulated eRNAs (Fig. 3B). In addition, decreased eRNAs produced from super-enhancers displayed a strong neuronal signature, whereas those produced from typical enhancers did not present any specific functional signature (Fig. 3C). This shows that neuronal signature of down-regulated eRNAs is essentially contributed by transcribed super-enhancers. We also crossed the list of decreased eRNA-associated genes regulated by a super-enhancer with the list of decreased eRNA-associated genes down-regulated in R6/1 striatum. Both lists were largely overlapping (Fig. 3D), showing that most down-regulated target genes associating with reduced eRNAs in R6/1 vs WT striatum are under the control of a super-enhancer. Genes within the resulting sub-list (36 genes) contained striatal markers (e.g. Drd1, Rgs9), including striatal voltage-gated potassium channels (e.g. Kcnj4/Kir2.3, ref. 24) (Fig. 3D). Finally, the size of H3K27ac enriched regions at enhancers associated with decreased eRNAs in R6/1 striatum were reduced in R6/1 when compared to WT mice (Fig. 3A), suggesting a link between reduction of H3K27ac-enriched regions and eRNA transcription in R6/1 striatum. Together, these data indicate that altered super-enhancer transcription in HD mouse striatum contributes to down-regulation of striatal identity genes.   RNAPII binding sites are lost at enhancers associated with decreased eRNAs in R6/1 striatum.
RNAPII is enriched at start sites of transcribed enhancers 16,19,25 . To investigate whether RNAPII signal varied between WT and R6/1 samples, we integrated RNAPII ChIPseq data previously generated 2 to the analysis, filtering eRNAs containing RNAPII peaks. Out of 6068 eRNAs, 3248 were retrieved in WT mice following this filtering. A similar number (2990) were identified from R6/1 mice (Fig. 4A). Thus, a substantial proportion (≈ 1/2) of eRNAs were associated with a RNAPII peak and the numbers of eRNAs with RNAPII peaks were similar between WT and R6/1 samples, though a slight decrease (8%) was observed in R6/1 vs WT striatum. 331 eRNAs with RNAPII peaks were decreased in R6/1 striatum when compared to WT. Out of these 331 eRNAs, 127 lost RNAPII peaks in the R6/1 condition, which corresponded to a 38% loss. Thus, RNAPII peaks were dramatically reduced at R6/1 enhancers associated with decreased eRNAs (Fig. 4A), suggesting that loss of RNAPII contributes to decreased transcription of eRNA in HD mouse striatum. In contrast, 351 eRNAs were found increased in R6/1 vs WT striatum, of which only 18 gained RNAPII peaks in R6/1 mice (e.g. 5%; Fig. 4A), suggesting that RNAPII gain may not be a major mechanism governing eRNA upregulation in R6/1 striatum.
Enhancers associated with deregulated eRNAs in R6/1 striatum are enriched in selective DNA motifs. We then asked whether deregulated eRNAs in R6/1 striatum were enriched in binding sites for transcriptional regulators. Remarkably, transcribed enhancer sequences, whether they led to deregulated eRNAs or not, were enriched in GC content (Fig. 4B). This was particularly true for up-regulated eRNAs (Fig. 4B), which were also enriched in GC-rich DNA motifs ( Fig. 4C and S1D). Within the list of DNA binding sites enriched in up-regulated eRNAs, none were recognized by transcriptional activators substantially up-regulated or transcriptional repressors down-regulated in the striatum of R6/1 mice 2 . However, Klf5 and Usf1 were slightly but significantly up-regulated in R6/1 striatum (adjusted p-value is 0.03 in both cases, ref. 2). Usf1 acts as a chromatin insulator element, promoting gene activation through inhibition of polycomb complex activity 26 , whereas Klf5 is a developmental transcription factor involved in maintenance of embryonic stem cells undifferentiated state 27,28 .  These DNA motif signatures are consistent with the functional signature of up-regulated eRNAs in R6/1 striatum, characterized by biological processes linked to both gene silencing and regulation of stem cell proliferation ( Fig. 2A).

Figure 4. Enhancers leading to decreased eRNAs in R6/1 vs WT striatum present reduced RNAPII signals in R6/1 striatum and are enriched in SRF binding sites. (A) Graphs showing the numbers of RNAPII peaks in
In contrast to increased eRNAs in R6/1 striatum, decreased eRNAs were exclusively enriched in a DNA motif recognized by SRF (Fig. 4D), a transcription factor that recruits RNAPII to transcription sites. Additionally, Srf expression, which is substantial in mouse striatum, is significantly decreased in HD mouse striatum (Fig. 4E). Since SRF is a key player of neuronal plasticity 29 , reduced Srf expression in HD striatum might contribute to decreased transcription at striatal enhancers, through impaired recruitment of RNAPII, which might in turn lead to down-regulation of striatal enhancer target genes.

Discussion
In this study, we have investigated enhancer transcription in normal and Huntington's disease mouse striatum. We show that eRNAs are widely expressed from striatal enhancers and many eRNAs are decreased in R6/1 vs WT mouse striatum. These decreased eRNAs in R6/1 striatum display a strong neuronal signature: their target genes are highly enriched in biological processes linked to synaptic plasticity, neuronal transmission and learning. This result is consistent with the fact that decreased eRNAs in R6/1 striatum are preferentially expressed from super-enhancers, which control genes that define cell type-specific identity and function. Our data also show that transcribed enhancers resulting in decreased eRNAs in R6/1 striatum are selectively enriched in a DNA motif recognized by SRF, a transcription factor altered in R6/1 striatum, display loss of RNAPII signals, and associate with genes that are down-regulated in R6/1 striatum. Together, these data indicate that repression of neuronal genes in HD striatum involves an epigenetic mechanism, eRNA dysregulation, originating from altered recruitment of RNAPII at super-enhancers, possibly due to reduced SRF levels.
Down-regulated genes in HD mouse striatum are preferentially regulated by super-enhancers 2 . As a result, they display a strong neuronal signature 2,4,9,10 . In this study, we showed that transcription at substantial numbers of super-enhancers was decreased in R6/1 striatum (Fig. 3). Target genes of super-enhancers showing decreased eRNAs were down-regulated in R6/1 striatum (Fig. 2). Since enhancer transcription is a marker of enhancer activity 16 and eRNAs positively regulate transcription of their target genes 18 , our results uncover that altered super-enhancer activity broadly contributes to down-regulation of markers of striatal identity in HD mice.
In contrast, increased (super-) enhancer activity had no broad impact on gene up-regulation in HD striatum, since up-regulated eRNAs in R6/1 striatum poorly overlapped with super-enhancers and their target genes were not enriched in up-regulated genes (Fig. 2B). However, our results showed that up-regulation of the developmental genes Onecut1 and Onecut2, expressed in neural stem cells 23 , correlated with augmented transcription of associated eRNAs (Fig. 2E). This suggests that selected increased eRNAs contribute to striatal neuron dedifferentiation. Thus, increased and decreased transcription at enhancers might both participate to loss of differentiated state of striatal neurons in response to the HD mutation, through re-activation of genes expressed in immature neurons and repression of genes defining striatal neuron identity, respectively. Using a cellular model of inflammation, Hah et al. reported extensive transcription within super-enhancers 19 . They also showed that super-enhancer transcription was required for the regulation of target genes involved in innate immunity. Our results showed that neuronal super-enhancer transcription was also extensive, which suggests that super-enhancer transcription is a general mechanism critical to the regulation of genes that determine cell type-specific identity and function.
We analyzed eRNA transcription in basal conditions -as opposed to stimulus-induced conditions -and found that many eRNAs (6068) were transcribed from mouse striatal enhancers. This result is consistent with previous study, which identified 8990 and 7779 eRNAs in mouse cortex and cerebellum, respectively 30 . Thus, many enhancers are constitutively active in a mouse brain tissue and our results indicate that repression of neuronal genes in HD striatum predominantly results from altered activity of constitutively active enhancers.
Using primary cortical cultures, Kim et al. previously identified neuronal stimulus-responsive enhancers using a genome-wide approach 16 . The epigenetic signature of these activity-regulated enhancers was characterized by increased H3K27ac signals and augmented eRNA transcription 16,17 . Our results showing that reduced eRNAs transcription correlates with decreased H3K27ac signals in R6/1 striatum (Fig. 3A) further support the idea of an interaction between mechanisms regulating enhancer transcription and H3K27 acetylation.
The transcription factor FOS was enriched at neuronal stimulus-responsive enhancers 17 . However, it is not likely to be the case for constitutively active enhancers, since basal expression of Fos is low in neurons. Remarkably, enhancers that associated with reduced eRNA transcription in HD mouse striatum were enriched in a DNA motif recognized by another transcription factor, SRF (Fig. 4D). SRF is a key regulator of synaptic plasticity and this function of SRF includes mechanisms that modulate constitutively expressed genes 31 . Particularly, basal gene expression changes resulting from inactivation of the Srf gene in the adult forebrain accounted for altered synaptic plasticity, specifically long term depression (LTD) 31 . Noticeably, it has been reported that synaptic plasticity, including LTD, is impaired in HD mice 32 . Moreover, Srf transcription is reduced in R6/1 striatum (Fig. 4E). Whether altered SRF regulation in HD mouse striatum affects constitutively active enhancers, including super-enhancers, thereby contributing to altered expression of synaptic plasticity genes, is an intriguing hypothesis. This possibility would be consistent with functional enrichment analysis showing that down-regulated eRNAs in R6/1 striatum were enriched in GO terms such as "regulation of synaptic plasticity", "long term depression" (LTD) and "learning" (Fig. 2A).
Voltage-gated potassium channels were enriched in the subsets of down-regulated genes associated with reduced enhancer transcription (Fig. 2D). Specific regulation of potassium currents is required to unique electrophysiological properties of striatal neurons, including maintenance of the hyperpolarized state 24,33,34 . Decreased expression of voltage-gated potassium channels in HD mice, which results in altered excitability properties of striatal neurons, due to depolarized resting state membrane potentials, was suggested to partially account for preferential vulnerability of striatal medium spiny neurons in HD 24 . Our data showing that altered transcribed enhancers in HD mouse striatum display a voltage-gated potassium channel signature might indicate that eRNA dysregulation underlies increased striatal vulnerability in HD.
We showed that transcribed enhancers display high GC content (Fig. 4B), which suggests that enhancers are subject to regulation by DNA methylation. Strikingly, decreased eRNAs in R6/1 vs WT striatum were less GC-rich, in comparison with up-regulated eRNAs or with unchanged eRNAs. This is consistent with results showing that super-enhancers are hypo-methylated, which might facilitate a chromatin state permissive to transcription 35,36 .
In conclusion, we propose that targeting striatal enhancers in an attempt to improve enhancer transcription might prevent repression of neuronal genes in HD. The development of therapeutic strategies permitting to target enhancers may thus represent a future challenge.

Materials and Methods
Animals. Hemizygous R6/1 (≈ 150 CAGs) 37  RNA extraction and qRT-PCR. Each sample was prepared from single striatum of R6/1 and WT littermate mice (for R6/1 vs WT comparisons) or from single striatum of Q140 and WT littermate mice (for Q140 vs WT comparisons). Biological replicates (4 to 6) were performed for each group. Tissues were finely cut with a razor blade and total RNA was extracted using TRIzol reagent (Invitrogen). An additional DNAse treatment (Euromedex) was added before RNA purification by phenol/chloroform extraction and ethanol precipitation. cDNA synthesis was performed either on 0.5 μ g of total RNA (iScript Reverse transcription Supermix for RT-qPCR kit; Bio-Rad) or on 2 μ g of total RNA for strand-specific reverse transcription (SuperScript II Reverse Transcriptase; Invitrogen or High-Capacity cDNA Reverse Transcription Kit; Appliedbiosystems). Gene-specific primers are available upon request. qRT-PCR analysis was performed on a Bio-Rad iCycler System (CFX) using SsoAdvanced SYBR Green Supermix (Bio-Rad). qRT-PCR conditions were 30 s at 95 °C, followed by 45 cycles of 5 s at 95 °C and 20 s at 63 °C. RT controls were performed by the omission of RNA template or RT enzyme. A specific standard curve was performed in parallel for each gene, and each sample was quantified in duplicate. Data were analyzed by gene regression using iCycler software and normalized to Gapdh or 36B4 levels.
ChIP-seq. ChIP-seq data was previously generated 2 : ChIP-seq reads were aligned to the mouse reference genome (GRCm38/mm10) using Bowtie v0.12.8 using the following parameters -m 1 -strata -best. Only uniquely aligned reads have been retained for further analyses. Peak calling was performed using either MACS v1.4.2 or SICER v1.1 38 . MACS was used to detect peaks into RNAPII data using default parameters except for -g mm. SICER was used to detect islands into H3K27ac data using the script SICER.sh with the following parameters: Species: mm10, Threshold for redundancy allowed for chip reads: 1, Threshold for redundancy allowed for control reads: 1, Window size: 200 bps, Effective genome size as a fraction of the reference genome of mm10: 0.77, Gap size: 600 bps, Evalue for identification of candidate islands that exhibit clustering: 1000, False discovery rate controlling significance: 10 −2 . Fragment size was set according to the value assessed by Homer v4.7.2 makeTagDirectory. Peaks/Islands were annotated using Homer v4.7.2 with annotations extracted from Ensembl v78.
RNAseq. RNA-seq data was previously generated 2 : RNAseq reads were aligned onto mouse rRNA sequences using bowtie 39 release 0.12.7. Reads that do not map to rRNA sequences were mapped onto the mouse reference genome (GRCm38/mm10) using Tophat2 release 2.0.10 40 . Only uniquely aligned reads have been retained for further analyses.
Gene expression was quantified using HTSeq 41 release 0.5.4p3 and gene annotations from Ensembl release 78. Read counts were normalized across libraries with the method proposed by Anders et al. 42 . Comparison between R6/1 and WT samples was performed using the method proposed by Love et al. 43 , implemented in the DESeq2 Bioconductor library (release 1.0.19). Adjustment for multiple testing was performed using a method previously described 44 . eRNA identification. In order to identify reads from putative eRNA, we removed split-mapped reads and reads that overlap (≥ 1 bp on the opposite strand) genes annotated by Ensembl (release 78) as: IG_C_gene, IGD_gene, IG_J_gene, IG_LV_gene, IG_V_pseudogene, Mt_rRNA, Mt_tRNA, polymorphic_pseudogene, pro-tein_coding, pseudogene, rRNA, TR_V_gene and TR_V_pseudogene. We extended the region corresponding to those genes to 3 kb upstream of the transcription start site and 10 kb downstream of the transcript end site, in order to minimize signal from polymerase read-through from genic transcripts (as previously described ref. 19). IntersectBed from BEDTools release 2.21.0 was used for this purpose 45 . This overlap was performed on the SCIENtIfIC REpoRTS | 7:42875 | DOI: 10.1038/srep42875 opposite strand as the library preparation protocol we used to construct these RNAseq libraries leads to sequence the strand generated during first strand cDNA synthesis.
On those filtered RNAseq reads, we then detected the location of putative eRNA as genomic regions enriched in RNAseq reads. We thus used a method used for ChIP-seq peak detection: we used MACS v1.4.2 46 with the following parameters -keep-dup = all -nomodel -nolambda -p 1e-4 -g mm. The parameter -shiftsize was set according to the value assessed by Homer v4.7.2 47 makeTagDirectory.
Motif analysis. Motifs searching of known motifs (Jaspar 2014 motif database) was made using FIMO 48 within sequences of up and down eRNA. FIMO results were then processed by a custom Perl script that computed the number of occurrence of each motif. To assess the enrichment of motifs within the regions of interest, the same analysis was done n times (n = 100) on randomly selected eRNA regions. We chose to use randomly selected eRNA regions as controls to compute an expected distribution of motif occurrence and to correct for the nucleotide composition bias that could occur specifically in eRNA sequences. Region size distribution of the randomly selected eRNA was the same as for the up and down eRNA of interest. The significance of the motif occurrence was estimated through the computation of a Z-score.
The Z-score was computed this way: where: − x is the observed value (number of motif occurrence) − μ is the mean of the number of occurrences (computed on randomly selected data) − σ is the standard deviation of the number of occurrences of motifs (computed on randomly selected data) From the Z-scores, P values for each motif were computed. The P values were corrected for multi-testing using a method previously described 44 . Statistical analysis was done with custom R scripts.
Gene ontology (GO) analysis. Functional enrichments analyses were performed using the tools DAVID and/or GREAT 22,49 . For analyses with GREAT, defaults setting were used. Whole Mouse genome was used as background. Top-enriched terms are shown (P values < 0.05 were considered).
Data access. RNA-seq and ChIP-seq data used for eRNA analysis are accessible through GEO (accession number GSE59572).