Genome-wide histone acetylation analysis reveals altered transcriptional regulation in the Parkinson's disease brain.

BACKGROUND: Parkinson's disease (PD) is a complex, age-related neurodegenerative disorder of largely unknown etiology. PD is strongly associated with mitochondrial respiratory dysfunction, which can lead to epigenetic dysregulation and specifically altered histone acetylation. Nevertheless, and despite the emerging role of epigenetics in age-related brain disorders, the question of whether aberrant histone acetylation is involved in PD remains unresolved.
METHODS: We studied fresh-frozen brain tissue from two independent cohorts of individuals with idiopathic PD (n = 28) and neurologically healthy controls (n = 21). We performed comprehensive immunoblotting to identify histone sites with altered acetylation levels in PD, followed by chromatin immunoprecipitation sequencing (ChIP-seq). RNA sequencing data from the same individuals was used to assess the impact of altered histone acetylation on gene expression.
RESULTS: Immunoblotting analyses revealed increased acetylation at several histone sites in PD, with the most prominent change observed for H3K27, a marker of active promoters and enhancers. ChIP-seq analysis further indicated that H3K27 hyperacetylation in the PD brain is a genome-wide phenomenon with a strong predilection for genes implicated in the disease, including SNCA, PARK7, PRKN and MAPT. Integration of the ChIP-seq with transcriptomic data from the same individuals revealed that the correlation between promoter H3K27 acetylation and gene expression is attenuated in PD patients, suggesting that H3K27 acetylation may be decoupled from transcription in the PD brain. Strikingly, this decoupling was most pronounced among nuclear-encoded mitochondrial genes, corroborating the notion that impaired crosstalk between the nucleus and mitochondria is involved in the pathogenesis of PD. Our findings independently replicated in the two cohorts.
CONCLUSIONS: Our findings strongly suggest that aberrant histone acetylation and altered transcriptional regulation are involved in the pathophysiology of PD. We demonstrate that PD-associated genes are particularly prone to epigenetic dysregulation and identify novel epigenetic signatures associated with the disease.

Background

Parkinson’s disease (PD) is the second most common neurodegenerative disorder, affecting ~ 1.8% of the population above the age of 65 years [1]. The neuropathological hallmark of PD is the progressive loss of dopaminergic neurons in the substantia nigra pars compacta (SNc) in the presence of α-synuclein-positive inclusions termed Lewy pathology. Additional neurodegenerative changes occur in multiple regions of the nervous system, including several brainstem nuclei, the olfactory bulb, hippocampus, amygdala and the neocortex, as well as the autonomic and enteric nervous systems [2].

The etiology of PD remains largely unknown. While monogenic forms of PD exist, they are rare, generally accounting for less than 5–10% of all cases in most populations. The vast majority of patients have “idiopathic PD”, which is of complex etiology [3]. Genome-wide association studies (GWAS) have identified multiple risk loci associated with idiopathic PD. These make a relatively small collective contribution to the disease risk, however, and have not been linked to specific molecular mechanisms that could be targeted therapeutically [4].

Epigenetic modifications, such as DNA methylation and histone acetylation, modulate gene expression without altering the underlying DNA sequence. Dysregulation of histone acetylation has been previously linked to aging [5] and neurodegeneration [6] and changes in histone acetylation at specific genomic regions have been reported in Alzheimer’s disease [7, 8]. Despite the emerging role of epigenetics in age-related brain disorders, only a few small-scale studies of histone acetylation in PD have been reported to date, providing inconclusive results [9-11]. Moreover, the genomic landscape of histone acetylation in PD has never been assessed.

Here, we report the first genome-wide study of histone acetylation in the PD brain, based on a population-based cohort from the Norwegian Park West study (PW) [12]. Our findings replicate in an independent cohort of subjects from the Netherlands Brain Bank (NBB) and reveal genome-wide dysregulation of H3K27 acetylation in the PD brain.

Methods

Cohorts and tissue

Brains were collected at autopsy and split sagittally along the corpus callosum. One hemisphere was fixed whole in formaldehyde for ~ 2 weeks, and the other coronally sectioned and snap-frozen in liquid nitrogen. Samples from the frozen hemisphere were included in this study. Prefrontal cortex (Brodmann area 9) was obtained from two independent cohorts. The discovery cohort (ParkWest, or PW cohort) comprised individuals with idiopathic PD (n = 17) from the ParkWest study, a prospective population-based cohort which has been described in detail [12], and neurologically healthy controls (n = 15) from the Neuromics Biobank for Aging and Neurodegeneration. For a subset of the PW cohort (PD: n = 7, Control: n = 7), tissue samples from the striatum (putamen) and cerebellar cortex were also included. All individuals of the PW cohort had undergone whole-exome sequencing [13] and pathogenic mutations in genes implicated in Mendelian PD and other monogenic neurological disorders had been ruled-out. The replication cohort comprised individuals with idiopathic PD (n = 10) and neurologically healthy controls (n = 11) of Dutch origin, from the Netherlands Brain Bank (https://www.brainbank.nl/). Demographic information and experimental allocation of all samples is shown in Supplementary Table S1. For each experiment, there were no significant differences for sex, age, or postmortem interval (PMI) between the compared groups. A complete list of the medication used by the individuals with PD during at least the last 12 months before death can be found in Supplementary Table S1. Individuals with PD fulfilled the National Institute of Neurological Disorders and Stroke [14] and the U.K. PD Society Brain Bank [15] diagnostic criteria for PD at their final visit. All cases showed neuropathological changes consistent with PD (Supplementary Table S1).

Immunoblotting

Three independently dissected samples of brain tissue (replicates) were analyzed from each individual. Immunoblots were repeated at least 2 times for each marker. To normalize the histone acetylation signal to histone quantity, each acetylation blot was stripped twice. Complete stripping was confirmed by incubation of the membranes with secondary antibodies only and ECL. To exclude potential bias due to incomplete stripping, we performed additional experiments where acetylated and total histone markers for each histone modification were assessed in parallel on different blots. To estimate the acetylated fraction of each histone marker, we calculated the signal ratio of acetylated marker to total histone protein. This ratio is referred to as the “acetylation fraction” for each histone marker. Detailed description is provided in the Supplementary material, methods section. Comparison of immunoblot densitometric values between groups was performed using linear regression for pan-acetylation and linear mixed-models for all other markers. Linear mixed-models analysis was chosen in order to account for intraindividual random variability. Subjects were specified as random effects and disease state, sex, age of death, PMI and oligodendrocyte MSP were used as fixed effects. The error distribution was assumed normal, and solution to the mixed models was obtained using the maximum likelihood estimation method. The mixed-models analyses were implemented in SPSS version 22.0.0.

Chromatin Immunoprecipitation sequencing (ChIP-seq)

Approximately 200 mg of frozen PFC tissue per sample was used for ChIP-seq. Briefly, the tissue was crosslinked using formaldehyde, lysed using a Dounce homogenizer and sonicated to obtained chromatin fragment of 300–500 bp. Chromatin was immunoprecipitated using an antibody against acetylated H3K27. Protein-DNA complexes were recovered using protein A agarose beads, de-crosslinked and DNA was purified by phenol-chloroform extraction and ethanol precipitation. Following quantitative PCR, libraries were built and sequenced on Illumina’s NextSeq 500 (75 base pair reads, single end). Detailed description is provided in the Supplementary material, methods section.

ChIP-seq quality control and filtering

The analysis of the ChIP-seq data was performed in accordance with ENCODE ChIP-seq standards https://www.encodeproject.org/chip-seq/histone/#standards. Raw FASTQ files were trimmed using Trimmomatic version 0.38 [16] with options ILLUMINACLIP:TruSeq2-SE.fa:2:30:10 LEADING:3 TRAILING:3

SLIDINGWINDOW:4:15 MINLEN:25. Assessment of FASTQ files was carried out with FastQC version 0.11.8 [https://www.bioinformatics.babraham.ac.uk/projects/fastqc] before and after trimming. Reads were aligned to the GRCh37 human genome using bwa version 0.7.17 [17]. Reads with mapping quality < 30 were discarded using samtools version 1.9 [https://www.htslib.org/] and duplicated reads removed using MarkDuplicates from Picard Tools version 2.18.27 [https://broadinstitute.github.io/picard/]. BAM files were converted into BED files and sorted.

Peak calling and read counting

Peak calling procedure was carried out in two complementary approaches: 1) Per group (PD or controls): the peaks are called separately for each group, aggregating the relevant samples and inputs. While this approach is not optimal for differential peak calling analysis as it requires harmonization of peaks between the two groups, it can be used to compare the total genome coverage as well as unique H3K27ac genomic regions in each group. 2) All samples combined: the peaks are called by aggregating all samples and input controls. This approach does not allow identification of unique genomic regions however, it provides a better identification of consensus peaks, and thus, was used for differential peak analysis. The resulting peak-sets were filtered using the ENCODE black listed regions and only peaks within canonical chromosomes were kept for downstream analyses. Sample-specific enrichment on the identified peaks was performed using featureCounts version 1.6.4 [18] with default parameters. For both the individual group and all samples combined approaches, peaks with non-adjusted p-value > 10− 7 were excluded from the analysis. Reads inside peaks (RiP) were quantified using “featureCounts” program from the “subread” package v2.0.0. Additional details are provided in the Supplementary material, methods section.

Calculation of marker site profiles (MSP)

In order to identify brain cell-type specific H3K27ac regions we first analyzed H3K27ac ChIP-seq data from NeuN+ and NeuN− brain cells [19]. Cell-type based broadPeaks (CellType_peak-set), count matrixes and metadata files were downloaded from https://www.synapse.org/#!Synapse:syn4566010. Differential regions between NeuN+ and NeuN− cells were calculated using “DESeq2” R package [20], including chromatin amount, library batch, sex, hemisphere, age and pH as covariates in the model. Peaks were defined as cell-type specific differentially acetylated regions (DARs) if they met the following criteria: a) |fold change| > 4 and b) mean count > 1000. Peaks were annotated to all genes for which they intersected a region between 5 kb upstream from their transcription start site (TSS) to the end of their 5’UTR using “annotatr” R package [21], based on UCSC hg19 genome assembly. Genes with DARs were then itersected with expression-based cortical marker gene lists based on NeuroExpresso database [22]. The DARs were next reassigned to specific cell-types if they were annotated to genes defined as cell-type specific based on NeuroExpresso. For example, all DARs annotated to MBP, defined as oligodendrocyte marker gene based on NeuroExpresso, were defined as oligodendrocyte marker sites (MSS). In the next step, we quantified the reads from our samples in the CellType_peak-set. The corresponding RiP were then converted to counts per million (CPM) and transformed to log2(CPM + 1). Next, for each cell-type specific MSS, we performed principal component analysis based on the relevant peaks as described in [22]. Marker Site Profiles (MSP) were defined as the scores of the samples in the first principal component, transformed to [0,1] range for visualization purposes. Detailed description of the analysis is provided in the Supplementary material, methods section.

Identification of differentially acetylated regions

Differential peak analysis was performed using “nbinomWaldTest” function from “DESeq2” R package. Dispersions were calculated by using “estimateDispertions” function from “DESeq2” R package, by setting the fitType parameter to “local”. Normalization factors calculated by the “estimateSizeFactors” function were replaced with the housekeeping normalization factor based on our manually selected housekeeping gene set (ACTB, GAPDH, UBC). Peak significance was determined using DESeq2 “results” function, setting the cutoff for adjusted p-value to 0.05, and independentFiltering = T. Detailed description of the analysis is provided in the Supplementary material, methods section.

Association of H3K27 hyperacetylation with p300 and non-SIRT1 HDAC binding

ChIP-seq data (wgEncodeRegTfbsClusteredV3.bed) was obtained from the Encyclopedia of DNA Elements (ENCODE) project [23]. The data was available for p300 and the following HDACs: HDAC1, HDAC2, HDAC6, HDAC8 and SIRT6 We first quantified the number of binding sites for each protein in our peaks using the “findOverlaps” function from the “GenomicRanges” R package [24]. For each peak, the HDAC binding sites were then summed up to represent the non-SIRT1 HDAC_binding. Peaks that had no binding sites for either of the proteins were excluded from downstream analyses. For the gene level analysis, the number of binding sites in each peak was collapsed to a gene level.

Difference between DARs and non DARs was assessed by comparing the value of of the peaks using Kolmogorov-Smirnov test (all genes) or Wilcoxson’s rank sum test (PD_implicated genes) using “ks.test” and “wilcox.test” functions from R “stats” package. For linear model, we first obtained one-sided hyperacetylation p-value (phyper) for each peak as followed: . The association of hyperacetylation with p300_binding and HDAC binding was then assessed by the linear model: − log 10(p)~p300 + HDAC + PDgene where PDgene indicates whether a peak is annotated to PD_implicated gene or not.

Association of PD_implicated genes and EGL with altered H3K27 acetylation

Effective gene length (EGL) was defined as sum of length of all regions annotated to the gene (including 1Kb upstream to TSS). We assessed the association between altered H3K27 acetylation defined as -log10(metaP) with being annotated to a PD_implicated gene (PDgene), EGL and the total number of common regions aligned to the gene using a linear model. The total number of common regions in a gene was included as a covariate in the model. Namely, for a common region i, annotated to gene j, the -log (metaP) value was modeled as: − log 10(metaP)~PDgene +  log 10(EGL) +  ∑ CommonRegions. For the analysis was performed in each cohort separately, metaP was replaced with cohort-specific p-value, and the total number of common regions was replaced with the total number of regions (peaks) annotated to the gene in the cohort-specific peak-set.

Transcription factor enrichment analysis

The genomic sequence of hypo- and hyper-acetylated regions were repeat masked, and the enrichment of transcription factor binding motifs (TFBM) was assessed using the Analysis Motif of Enrichment (AME tool) [25] from the MEME Suite [26] v5.3.2, by setting the fraction of maximum log-odds to 0.75, and using human (HOCOMOCO v11) motif database [27]. The enrichment was assessed using Fisher’s exact test, with length-matched non-differentially acetylated regions set as the background.

RNA-seq analysis

Total RNA was extracted from prefrontal cortex tissue homogenate for all samples for which ChIP-seq was performed and sequenced following ribosomal RNA depletion A detailed description of the RNA-seq analysis is provided in the Supplementary material, methods section and is reported elsewhere [28]. Relative cell abundance was estimated using cell-type specific marker genes, as previously described [22, 29]. Since significant group effect was found for oligodendrocytes and microglia estimates in PW (Wilcoxon’s p = 2.6 × 10− 4) and NBB (Wilcoxon’s p = 0.009) cohorts respectively, these cell types were included at covariates in the analysis of both cohorts. Count matrix and the code necessary to reproduce these results are available on https://git.app.uib.no/neuromics/cell-composition-rna-pd.

Promoter H3K27ac - expression correlation

ChIP-seq and RNA-seq counts were first adjusted to the relevant covariates (age, sex, PMI, batch, oligodendrocyte estimates and microglia estimates). The adjustment was done by re-fitting the counts using the coefficients calculated by DESeq2 R function with adjusted model matrix in which categorical covariates were set to 0, and the continuous covariates to the mean of the original values. Next, for each gene, Pearson’s correlation between the adjusted ChIP-seq/RNA-seq seq counts was computed using “cor.test” R function, in each group (controls or PD subjects) separately. Only peaks annotated to promoter regions of the genes were included in the analysis. When multiple promoter peaks were defined for a gene, the correlation was calculated for each promoter separately, and the highest correlation was kept. This step was performed independently for each group. Since not all genes are expected to exhibit a good correspondence between promoter H3K27ac state and expression level due to existence of multiple levels of regulation, we compared the distribution of correlations in each group for different thresholds of absolute correlation in either of the groups (e.g., cor > |0.1| means that only genes with absolute correlation > 0.1 in at least one of the groups are included in the analysis). Group comparison was performed using Wilcoxon rank sum test, using “wilcox.test” function from R “stats” function. Functional enrichment analysis was performed using the gene score resampling method implemented in the ermineR package version 1.0.1 [30], an R wrapper package for ermineJ [31], setting the value of “BigIsBetter” to FALSE, with the complete Gene Ontology (GO) database annotation [32]. For the enrichment of genes based on the decoupling levels, gene scores were defined as ∆Cor = Cor − Cor. For enrichment of genes based on the lowest correlation in PD group, gene scored were defined as the correlation values in the PD group.

Data and code availability

The R code necessary to reproduce the results presented in the paper is available on https://github.com/ltoker/ChIPseqPD. Raw counts for the ChIP-seq analyses are available in the same repository (https://github.com/ltoker/ChIPseqPD/tree/master/data). UCSC genome browser tracks showing the H3K27ac peak-sets based on PW and NBB cohorts can be accessed at https://genome-euro.ucsc.edu/s/ltoker/H3K27ChIP_PD. Additional R code to reproduce the RNA-seq analyses is available on https://git.app.uib.no/neuromics/cell-composition-rna-pd. The corresponding raw data counts are included in the repository (https://git.app.uib.no/neuromics/cell-composition-rna-pd/tree/master/Data).

Results

Increased histone acetylation and histone protein levels in the PD brain

We first assessed global protein acetylation in the prefrontal cortex (PFC) of 13 individuals with PD and 13 demographically matched controls from the PW cohort by immunoblot analysis using an antibody against acetylated lysine (Fig. 1a, Supplementary Table S1). We concentrated on PFC because, while this area is typically affected in PD [33-35], it is generally less confounded by altered cell composition due to neurodegeneration and neuroinflammation [2]. Linear regression analysis of the densitometric measurements, adjusted for age, sex and PMI, indicated an increased acetylation signal at ~ 17 kDa, consistent with the molecular weight of core histone proteins (p = 0.014, Fig. 1a).

Fig. 1

Hyperacetylation of multiple histone sites in the PD prefrontal cortex. Representative Western blots and the ratios of their quantified densitometric measurements of the indicated proteins (or their modifications) in individuals with PD (orange) and controls (blue) are shown. GAPDH served as loading control. Values were rescaled [0,1]. Boxplots show the mean of replicates for each individual. The indicated p-values are based on a linear mixed model (see Methods). a Global protein acetylation detected by a pan-acetyl lysine antibody. b Relative acetylation of histone lysine residues and quantification of the total histone proteins. c Relative acetylation of α-tubulin on lysine K40. d Immunodetection and densitometric quantification of relative protein levels of the nuclear marker (NUP188) and cell type-specific marker proteins: NeuN – neurons, GFAP – astrocytes, CNP1 – oligodendrocytes, CX3CR1 - microglia

To identify histone sites contributing to the hyperacetylation observed in our samples and select a specific site for downstream analyses, we assessed the acetylation state at multiple sites on the four core histones in 7 PD and 7 control samples showing a prominent difference in the pan-acetylation analysis (Fig. 1b, Supplementary Table S1). The acetylated fraction of each histone site and the total histone quantity were analyzed using linear mixed models with gender, age and PMI as fixed effects and subject and blot number as random effects. These analyses indicated a significant increase in H2BK15ac (p = 1.4 × 10− 5), H3K9/14 ac (p = 3.9 × 10− 4), H3K27ac (p = 4.5 × 10− 8), H3K56ac (p = 0.01) and H4K12ac (p = 5.0 × 10− 3) in PD. No significant difference was found for H2AK5, H4K5, H4K16 or α-tubulin, used as a non-histone control (Fig. 1b-c, Supplementary Table S2). In addition, we found increase in the total quantity of the four histone proteins (H2A: p = 2.7 × 10− 5, H2B: p = 1.1 × 10− 8, H3: p = 2.0 × 10− 7; H4: p = 3.9 × 10− 18).

To control that our findings were not driven by differences in cell type composition, we carried out immunoblot analyses of protein markers for neurons (NeuN), astrocytes (GFAP), oligodendrocytes (CNP1), microglia (CX3CR1) and total number of nuclei (NUP188), in the same samples used for acetylation assessment. Statistical analyses of these data indicated increase in CNP1 (p = 3.9 × 10− 5, Fig. 1d, Supplementary Table S2), but no changes in the other markers. We thus repeated the histone acetylation analyses adding CNP1 as a covariate to the statistical model. As shown in Supplementary Table S2, all of the indicated hyperacetylated sites remained significant after CNP1 adjustment.

Similar to findings from PFC, H3K27 hyperacetylation was also observed in the PD striatum (p = 0.034) and cerebellum (p = 0.016), areas typically spared by neuronal loss [2], but exhibiting disease-related pathology, including Lewy pathology [36, 37] and mitochondrial respiratory deficiency [38] (Fig. 2, Supplementary Table S2).

Fig. 2

Altered H3 acetylation occurs in multiple brain regions, is accompanied by sirtuin upregulation and is not induced by anti-Parkinson drugs in-vitro. a-c Representative Western blots and quantification plots showing the ratios of quantified densitometric measurements of the indicated proteins in individuals with PD (orange) and controls (blue). GAPDH and β-tubulin serve as loading controls. Values were rescaled [0,1]. Boxplots show the mean of replicates for each individual. The indicated p-values are based on a linear mixed model (see Methods). a, b Acetylated H3K9/14 and H3K27 and total histone protein H3 levels in the cerebellum (a) and striatum (b) from individuals with PD and controls. c Sirtuin proteins in the prefrontal cortex of PD and controls. d Acetylated H3K27 in SH-SY5Y cells, untreated (UT) or treated with different concentrations of the sirtuin inhibitor Sirtinol, or the solvent (DMSO) for 6 h. A representative blot is shown. Boxplots show the quantified ratios of three biological replicates (indicated in different colors) carried out in duplicates, rescaled [0,1]. e Acetylated H3K9/14, H3K27 and total H3 protein in differentiated SH-SY5Y-derived dopaminergic neuron-like cells treated with the histone deacetylase inhibitor Trichostatin A (TSA) or anti-Parkinson drugs used by the PW patients during their last year of life. The numbers indicate the concentration (in μM) of the drug. Carb: Carbidopa; Enta: Entacapone; Ropi: Ropinirole; DOPA: L-Dopa. UT: untreated control. A representative blot is shown. Bar plots show the mean quantified ratios of three biological replicates (rescaled [0,1]). Dots indicate the values of each of the individual replicates

Histone hyperacetylation in the PD brain is accompanied by upregulation of sirtuin protein levels

Histone hyperacetylation can be induced by mitochondrial respiratory complex-I deficiency [10, 39, 40], a pathological feature found throughout the PD brain, including the PFC, striatum and cerebellum [38]. This effect is thought to be mediated by decreased NAD+/NADH ratio, attenuating the activity of NAD+-dependent class III histone deacetylases known as sirtuins [39, 41–44]. Direct confirmation of this mechanism in our samples would require a reliable determination of their intracellular NAD+, which was not feasible due to rapid postmortem degradation [45]. However, it was previously reported that decreased sirtuin activity caused by NAD+ deficiency is accompanied by increased protein quantity [46, 47]. Thus, we assessed the protein levels of SIRT1, SIRT2 and SIRT3, localized to the nucleus, cytosol and mitochondria, respectively, as surrogates for their activity. We were specifically interested in these proteins, because they have been linked to aging and neurodegeneration [41]. Immunoblot analyses of the 14 samples assessed for histone acetylation, indicated a significant upregulation of SIRT1 (p = 2.0 × 10− 3) and SIRT3 (p = 3.0 × 10− 3) in PD. No difference was detected for SIRT2 (Fig. 2c, Supplementary Table S2), in line with the unaltered acetylation state of its target, α-tubulin (Fig. 1c). Moreover, treatment of SH-SY5Y cells with the sirtuin inhibitor sirtinol led to a substantial hyperacetylation of H3K27, in comparison to control cells (Fig. 2d), indicating that sirtuins can indeed modulate the acetylation of this residue.

Histone hyperacetylation is not induced by anti-Parkinson drugs in cellulo

We next investigated whether the observed altered histone acetylation in PD could result from treatment with anti-Parkinson agents. To test this, we exposed fully differentiated SH-SY5Y cells to the drugs used by most individuals with PD from the PW cohort during their last year of life (L-DOPA, carbidopa, entacapone and ropinirole), and monitored acetylation of H3K9/14 and H3K27 by immunoblot (Supplementary methods). These analyses indicated no effect for any of the tested anti-Parkinson drugs on histone acetylation in our model (Fig. 2e). Nevertheless, these findings cannot completely rule out the possibility that altered histone acetylation can be, at least partially, induced by drug-related effects in the human brain over a prolonged period of time.

Acetylated H3K27 exhibits higher genomic occupancy in PD

We next carried out ChIP-seq analysis of acetylated H3K27 (H3K27ac) in all PFC samples from the PW cohort (n = 28), Supplementary Table S1). First, we characterized the genomic distribution of H3K27ac in the PD and control groups. For this purpose, we identified H3K27ac enriched genomic regions (peaks) in each group separately (see Methods). In line with the increased H3K27ac signal observed in the immunoblot analysis, a higher number of peaks (146,763 vs. 135,662), as well as a higher proportion of group-unique peaks (32.2% vs. 13.2%), and wider genome coverage (10.6% vs. 9.11%) were found in the PD- compared to the control group (Fig. 3a, Supplementary Table S3). We then compared the distribution of genomic annotations and detection p-values of overlapping and unique peaks between the PD and control groups. In both groups, unique peaks occurred more frequently in intergenic and intronic regions than in promoter and exonic regions (Supplementary Fig. S1a, b). Since unique peaks were also characterized by significantly higher detection p-values, compared to the overlapping peaks (Fig. S1c, d), we believe that most of these represent “noisy” H3K27 acetylation.

Fig. 3

H3K27ac ChIP-seq analysis of PD and control samples. a Different peak-sets used for the different stages of ChIP-seq analysis. Shown, as an example, is a genomic region around the oligodendrocyte marker MBP. Peaks annotated to this region, in each peak-set are shown by colored boxes. The CellType_peak-set was used for the MSP calculation. Group specific peak-sets were used to characterize the distribution of H3K27ac bound regions in cases and controls. Barplots show the comparison between the two groups (Supplementary Table S3). PD group exhibited higher number of peaks, higher percentage of unique peaks and higher percentage of genomic coverage, consistent with genome-wide H3K27 hyperacetylation. The PW_peak-set was used for the differential peak analysis. Shown are group overlays of sample fold of change compared to input control. b Hierarchical clustering of the samples indicated that samples cluster based on their cellular composition. Supplementary Fig. S6 shows the association of each of the variables with the main principal components of the data. c MA plot based on the PW_peak-set indicates that increased H3K27ac in PD is observed genome-wide, rather than being restricted to specific regions. d Manhattan plot showing the distribution of genomic locations and differential p-values of the H3K27ac PW_peak-set. The red dashed line indicates the -log10 of the highest p-value with < 5% false positive rate (Benjamini-Hochberg adjusted p-value < 0.05). Corresponding plots for the NBB cohort are shown in Supplementary Fig. S5

ChIP-seq data corroborates H3K27 hyperacetylation in the PD brain

Next, we redefined H3K27ac peaks by identifying H3K27ac enriched regions by aggregating all PW samples (see Methods). The final peak-set consisted of 160,521 peaks, covering 11.6% of the genome and distributed over all 22 autosomes and both sex chromosomes (median peak size, 817 bp, Fig. 3d). All downstream analyses were performed based on this peak-set, henceforth referred to as PW_peak-set.

ChIP-seq data reflects the abundance and genomic distribution of the target protein (or protein modification) which can vary substantially between samples [48]. Based on the observed increase in H3K27ac fraction in the immunoblot analysis, we anticipated an increase in overall H3K27ac ChIP-seq signal in PD. Thus, we reasoned that the most suitable approach for normalization of our data would be using reads in peaks (RiP) annotated to house-keeping genes, as previously proposed [48]. We modified this approach by manually selecting peaks annotated to known housekeeping genes that exhibit a saturated peak signal covering most of the gene body (Supplementary Fig. S2). The rationale here was that these genes should be saturated with histone acetylation and are, therefore, expected to exhibit little biological variation.

To validate our normalization approach, we performed additional H3K27ac immunoblot and ChIP-seq analyses in samples from 55 individuals (including the study subjects, Supplementary Table S1) and tested the correlation between the immunoblot and ChIP-seq data normalized by different approaches. As shown in Supplementary Fig. S3, the highest correlation between immunoblotting and ChIP-seq data was observed when the total number of RiP was normalized to our manually chosen peaks. Importantly, comparison of normalized RiP between PD and control samples corroborated the finding of general H3K27 hyperacetylation in PD (Supplementary Fig. S3).

Cell composition is a major source of variation in bulk tissue H3K27ac ChIP-seq data

Since different cell-types exhibit distinct epigenetic landscapes [19], bulk tissue data must be evaluated in the context of the underlying cellular composition [22, 49, 50]. To estimate the cellular composition in our ChIP-seq data, we used an approach similar to the one described by Mancarci et al. [22]. Briefly, we intersected genomic H3K27ac binding sites differentiating between NeuN-positive and NeuN-negative brain cells [19] with the list of brain cell-type specific marker-genes [22], in order to obtain genomic H3K27ac cell-type marker sites. The first principal component of reads aligned to the marker sites was used to obtain Marker Site Profiles (MSPs) and was used as a proxy for the relative abundance of the corresponding cell-types across the samples. We validated this approach in publicly available ChIP-seq data from cortical and cerebellar bulk tissue samples [51], and ChIP-seq data from entorhinal cortex of individuals with Alzheimer’s disease and healthy controls [8]. MSP analysis of these data successfully recapitulated the well-known increased abundance of glial cells in cortical vs cerebellar samples [52, 53], and decreased number of neurons in the entorhinal cortex of subjects with Alzheimer’s disease compared to controls [54, 55] (Supplementary Fig. S4a, b).

Next, we performed MSP analysis in the PW_peak-set data. This analysis indicated no significant differences in cell-type abundance between the PD and control groups (Supplementary Fig. S5). Nevertheless, since the cellular composition of cortical samples can be affected by variable white matter contamination introduced during dissection [28], we assessed the effect of oligodendrocyte MSPs and other experimental and demographic covariates on the ChIP-seq data. Sample-to-sample correlation indicated that the samples clustered mainly based on oligodendrocyte MSPs (Fig. 3b). Accordingly, the first principal component of the data was mostly associated with oligodendrocyte MSP (Supplementary Fig. S6), confirming that oligodendrocyte composition is a major source of variance in our data. Therefore, oligodendrocyte MSPs were accounted for in all downstream analyses.

Hyperacetylation of H3K27 in the PD brain is a genome-wide phenomenon

We next sought to identify differentially acetylated regions (DARs) for H3K27 between PD and control samples. After filtering (see Methods), 133,716 peaks annotated to 17,182 genes remained for downstream analyses. Five samples (one control and four PD samples) were detected as outliers and excluded from differential peak analysis (Supplementary Table S1). Differential peak analysis was performed using the “DESeq2” R package, adjusting for age, sex, batch, PMI, oligodendrocyte MSP and normalization peak ratio. DARs were defined as differentially acetylated peaks with adjusted p-value < 0.05.

Our analysis identified 2877 hyperacetylated and only 14 hypoacetylated PD-associated DARs, annotated to 1434 and 9 genes, respectively. The top-ranked hyperacetylated and hypoacetylated DARs resided within the DLG2 (adjusted p = 9.8 × 10− 4) and PTPRH (adjusted p-value 9.0 × 10− 3) genes, respectively, both of which have been previously associated with PD [4, 56–58]. The full list of DARs with their annotated genes is provided in Supplementary Table S4. An MA-plot (the log ratio vs. mean count) of all regions included in the analysis indicated a genome-wide trend for H3K27 hyperacetylation in PD (Fig. 3c).

H3K27 hyperacetylation in the PD cortex is replicated in an independent cohort

To validate our findings, we carried out ChIP-seq in an independent cohort from the Netherlands Brain Bank (NBB, Supplementary Table S1). In agreement with our finding from the PW cohort, we observed an increase in the percentage of unique peaks as well as increased total genomic coverage in the PD group (Supplementary Table S3). The distribution of genomic annotations among the unique and overlapping peaks and the distribution of peak p-values were similar between the two cohorts (Supplementary Fig. S1). As in the PW cohort, cellular composition accounted for most of the interindividual variance in this cohort (Supplementary Fig. S6b, S7a).

Differential acetylation analysis corroborated the finding of PD-associated genome-wide H3K27 hyperacetylation observed in the PW cohort (Supplementary Fig. S5b). In total, the analysis identified 2486 hyperacetylated DARs, mapped to 946 genes, and 227 hypoacetylated DARs, mapped to 253 genes (Supplementary Table S5). Strikingly, 275 out of the 946 genes with hyperacetylated DARs and 2 out of 253 genes with hypoacetylated DARs replicated across both cohorts (Fig. 4c, hypergeometric test: p = 2.49 × 10− 85 and p = 8.2 × 10− 3, respectively). These 277 genes are henceforth referred to as “replicated genes”. The top replicated genes (ranked based on their adjusted p-value in the PW cohort) are shown in Table 1 and a complete list is provided in Supplementary Table S6.

Fig. 4

Replicated genes and replicated DARs. a, b Schematic representation of two levels of replication between the PW and NBB cohorts. H3K27ac peaks are illustrated by colored lines, representing control (blue) and PD (orange) samples, with the height of the peak representing the peak intensity. Shown is a hypothetical gene X, with exons represented as grey boxes and introns represented by lines. In both cases, the gene has two pairs of common peaks: Peak_a/A, mapped to the first exon of the gene and Peak_b/B, mapped to the third exon of the gene. a Altered H3K27 acetylation occurs at different regions of gene X in each cohort. Thus, while gene X harbors altered H3K27 acetylation in both cohorts (i.e. gene X is a replicated gene), it does not harbor replicated DARs. b Altered H3K27 acetylation occurs at the same region of gene X in both cohorts. In this case, gene X is a replicated gene with a replicated DAR. c Venn diagram showing the number of replicated genes (with or without replicated DARs) between the two cohorts, out of all genes represented in our ChIP-seq data from both cohorts. d Venn diagrams showing the number of replicated DARs between the two cohorts, out of all common peak pairs (left) and their corresponding annotated genes (right). e Venn diagram showing the overlap of PD-implicated genes with DARs between cohorts. All DARs exhibited hyperacetylation. f Venn diagram showing the number of PD-implicated genes with adjusted metaP < 0.05

Table 1

Top replicated genes between PW and NBB cohorts

GeneSymbol	log2(FC)_PW	padj_PW	log2(FC)_NBB	padj_NBB
Hypoacetylated genes
PTPRH	−0.81	9.1E-03	−0.94	0.03
JUP	−0.43	0.05	−0.44	0.02
Hyperacetylated genes
DLG2	0.68	9.8E-04	0.77	2.0E-03
TNRC6B	0.74	9.8E-04	0.9	0.03
SIK3	0.66	3.2E-03	0.52	0.01
SGK1	0.65	3.2E-03	0.72	0.01
PPP2R2B	1.3	3.5E-03	0.7	8.0E-03
FCHSD2	0.83	3.5E-03	1	0.01
CTNND2	0.52	3.5E-03	0.63	0.01
PTK2	0.71	3.5E-03	0.55	0.02
SLCO1A2	0.58	3.5E-03	0.68	0.02
LPAR1	0.99	3.5E-03	0.55	0.04
FRMD4B	0.64	3.5E-03	0.53	0.04
RNF216	0.6	5.1E-03	0.59	2.3E-03
TMEM132B	0.88	5.1E-03	0.59	0.02
PIP4K2A	0.57	5.1E-03	0.52	0.03
MAP7	0.55	5.9E-03	0.52	0.04
PDE4B	0.59	6.5E-03	0.67	2.4E-03
FRMD5	0.55	6.5E-03	0.67	0.01

Meta-analysis of ChIP-seq data from PW and NBB cohorts

In order to perform meta-analysis of the ChIP-seq data we first harmonized the data from the two cohorts by pairing peaks with the same structural and functional genomic annotation (e.g. promoter, exon, UTR). These peaks are henceforth referred to as “common peaks”. For each common peak pair, we then calculated the Fisher’s meta p-value (metaP). A toy example of a common peak pair is shown in Fig. 4a, b.

In total, the harmonization process resulted in 52,286 pairs of common peaks. Out of these, 7288 (mapped to 2972 genes) had adjusted metaP < 0.05 (Supplementary Table S7, Fig. 4d, f). Moreover, 300 pairs represented peaks defined as DARs in both cohorts (Fig. 4d, hypergeometric p-value = 5.29 × 10− 37). These 300 common peak pairs, mapped to 197 genes, are henceforth referred to as “replicated DARs”. Nearly all (295/300) replicated DARs were hyperacetylated in both cohorts and only a single DAR, mapped to the PTPRH gene, was hypoacetylated. Genes annotated to the top replicated DARs, ranked based on their adjusted metaP, are provided in Table 2.

Table 2

Genes annotated to the top replicated DARs

Symbol	Region Annotation	CHR	RegionLocation	PeakLocation (PW, NBB)	MetaP	LFCPW, NBB	PeakNamePW, NBB	ProportionSignificant
Common hypoacetylated regions
PTPRH	introns	19	55,718,568–55,720,750	55,720,288–55,722,414, 55,717,912–55,719,103	1.7E-07	−0.81, −0.94	Peak_20188, Peak_15767	3p/4 (0.8)
JUP	introns	17	39,681,526–39,923,630	39,805,764–39,808,133, 39,807,092–39,808,006	1.1E-05	−0.43, − 0.99	Peak_53983, Peak_81536	5p/10p (0.5)
Common hyperacetylated regions (top 15)
DLG2a	exons	11	83,874,504–83,874,554	83,873,551–83,879,844, 83,873,666–83,877,544	7.7E-12	0.68, 1.05	Peak_21808, Peak_40741	30p/49p(0.6)
KIT	exons	4	55,604,595–55,606,881	55,606,747–55,609,274, 55,606,782–55,609,553	8.1E-11	0.63, 1.26	Peak_21777, Peak_37103	2/3p (0.7)
CNTNAP2	exons	7	147,336,198–147,336,398	147,333,973–147,338,260, 147,333,944–147,337,213	6.2E-09	0.51, 1.07	Peak_44999, Peak_47588	11/18p (0.6)
RNF216	introns	7	5,681,008–5,692,043	5,679,685–5,685,127, 5,688,370–5,705,729	1.2E-08	0.6, 0.53	Peak_38631, Peak_12618	11p/16 (0.7)
FRMD5	introns	15	44,216,508–44,393,656	44,260,324–44,263,821, 44,244,249–44,247,825	1.6E-08	0.68, 0.67	Peak_18018, Peak_18744	12p/17p(0.7)
SOX2-OT	introns	3	181,328,718–181,417,385	181,358,880–181,360,150, 181,341,154–181,348,014	1.7E-08	0.55, 0.91	Peak_43667, Peak_31549	9/12 (0.8)
MAGI2	1to5kb	7	78,120,338–78,124,337	78,121,257–78,128,843, 78,115,190–78,121,170	2.2E-08	0.4, 0.93	Peak_18816, Peak_24042	9p/23 (0.4)
PPP2R2B	introns	5	146,080,706–146,460,621	146,383,635–146,384,203, 146,459,075–146,469,815	2.8E-08	1.3, 0.52	Peak_135972, Peak_17645	17p/27p(0.6)
KCNH1	introns	1	210,948,887–210,970,849	210,966,438–210,969,215, 210,946,491–210,950,181	3.8E-08	0.55, 1.0	Peak_51838, Peak_38212	9/18p (0.5)
CDH1	introns	16	68,835,797–68,842,326	68,840,159–68,841,130, 68,782,065–68,837,136	4.8E-08	0.84, 0.37	Peak_66179, Peak_3389	5/7p (0.7)
MAP 4 K4a	introns	2	102,315,001–102,407,181	102,330,963–102,334,475, 102,367,013–102,375,949	5.2E-08	0.54, 0.62	Peak_30093, Peak_9009	10p/24p(0.4)
METTL24	introns	6	110,644,076–110,729,536	110,656,829–110,661,962, 110,685,537–110,688,887	5.4E-08	0.48, 0.97	Peak_22739, Peak_27591	8/13p(0.6)
PIP4K2A	introns	10	22,898,647–23,003,111	22,934,570–22,952,399, 22,957,818–22,962,517	5.6E-08	0.58, 0.52	Peak_22480, Peak_27954	9p/18p(0.5)
SIK3	introns	11	116,827,781–116,968,858	116,825,819–116,835,556, 116,846,923–116,884,577	5.8E-08	0.6, 0.49	Peak_21504, Peak_11226	13p/16p(0.8)
DOCK10	introns	2	225,796,386–225,906,968	225,798,681–225,800,264, 225,845,772–225,849,871	8.1E-08	0.61, 0.77	Peak_94440, Peak_22866	6/10p(0.6)
Common hyperacetylated regions (additional GWAS hits)
CRHR1a	introns	17	43,707,525–43,750,192	43,726,544–43,728,643, 43,741,886–43,747,152	2.0E-07	0.43, 0.83	Peak_74390, Peak_30944	6/14p(0.4)
MBNL2a	introns	13	97,873,814–97,927,885	97,919,075–97,922,481, 97,906,667–97,916,758	1.2E-05	0.44, 0.44	Peak_29656, Peak_16279	14p/22p(0.6)

For hyperacetylated regions, the top 15 genes with DARs common to the two cohorts and common DARs in genes previously associated with PD are shown. The genes were ranked based on Fisher’s meta p-value. The full list of genes is provided in Supplementary Table S6

aPrevious association with PD based on GWAS; Region Location – genomic location (on the indicated chromosome) of the annotated region; Peak location – genomic location of the relevant peak from each cohort (upper – PW, lower - NBB); LFC – log fold change compared to controls; Peak Name – the peak identifier in each cohort; Proportion significance – DARs/total common regions annotated to the gene; p- peaks annotated to promoter region of the gene

Altered H3K27 acetylation has a predilection for genes with an established link to PD

We next wondered whether genes with an established involvement in PD are over-represented among DAR-harboring genes. To assess this, we compiled a list of genes that fulfilled at least one of the following criteria: (a) are associated with PD based on the most recent and largest GWAS meta-analysis [4], (b) are implicated in monogenic PD [59, 60], or (c) encode proteins with a central role in PD-related neuropathology [2, 61, 62]. Among the 92 genes meeting these criteria (Supplementary Table S8), 83 could be assessed in ChIP-seq data from both cohorts and are henceforth referred to as “PD-implicated genes”. Out of these, six genes: DLG2, MAP 4 K4, CRHR1, MBNL2, SH3GL2 and APP, harbored DARs in both cohorts (hypergeometric p-value = 2.9 × 10− 3, Fig. 4e). In addition, 24/83 PD-implicated genes had adjusted metaP < 0.05 (hypergeometric p-value = 6.9 × 10− 3, Supplementary Table S8, Fig. 4f). These included genes with key-roles in the pathogenesis of idiopathic and/or monogenic PD: SNCA, MAPT, APP, PRKN and PARK7. Interestingly, the significantly hyperacetylated peak in SNCA was in an enhancer region of the gene previously shown to be affected by both genetic variation [63] and drug exposure associated with PD [64] (Fig. 5a).

Fig. 5

H3K27 acetylation is associated with the relative proportion of p300 and non-SIRT1 HDAC binding and overrepresented in PD-implicated genes. a, b Examples of two PD_implicated genes, SNCA (a) and VPS35 (b). SNCA harbors a DAR (Peak_35531/27168 in PW/NBB respectively) in an enhancer element whose H3K27 acetylation status is associated with both genetic variation (SNP rs356168) and drugs influencing the risk of PD. VPS35 does not exhibit DARs in our data. Purple bars indicate H3K27ac peaks in the PW and NBB peaksets. H3K27ac ChIP-seq data from ENCODE is shown for reference. Pink, blue and yellow bars indicate the binding sites of p300, HDAC1 and HDAC2 respectively, based on ENCODE ChiP-seq data wgEncodeRegTfbsClusteredV3.bed. c, d Difference in the relative number of p300 and non-SIRT1 HDAC bindings sites between significantly hyperacetylated and non-significantly hyperacetylated peaks for all (c, Kolmogorov–Smirnov test) or only PD_implicated (d, two-sided Wilcoxon rank sum test) genes. e Hyperacetylated – regions showing hyperacetylation in both cohorts; Mixed – regions with opposite acetylation trends in the two cohorts. The dashed line represents the -log10 of the highest metaP corresponding to adjusted metaP < 0.05. Multiple common regions within the same gene are shown. Highlighted are regions identified as DAR in at least one of the cohorts within PD-implicated genes. metaP – Fisher’s meta-p value of the corresponding two peaks from PW and NBB cohorts. f Association between -log(p) and each of the indicated covariates based on linear model. Shown are the coefficients and 95% confidence intervals. The results are shown for the common regions (metaP) and each of the cohorts separately. Peaks annotated to miRNAs and snRNAs were excluded from the analysis. PDgene – a gene is a PD-implicated gene. EGL – effective gene length. TotalPeaks – total number of common regions or peaks annotated to the gene

H3K27 hyperacetylated regions are associated with an increased ratio of p300 binding sites relative to non-SIRT1 HDAC binding sites

Our immunoblot analyses suggested that H3K27 hyperacetylation may be accompanied by decreased sirtuin activity, primarily SIRT1. It is well established that SIRT1 can also deacetylate and thus suppresses the activity of the histone acetyl transferase p300 [65], the primary mediator of H3K27 acetylation [66]. Thus, we reasoned that if H3K27 hyperacetylation in PD is indeed mediated by decreased SIRT1 activity, it should be associated with regions exhibiting increased binding of p300 and decreased binding of non-SIRT1 HDACs. To test this, we integrated our data with the relevant ChIP-seq data derived from the ENCODE project [23] (see Methods). Specifically, we quantified the number of p300 bound regions (p300_binding) or non-SIRT1 HDAC bound regions (HDAC_binding) from ENCODE that overlap with H3K27ac peaks in our data and assessed the association of the two measures with H3K27 hyperacetylation level of each region. In both cohorts, the ratio of p300 binding sites over non-SIRT1 HDAC binding sites was significantly higher in regions overlapping with hyperacetylated DARs compared to non-DARs (Fig. 5c) . Notably, this was also true when restricting the analyses to the PD-implicated genes (Fig. 5d). Linear model analysis showed positive association of H3K27 hyperacetylation with p300_binding (PW: β = 0.02, p < 2 × 10− 15; NBB: β = 0.09, p < 10− 26) and negative association with HDAC_binding (PW: β = − 0.07, p < 10− 26; NBB: β = − 0.14, p < 10− 26). Similar results were obtained when peaks were collapsed to a gene level (Supplementary Table S9).

Differential H3K27 acetylation is associated with effective gene length

We noted that many of the genes with DARs, including PD-implicated genes, are among the longest in the genome. We reasoned that this may be because longer genes have higher histone binding capacity and therefore are more likely to show a differential signal in the event of a generalized, unspecific H3K27 hyperacetylation. To test this hypothesis and ensure that the over-representation of PD-implicated genes among DAR-harboring genes is not driven by longer gene length, we used linear models to assess whether H3K27ac, defined as -log10(metaP), is associated with PD-implicated genes while accounting for the effective gene length (EGL), defined as the total length of all regions annotated to a gene (see Methods). The total number of common regions in a gene was included as a covariate in the model. For hyperacetylated regions, this analysis indicated a significant positive association between altered acetylation and PD-implicated genes (p = 1.0 × 10− 7) as well as EGL (p < 2.0 × 10− 26). For the hypoacetylated regions, however, the analysis indicated a significant negative association with both PD-implicated genes (p < 2.0 × 10− 26) and EGL (p < 2.0 × 10− 26). Analysis of each cohort independently corroborated these findings, although the negative association with PD-implicated genes for the hypoacetylated regions was only significant in the NBB cohort (Fig. 5b, Supplementary Table S10).

DARs differ in their GC content and are enriched in multiple transcription factor motifs

In both cohorts, the hyperacetylated DARs exhibited low GC content compared to non-differentially acetylated regions, matched for region length (NBB: 44% vs 51%; PW: 46% vs. 50%, DARs vs. non-DARs, respectively, Supplementary Fig. S8). The opposite phenomenon was observed for the hypoacetylated DARs (NBB: 61% vs 51%; PW: 55% vs. 50%, DARs vs. non-DARs, respectively, Supplementary Fig. S8). Similar results were obtained when comparison was performed only for peaks annotated to either exonic or intronic regions. Motif enrichment analysis of the hyper- and hypoacetylated DARs indicated enrichment in multiple TFBMs. For the hyperacetylated DARs, significant enrichment was found for 71 and 94 TFBMs, in the PW and NBB cohorts respectively, with 64 TFBMs overlapping between the two cohorts. For the hypoacetylated DARs, significant enrichment was found for 3 and 49 TFBMs, in the PW and NBB cohorts respectively, two of which were enriched in both (Supplementary Table S11).

Promoter H3K27 acetylation is decoupled from gene expression in PD

While it is well accepted that modifications of H3K27ac are predictive of gene expression levels [67], it was previously shown that widespread histone hyperacetylation induced by potent HDAC inhibition does not lead to major changes in gene expression [68, 69]. To investigate whether a similar decoupling occurs in the PD brain, we analyzed RNA sequencing (RNA-seq) data from the PFC of the same individuals included in the ChIP-seq analysis (Supplementary Table 1, Supplementary Methods). For each gene, we then calculated the Pearson’s correlation between the level of its expression and promoter H3K27ac across controls or individuals with PD.

In the control group from both cohorts, most genes exhibited the expected positive correlation between promoter H3K27ac and expression. In the PD group, however, the distribution of correlations was centered closer to 0, regardless of the absolute correlation threshold chosen (Fig. 6a-e). The deviation between the two groups progressively increased as only genes with higher correlation between promoter H3K27ac and expression were included in the analysis (Fig. 6b-e). Importantly, similar results were obtained based on RLE-normalized ChIP-seq data (Supplementary Fig. S9), indicating that the observed decreased correlation between H3K27ac and transcription levels in PD was not driven by the normalization method. These results suggest that the coupling between promoter H3K27 acetylation state and proximal gene transcription is impaired in the PFC of individuals with PD.

Fig. 6

Decreased correlation between promoter H3K27 acetylation state and gene expression in PD Pearson’s correlation between the adjusted promoter H3K27ac counts and the adjusted RNA-seq counts was calculated for each gene across control or PD subjects in each cohort separately. a Representative scatter plots showing the adjusted promoter H3K27ac counts (X-axis), the adjusted RNA-seq counts (Y-axis) and the calculated Pearson’s correlation for three genes: MED13, PTPRH and DLG2 in control (blue) or PD (orange) individuals from the PW or NBB cohorts. Each dot represents an individual. Cor - Pearson’s correlation. b-e The distribution of correlations for all compatible genes was compared between the groups for various thresholds of minimal absolute correlation in either of the groups. Results are shown for PW (b, c) and NBB (d, e) cohorts. The number of genes that remained after applying the correlation threshold is shown in parentheses. The correlations remained closely distributed around 0 in PD subjects from both cohorts, regardless of the correlation threshold. c, e Median correlations in each group and the calculated Wilcoxon’s delta shift for the different correlation thresholds. f Nucleus-encoded mitochondrial genes exhibit negative correlation between promoter H3K27 acetylation and gene expression. Shown are examples for three mitochondria related GO terms

Functional enrichment analysis based on the decoupling level of the genes, indicated that nucleus-encoded mitochondrial genes are the most decoupled genes in both cohorts (Supplementary Table S12). Moreover, closer examination of these genes revealed that most of them exhibit negative correlation between the expression level and promoter H3K27 acetylation (Fig. 6f, Supplementary Table S12).

Discussion

Aberrant H3K27ac occurs genome-wide in the PD brain

We report evidence of increased histone acetylation in the PFC of individuals with PD with the strongest change observed for H3K27, a marker of active promoters and enhancers with a fundamental role in regulating gene expression [70, 71]. H3K27 hyperacetylation in PD is supported by multiple levels of evidence in our study: (a) immunoblot analyses showed a highly significant increase in the acetylated fraction of H3K27 across three brain regions in PD. Importantly, we show that this is not driven by underlying differences in cell composition of the samples and is unlikely to be induced by anti-Parkinson medication. (b) ChIP-seq analysis, replicated in two independent cohorts, showed that the total number of peaks, percentage of unique peaks, genome coverage and the total number of RiPs are all increased in the PD samples. The analysis further revealed that increased H3K27ac occupancy has a genome-wide distribution. (c) We found a strong association between differential H3K27 acetylation and the effective length of the targeted gene. This association was positive for hyperacetylated regions and negative for hypoacetylated regions, in line with a hyperacetylation-predisposing cellular environment in PD.

H3K27 hyperacetylation may be induced by altered sirtuin activity

We hypothesized that H3K27 hyperacetylation in PD may be mediated by altered sirtuin activity, due to a decreased NAD+/NADH ratio [39, 42, 72] resulting from complex-I deficiency, a pathological hallmark of the PD brain [38, 73]. Hyperacetylation of several histone sites, including H3K27, was indeed observed in animal and cell studies, following chemical complex-I inhibition with PD-associated neurotoxins [10, 74]. In line with our hypothesis, we observed increased levels of SIRT1 and SIRT3 proteins in PD samples and no change in SIRT2 quantity, or in the acetylation state of its substrate, α-tubulin. This is consistent with the notion that SIRT1 and SIRT3, but not SIRT2 activity are susceptible to changes in physiological NAD+ concentrations [72]. In line with our findings, reduced sirtuin activity accompanied by increase in their protein levels was reported in iPSC-derived dopaminergic cells carrying a LRRK2 mutation [75]. Moreover, complex-I deficiency, decreased NAD+/NADH ratio, decreased sirtuin activity (accompanied by increased protein levels), and increased histone acetylation have all been shown to occur with aging [5, 41, 44, 46, 47], which is the strongest known risk factor for PD.

The putative role of SIRT1 in mediating H3K27 hyperacetylation in PD is further supported by the observation that hyperacetylation is more pronounced among regions with high p300/non-SIRT1 HDAC binding site ratio. This is because decreased SIRT1 activity increases p300-mediated H3K27 acetylation, and can be compensated by other HDACs. Namely, H3K27 hyperacetylation induced by decreased SIRT1 activity is expected be more pronounced in regions exhibiting multiple p300 binding sites and fewer non-SIRT1 HDAC binding sites. While our findings strongly suggest the sirtuins may be involved in H3K27 hyperacetylation in PD, further studies are required to validate this effect, and to establish whether this involvement is direct (i.e. by deacetylation of H3K27) or indirect (e.g., via reduced deacetylation and thereby activation of p300).

H3K27 hyperacetylation is over-represented in PD-implicated genes

While H3K27 hyperacetylation appears to be a genome-wide phenomenon in the PD brain, it also shows a highly significant predilection for genes implicated in PD. Meta-analysis of our two cohorts identified H3K27 hyperacetylated regions in 24/83 genes linked to idiopathic and/or monogenic forms of PD (Supplementary Table S8). These include 21 genes associated with PD by GWAS [4], genes encoding proteins involved in PD pathology (SNCA, MAPT, APP) [2, 61, 62] and genes in which mutations cause Mendelian PD (SNCA, PRKN, PARK7) [76, 77]. In addition, H3K27 hyperacetylation was observed in genes associated with more complex forms of parkinsonism and/or vulnerability to PD-related pathology such as FBOX7 [78] and POLG [59].

Several mechanisms may underlie the observed over-representation of PD-implicated genes. First, it is plausible that H3K27 acetylation at promoter and/or enhancer regions of these genes is modulated by genetic variation and/or environmental factors associated with PD. For instance, the H3K27 acetylation state of the SNCA enhancer region, significantly hyperacetylated in our data, was shown to be affected by both genetic variation [63] and drug exposure associated with PD [64]. Second, it is likely that pathogenic processes involved in PD induce the expression of genes implicated in the disease. For example, the expression levels of PRKN and PARK7 increase in response to oxidative stress [79-82]. Since, the regulation of transcription often involves histone acetylation, impairment of the deacetylation machinery, would over time lead to excessive H3K27 acetylation of these genes. Finally, by relaxing chromatin structure, H3K27 hyperacetylation in genomic regions associated with PD, predisposes these areas to somatic mutations, including genomic rearrangements [83]. As previously suggested [84], this could act synergistically with other predisposing factors, such as age and oxidative injury, eventually leading to PD-specific pathology.

H3K27 hypoacetylation in the PTPRH gene may contribute to PD pathophysiology

Considering the general H3K27 hyperacetylation in the PD brain, the robust hypoacetylation of H3K27 in the PTPRH gene, encoding protein tyrosine phosphatase receptor type H, is of particular interest. While little is known about the physiological function of PTPRH, it was previously found to be enriched in loss-of-function mutations in PD [57]. Furthermore, it was shown that RNAi-mediated downregulation of the PTPRH homolog in Drosophila enhances α-synuclein neurotoxicity [57]. Our finding of selective H3K27 hypoacetylation of the PTPRH region in PD, corroborates the hypothesis that decreased function of this gene may contribute to the pathophysiology of PD. Additional studies are required to establish the physiological function of PTPRH and its involvement in the pathophysiology of PD.

Differential H3K27ac is associated with GC content

The observation that hyperacetylated DARs exhibit low GC content, whereas hypoacetylated DARS exhibit high GC content, is intriguing. While the reason behind this phenomenon is unclear, the existence of enhancers with high and low GC content was previously described, linking low GC content-enhancers to the immune response [85]. Additional studies are required to understand the mechanism underlying the selectivity of differential H3K27ac towards different sequence composition and its link to PD.

Decoupling between promoter H3K27 acetylation and gene expression in the PD brain

The attenuated correlation between promoter H3K27 acetylation and gene expression is intriguing, since it indicates that a fundamental mechanism of gene regulation is defective in the PD brain. A plausible explanation for this phenomenon is a global hyperacetylation of enhancer regions, compromising the regulation of gene expression. Indeed, in both cohorts, about half of the identified H3K27ac peaks as well as DARs were annotated to non-coding regions. Furthermore, since histone acetylation is tightly linked to the metabolic state of the cell [72], this decoupling suggests that cellular responses to metabolic stress may be severely compromised in PD. This putative compromised response to metabolic stress is further supported by the finding that, in both cohorts, the highest degree of decoupling is observed among nuclear-encoded mitochondrial genes. Moreover, the negative correlation between H3K27 acetylation and gene expression of these genes supports the involvement of sirtuins, in particular SIRT1, in the pathophysiology of PD. This is because in addition to histones, SIRT1 deacetylates and activates PGC-1α, the key regulator of mitochondrial biogenesis [42, 72, 86]. Thus, decreased SIRT1 activity would simultaneously lead to increased histone acetylation and reduced expression of nuclear-encoded mitochondrial genes, resulting in negative correlation between their expression level and promoter H3K27 acetylation state. Moreover, decreased expression of mitochondrial genes would initiate a vicious circle by decreasing NAD+/NADH ratio and further reducing SIRT1 activity.

Conclusions

We show that H3K27 acetylation, a fundamental mechanism modulating gene expression, is severely dysregulated in the PFC of individuals with idiopathic PD. The highly significant concordance between the two cohorts and the clear predilection for genes strongly linked to PD, validate our findings and suggest that epigenetic dysregulation plays an important role in the pathogenesis of the disease. Importantly, since histone acetylation and NAD metabolism can be pharmaceutically modulated, our results open the possibility that agents targeting these biological processes may have therapeutic potential for PD and should be further explored.

Additional file 1.