A histone acetylome-wide association study of Alzheimer’s disease: neuropathology-associated regulatory variation in the human entorhinal cortex

Alzheimer’s disease (AD) is a chronic neurodegenerative disorder characterized by the progressive accumulation of amyloid-β (Aβ) plaques and neurofibrillary tangles in the neocortex. Recent studies have implicated a role for regulatory genomic variation in AD progression, finding widespread evidence for altered DNA methylation associated with neuropathology. To date, however, no study has systematically examined other types of regulatory genomic modifications in AD. In this study, we quantified genome-wide patterns of lysine H3K27 acetylation (H3K27ac) - a robust mark of active enhancers and promoters that is strongly correlated with gene expression and transcription factor binding - in entorhinal cortex samples from AD cases and matched controls (n = 47) using chromatin immunoprecipitation followed by highly parallel sequencing (ChIP-seq). Across ~182,000 robustly detected H3K27ac peak regions, we found widespread acetylomic variation associated with AD neuropathology, identifying 4,162 differential peaks (FDR < 0.05) between AD cases and controls. These differentially acetylated peaks are enriched in disease-specific biological pathways and include regions annotated to multiple genes directly involved in the progression of Aβ and tau pathology (e.g. APP, PSEN1, PSEN2, MAPT), as well as genomic regions containing variants associated with sporadic late-onset AD. This is the first study of variable H3K27ac yet undertaken in AD and the largest study investigating this modification in the entorhinal cortex. In addition to identifying molecular pathways associated with AD neuropathology, we present a framework for genome-wide studies of histone modifications in complex disease, integrating our data with results obtained from genome-wide association studies as well as other epigenetic marks profiled on the same samples.

Alzheimer's disease (AD) is a chronic neurodegenerative disorder characterized by cognitive decline and memory loss that contributes substantially to the global burden of disease, affecting in excess of 26 million people worldwide 1 . The symptoms of AD are associated with progressive neuropathology in the neocortex, with regions surrounding the entorhinal cortex being particularly affected early in the disease 2 . These neuropathological hallmarks of AD include the extracellular deposition of neurotoxic amyloid-β (Aβ) in the form of amyloid plaques and an accumulation of intracellular neurofibrillary tangles composed of hyperphosphorylated tau 3 . Despite progress in understanding risk factors contributing to AD progression, the mechanisms involved in disease progression are not fully understood and long-term treatments, reversing the cellular disease process in the cortex, are elusive. There has been considerable success in identifying genetic risk factors for AD 4 . While autosomal dominant mutations in three genes (APP, PSEN1, and PSEN2) can explain earlyonset (< 65 years) familial AD, these account for only 1-5% of the total disease burden 5 .
Most cases of AD are late-onset (> 65 years), non-Mendelian and highly sporadic, with susceptibility attributed to the action of highly prevalent genetic variants of low penetrance.
In addition to the well-established risk associated with the APOE locus 6, 7 there has been notable success in identifying novel AD-associated variants capitalising on the power of genome-wide association studies (GWAS) in large sample cohorts; a recent large GWAS meta-analysis of AD, incorporating > 74,000 samples, identified 19 genome-wide significant risk loci for sporadic AD 8 . Despite these advances, little is known about the functional mechanisms by which risk variants mediate disease susceptibility.
Increased understanding about the functional complexity of the genome has led to growing recognition about the likely role of non-sequence-based regulatory variation in health and disease. Building on the hypothesis that epigenomic dysregulation is important in the aetiology and progression of AD neuropathology 9 , we and others recently performed the first genome-scale cross-tissue analyses of DNA methylation in AD identifying robust DNA methylation differences associated with AD neuropathology across multiple independent human post-mortem brain cohorts 10,11 . To date, however, no study has systematically examined other types of regulatory genomic modifications in AD. In this study, we focus on lysine H3K27 acetylation (H3K27ac), a robust mark of active enhancers and promoters that is strongly correlated with gene expression and transcription factor binding 12 . Interestingly, histone deacetylase (HDAC) inhibitors have been shown to ameliorate symptoms of cognitive decline and synaptic dysfunction in mouse models of AD 13,14 and are promising targets for novel human AD treatments 15,16 . Despite this, investigations into global levels of histone acetylation in AD have thus far been inconclusive [17][18][19][20] and no study has taken a genome-wide approach. In fact, few studies have systematically profiled H3K27ac across large numbers of samples in the context of complex disease, and optimal methods for these analyses are still being developed 21 .
In this study, we used chromatin immunoprecipitation combined with highly-parallel sequencing (ChIP-seq) to quantify levels of H3K27ac across the genome in post-mortem entorhinal cortex samples from AD patients and matched controls. We identify regulatory genomic signatures associated with AD, including variable H3K27ac across discrete regions annotated to genomic loci mechanistically implicated in the onset of both tau and amyloid pathology. This is the first study of variable H3K27ac yet undertaken for AD; in addition to identifying molecular pathways associated with AD neuropathology, we present a framework for genome-wide studies of this modification in complex disease.

Genome-wide profiling of inter-individual variation in H3K27ac in the entorhinal cortex
We generated H3K27ac ChIP-seq data using post-mortem entorhinal cortex tissue dissected from 47 elderly individuals (average age = 77.43, SD = 9.66, range = 58-97) comprising both AD cases (n = 24, mean Braak stage = 6.00, SD = 0.00) and age-matched low pathology controls (n = 23, mean Braak stage = 1.30, SD = 1.11) (Supplementary Table 1). Raw H3K27ac ChIP-seq data is available for download from the Gene Expression Omnibus (GEO) (accession number GSE102538). After stringent quality control (QC) of the raw data (see Methods), we obtained an average of 30,032,623 (SD = 10,638,091) sequencing reads per sample, with no difference in read-depth between AD cases and controls (P = 0.93, Supplementary Fig. 1). This represents, to our knowledge, the most extensive analysis of H3K27ac in the human entorhinal cortex yet undertaken. Using combined data from all 47 samples (see Methods) we identified 182,065 high confidence H3K27ac peaks; these are distributed across all 24 chromosomes (Supplementary Table 2) spanning a mean length of 983bp (SD = 682bp) with a mean distance between neighbouring peaks of 15,536bp (SD = 116,040bp). We validated the identified peaks using two independent ChIP-seq datasets: first, we obtained locations for frontal cortex (BA9) and cerebellum H3K27ac peaks from a recent analysis of autism and control brains 21 ; second, we downloaded peak profiles for multiple cell-and tissue-types from the NIH Epigenomics Roadmap Consortium 22 (see Methods). As expected, there was a very high overlap between H3K27ac peaks called in other neocortical datasets and our ChIP-seq data, with a notably lower overlap with H3K27ac data from non-cortical tissues (Supplementary Fig. 2).

AD-associated differential acetylation in the entorhinal cortex
We quantified read counts across every peak in each individual sample using HTSeq 23 and employed a quasi-likelihood F test, implemented in the Bioconductor package EdgeR 24 , to test for differences in H3K27ac between AD cases and low pathology controls (see Methods). Our analysis model controlled for age at death and neuronal cell proportion estimates derived from DNA methylation data generated on the same samples using the CETS R package 25 (Supplementary Table 1, Supplementary Fig. 3). Subsequent sensitivity analyses excluded an effect of sex on our results, but revealed widespread effects of cell-type proportion estimates and age at death highlighting the importance of including these covariates (Supplementary Fig. 4). A total of 4,162 (2.3%) of the 182,065 peaks were characterized by AD-associated differential acetylation at a false discovery rate (FDR) < 0.05 ( Fig. 1), with a significant enrichment of hypo-acetylated AD-associated peaks (2,687 (1.5%)) compared to hyper-acetylated AD-associated peaks (1,475 (0.8%)) (P = 9.89E-80, exact binomial test) (Fig. 1). UCSC Genome Browser tracks showing H3K27ac levels in AD cases and controls, in addition to association statistics, across the genome can be accessed at https://tinyurl.com/y8f3tmoa. The ten top-ranked hyper-and hypo-acetylated peaks associated with AD are shown in Table 1, with a complete list given in Supplementary Peaks were subsequently annotated to genes using the Genomic Regions Enrichment of Annotations Tool (GREAT) 26 , which takes into account the strength of proximal and distal DNA-binding events. The most significant AD-associated hyperacetylated peak (chr13: 112101248-112102698; P = 2.04E-08; log fold change = 0.93) is annotated to both SOX1 and TEX29 on chromosome 13 (Fig. 2, Table 1). Of note, H3K27ac data from the Epigenomics Roadmap Consortium show that this region is characterized by brain-specific enhancer activity (Supplementary Fig. 5). The most significant AD-associated hypoacetylated peak (chr7: 64011549-64012825; P = 1.66E-08; log fold change = -0.86) is located within intron 1 of ZNF680 on chromosome 7 (Fig. 3, Supplementary Fig. 5, Table   1). Global clustering of samples by normalized read counts across all hyper-and hypoacetylated peaks (FDR < 0.05) indicated that, as expected, samples group primarily by disease status (Fig. 4). AD-associated differentially acetylated peaks (FDR < 0.05) are significantly longer (1295bp vs 975bp, P < 1.00E-50) and characterized by higher readdepths (2.56 log counts per million (CPM) vs 1.59 log CPM, P < 1.00E-50) than nonsignificant peaks (Supplementary Fig. 6). Of note, within AD-associated peaks, hypoacetylated peaks are significantly longer (1412bp vs 1081bp, P = 5.66E-31) and have higher read-depths (2.21 log CPM vs 1.77 log CPM, P = 2.69E-50) compared to hyperacetylated peaks. We used RSAT 27,28 to identify enriched transcription factor binding motifs located within AD-associated differentially acetylated peaks (see Methods), observing a significant enrichment of binding motifs for specificity protein 1 (Sp1) (P < 1.0E-50) amongst AD-hyperacetylated peaks (FDR < 0.05). Of note, previous publications have reported dysregulated expression of Sp1 and its co-localization with neurofibrillary tangles in AD 29,30 .
Increased H3K27ac is observed in regulatory regions annotated to genes previously implicated in both tau and amyloid neuropathology One of the top-ranked AD-associated hyper-acetylated peaks is located proximal to the gene encoding microtubule associated protein tau (MAPT) (chr17: 43925717-43927482; P = 7.01E-07; log fold change = 0.71; Table 1), which is widely expressed in the nervous system where it functions to promote microtubule assembly and stability. Tau is believed to play a key role in AD neuropathology, with hyperphosphorylation of the tau protein precipitating the neurofibrillary tangles associated with the pathogenesis of AD 31,32 . Closer inspection of the region around this AD-associated peak highlighted an extended cluster of six hyper-  Fig. 7). Strikingly, AD-associated differentially-acetylated peaks were also found in the vicinity of other genes known to play a direct mechanistic role in AD. We identified a significantly hypoacetylated peak (chr21: 27160993 -27161475; P = 3.94E-04; log fold change = -0.72) on chromosome 21, located ~100kb downstream of the amyloid precursor protein gene (APP), which encodes the precursor molecule to Ab, the main component of amyloid plaques [34][35][36] (Supplementary Fig. 8). We also identified significant hyperacetylation in the vicinity of the presenilin genes PSEN1 and PSEN2, which encode integral components of the gamma secretase complex and play a key role in generation of Aβ from APP 37 . In PSEN1 we found significantly elevated H3K27ac across a peak within intron 6 (chr14: 73656445 -73656860; P = 3.44E-04; log fold change = 0.68; Supplementary Fig. 9). In PSEN2 we identified consistent hyperacetylation in AD cases across nine H3K27ac peaks (FDR < 0.05) spanning a ~57 kb region upstream of the  Table 6). Of note, highly-penetrant mutations in APP, PSEN1, and PSEN2 are associated with familial forms of early-onset AD 38 . The identification of altered regulation of these loci in late-onset sporadic AD brain further supports a key role for altered amyloid processing in the onset of neuropathology.

Specific differentially-acetylated peaks also overlap known AD GWAS regions
Using data from a large GWAS meta-analysis of AD 8 we tested for instances where there is an overlap between AD-associated differential H3K27ac and genomic regions harbouring risk variants. Briefly, we defined linkage-disequilibrium (LD) blocks around the genome-wide significant (P < 5.0E-08) GWAS variants identified by the stage 1 meta-analysis by Lambert and colleagues 8 (Supplementary Table 7), which contained a total of 292 overlapping entorhinal cortex H2K27ac peaks (see Methods). Two of the 11 GWAS LD blocks contained significant AD-associated H3K27ac peaks (FDR < 0.05), although there was no overall enrichment of AD-associated differential acetylation at the 292 peaks (P = 0.364, Wilcoxon rank-sum test). Two peaks of AD-associated hyperacetylation were located within a GWAS

Integrative analysis of DNA and histone modifications reveal unique distributions of DNA modifications across regions of differential acetylation
Our previous work identified cortex-specific variation in DNA methylation (5mC) robustly associated with AD pathology 10,11 . We were therefore interested in exploring the relationship between H3K27ac and both 5mC and another DNA modification -DNA hydroxymethylation (5hmC), which is enriched in the brain and believed to play an important role in neuronal function, learning and memory 50-52 -in AD. Both modifications were profiled using DNA isolated from the same entorhinal cortex samples using oxidative bisulfite (oxBS) conversion in conjunction with the Illumina 450K HumanMethylation array ("450K array") (see Methods). 6,838 loci mapped to within 1kb of differentially acetylated peaks (FDR < 0.5; mapping to 1,649 unique peaks). We tested differential 5mC and 5hmC associated with AD at these probes, controlling for age at death and cell-type proportion estimates. None of the differences in 5mC (minimum P = 2.47E-03) or 5hmC (minimum P = 1.53E-03) were significant when correcting for multiple testing (n = 6,838 tests; P < 7.31E-05), indicating that there is little direct overlap in AD-associated variation in H3K27ac and DNA modifications.
Furthermore, comparing effect sizes at these 6,838 peak-probe pairs identified no evidence for a correlation between AD-associated H3K27ac and 5mC differences (r = 0.009, P = 0.443; Supplementary Fig. 13) with a small, but significant, negative correlation for 5hmC (r H3K27ac being localized at active enhancers and promoters.

Discussion
We quantified H3K27ac across the genome in post-mortem entorhinal cortex tissue samples, identifying widespread AD-associated acetylomic variation. Strikingly, differentially acetylated peaks were identified in the vicinity of genes implicated in both tau and amyloid neuropathology as well as genomic regions containing variants associated with sporadic late-onset AD. This is the first study of variable H3K27ac yet undertaken for AD; in addition to identifying molecular pathways associated with AD neuropathology, we introduce a framework for genome-wide studies of this modification in complex disease.
Given its close relationship with transcriptional activation, for example via the mediation of transcription factor binding, the identification of AD-associated variation in H3K27ac highlights potential novel regulatory genomic pathways involved in disease etiology. We find widespread alterations in H3K27ac associated with AD, including in the vicinity of several genes known to be directly involved in the progression of Ab and tau pathology 31, 53 (APP, PSEN1, PSEN2, MAPT), supporting the notion that dysregulation of both pathways is involved in the onset of AD. Interestingly, although our study assessed brains from donors affected by sporadic late-onset AD, we identify widespread altered H3K27ac in the vicinity of genes implicated in familial early-onset AD. This indicates that these two forms of the disease may share common pathogenic pathways and mechanisms. Given that histoneacetylation modifiers are amongst the most promising target pharmacological treatments of AD 15,54 , the identification of altered H3K27ac in AD is important, giving clues as to which genes and pathways may be involved.
Our study has a number of limitations, which should be considered when interpreting these results. First, we undertook ChIP-seq using bulk entorhinal cortex samples comprising a mix of neuronal and non-neuronal cell-types. This is an important limitation in epigenomic studies of a disease characterized by cortical neuronal loss. However, we were able to control for variation in neuronal proportions in our samples by deriving neuronal proportion estimates for each sample using DNA methylation data generated on the same tissue samples 25 . Second, our cross-sectional analysis of post-mortem brain tissue makes direct causal inference difficult, and it is likely that many of the changes in H3K27ac we observe result from the AD pathology itself. In this regard, however, it is interesting that we see disease-associated H3K27ac in the vicinity of genes causally implicated in familial forms of AD. Third, we have assessed a relatively small number of samples. In this light, it is notable that we identify substantial differences between AD cases and controls, with diseaseassociated regulatory variation in genes and functional pathways known to play a role in the onset and progression of neuropathology. The clear clustering between patients and controls at our differentially acetylated peaks suggests that, despite a complex and heterogeneous etiology, AD may be characterized by a common molecular pathology in the entorhinal cortex, reflecting neuropathological analyses. Fourth, chromatin architecture and transcriptional regulation is influenced by a multitude of epigenetic mechanisms. Although profiling H3K27ac can provide relatively robust information about transcriptional activity, it represents only one of perhaps ~100 post-translational modifications occurring at > 60 histone amino-acid residues regulating genomic function.
In summary, we provide compelling evidence for widespread acetylomic dysregulation in the entorhinal cortex in AD. Our data suggest that regulatory variation at multiple loci, including in the vicinity of several known AD risk genes -APP, CR1, MAPT, PSEN1, PSEN2 and TOMM40 -is robustly associated with disease, supporting the notion of common molecular pathways in both familial and sporadic AD. In addition to identifying molecular pathways associated with AD neuropathology, we present a framework for genome-wide studies of histone modifications in complex disease, integrating our data with results obtained from genome-wide association studies as well as other epigenetic marks profiled on the same samples.

Samples
Post-mortem brain samples from 54 individuals -27 with advanced AD neuropathology and  Table 1). Subjects were approached in life for written consent for brain banking, and all tissue donations were collected and stored following legal and ethical guidelines (NHS reference number 08/MRE09/38; the HTA license number for the LBBND brain bank is 12293). All samples were dissected by trained specialists, snap-frozen and stored at −80 °C. A detailed list of demographic and sample data for each individual included in the final analyses is provided in Supplementary Table 1.

Chromatin immunoprecipitation
Chromatin immunoprecipitation was performed using the iDeal ChIP-Seq kit for Histones (Cat# C01010051, Diagenode, Seraing, Belgium) as detailed below, using the standard kit components unless otherwise stated. 30 mg of entorhinal cortex tissue was homogenized with a dounce homogenizer in 1 mL ice-cold phosphate buffered saline (PBS) buffer with protease inhibitor cocktail (PIC). The suspension was centrifuged at 4,000 rpm for 5 minutes at 4°C, discarding the supernatant. The pellets were resuspended in 1 mL PBS containing 1% formaldehyde, rotating at room temperature for 8 minutes. The cross-linking process was terminated by adding 100 μL glycine solution, followed by 5 minutes of rotation. After 5 minutes of centrifugation at 4,000 rpm and 4°C, the pellet was washed twice with ice-cold PBS (suspending the pellet in 1 mL PBS with PIC, centrifuging for 5 minutes at 4,000 rpm and 4°C, and discarding the supernatant), then lysed in 10 mL ice-cold lysis buffer iL1 and iL2, sequentially (re-suspending the pellet in 10 mL lysis buffer, mixing gently for 10 minutes at 4°C, centrifuging for 5 minutes at 4,000 rpm and 4°C, and discarding the supernatant).

Data pre-processing and quality control
Global sample anomalies were ruled out using fastqc 56  and individual read counts did not associate with disease status (P = 0.93).

Peak calling and read counts
All filtered BAM files were merged into one grouped file and converted to tagAlign format using bedtools 59 . Peaks were called on this merged file using MACS2 60 , keeping all duplicates, since duplicates were removed from each sample previously and any remaining duplicates would result from the same read occurring in more than one sample. From the resulting peaks those located in unmapped contigs and mitochondrial DNA were filtered out as well as peaks that did not meet a significance threshold of P < 1.0E-07 for peak calling.
The bed file of peaks was converted to gff format using awk and R, and reads for each individual sample were generated using HTSeq 23 . Final filtering was performed using the Bioconductor package EdgeR 24 , excluding peaks with fewer than 2 samples showing at least 1 read per million, resulting in a total of 182,065 peaks to be tested. Principal components analysis (PCA) in R using DESeq2 61 confirmed that the epigenetically predicted gender was identical to the recorded one (Supplementary Fig. 15), with load on the first two principal components not related to disease status.

Peak validation
We validated the 182,065 union peaks in two ways. First, we obtained the locations of H3K27ac peaks called in the cortex (BA9) and cerebellum from a recent paper by Sun and colleagues 21  https://www.ncbi.nlm.nih.gov/geo) for multiple cell-/tissue-types including several brain regions (mid frontal lobe (GSM773015), inferior temporal gyrus (GSM772995), middle hippocampus (GSM773020), substantia nigra (GSM997258), cingulate gyrus (GSM773011), H1-derived neuronal progenitor cells (HDNPs, GSM753429), lung (GSM906395), liver (GSM1112808) and skeletal muscle (GSM916064)). The downloaded files were in bed format, on which we performed peak calling using MACS2 and the same specifications as described for our own samples, discounting any duplicate reads. We calculated the overlap between each peak set and our peaks by quantifying the percentage of peaks from the external sample overlapping our peaks using the Bioconductor package GenomicRanges 62 .

Differential peak calling
We used the quasi-likelihood F test 63 in EdgeR 24 to analyse peak differences between ADcases and controls, allowing us to correct for potential confounders in the analysis of differential peaks. Our analyses accounted for additional phenotypic variation across the samples, including age at death and neuronal proportion estimates based on DNA methylation profiles from the Illumina 450K HumanMethylation Array from the same samples, which were calculated using the CETS R package 25 . To test for the effects of these two covariates and sex we ran three further analyses, omitting age at death and CETS and adding sex as a covariate to our main model. We imputed the median CETS estimate for one individual with missing DNA methylation data. To include age at death and CETS estimates in the EdgeR differential peak calling method, these variables were converted to five-level factors using the R function cut(). The distribution of the age and CETS variable (including the imputed individual) with the respective bins of the factor variables are shown in Supplementary Fig. 3. We next calculated normalization factors based on samplespecific library compositions and estimated both sample and peak-specific dispersions, specifically for a generalized linear model controlling for factorized CETS estimates and age at death. The quasi-likelihood F-test was conducted after fitting a quasi-likelihood model 63 using the glmQLFit() and glmQLFTest() functions respectively. Effect sizes are reported as log fold change, a standard measure for quantifying sequencing read count differences between different conditions. Log fold change refers to the log 2 -transformed ratio of normalized read counts between cases and controls, with positive values indicating higher normalized read counts in AD samples. The bedtools program genomecov was used to generate coverage value scaled by library size and the number of samples per group, for each sample. These were then joined using unionbedg and summed using a Perl script to produce a weighted mean for each variable sized interval defined by read overlaps.

Genomic annotation and enrichment analyses
Peaks were annotated to genes using the Genomic Region Enrichment and Annotation Tools (GREAT) 26

Motif enrichment analysis
Motif analysis was performed using the Regulatory Sequence Analysis Tools suite (RSAT) 27,28 , available at http://rsat.sb-roscoff.fr. Peak sequences were reduced to 500bp on each side of the peak centre, and motif discovery was conducted on 6 and 7mer oligonucleotides, comparing the statistically enriched sequences with known transcription factor motifs from JASPAR 64 (core nonredundant vertebrates) and Homer 65 (Human TF motifs). Enrichments were computed relative to the background peak sequences (n = 182,065 peaks) for significantly hyper-and hypoacetylated peaks (FDR < 0.05).

Analysis of differential H3K27ac across AD regions from genome-wide association studies (GWAS)
The summary statistics for the stage 1 GWAS from Lambert and colleagues 8 were downloaded from http://web.pasteur-lille.fr/en/recherche/u744/igap/igap_download.php.
These results were clumped (p1 = 1e-4; p2 = 1e-4, r2 =0.1, window = 3000kb) using plink 66 , which collapses multiple correlated signals (due to linkage disequilibrium (LD)) into regions which represent independent signals. LD relationships were inferred from a reference GWAS dataset (Phase 1) from another study 67 . Neighbouring regions located within 250kb of each other on the same chromosome were subsequently merged. After clumping, each region was assigned the minimum P value for all SNPs contained in the region (from Lambert et al), and regions were then filtered to the genome-wide significance threshold (P < 5.0E-08). This yielded 11 LD blocks for the genome-wide significant findings from Lambert et al., which were then overlapped with our AD-associated differentially acetylated peaks using the Bioconductor package GenomicRanges 62 .

Integrative analysis with DNA methylation and hydroxymethylation
DNA methylation and hydroxymethylation data was available (unpublished) from entorhinal cortex DNA for 42 of the samples profiled in this study. DNA methylation and hydroxymethylation profiles were generated on the Illumina Infinium HumanMethylation450 BeadChip (Illumina Inc., CA, USA) ("Illumina 450K array") using the TrueMethyl Array kit (Cambridge Epigenetix, Cambridge, UK). Profiles for both modifications were pre-processed, normalized and filtered according to a stringent standardised quality control pipeline, as described previously 50 using the WateRmelon 68 package in R. We identified probes on the array within 1kb of differentially acetylated peaks (FDR < 0.05) using the Bioconductor package GenomicRanges 62 . A total of 1,659 of the 4,162 FDR significant differentially acetylated peaks were located within 1kb of at least one CpG probe on the array, with a total of 6,838 probes mapping to the 1kb neighbourhood of these 1,649 peaks. For each CpGpeak pair we correlated the log fold change in H3K27ac between AD cases and controls to the difference in DNA methylation or hydroxymethylation in a linear model controlling for the same covariates as in the differential acetylation analysis. Moreover, we compared DNA methylation and hydroxymethylation between probes in vicinity of AD hyper-and hypoacetylated peaks, as well as those in vicinity of all background peaks and the whole microarray background using two sample t-tests.

Accession codes.
Raw data has been deposited in GEO under accession number GSE102538.

Competing financial interests
The authors declare no competing financial interests.     divides the samples into two main groups: group 1 is composed mainly of controls, whereas group 2 contains more cases than controls. Interestingly, controls classified into group 2 are characterized by lower neuronal proportion estimates than those in group 1 (P = 0.004, ttest). The clustering defined by hyper-or hypoacetylated peaks is not significantly associated with sex (P > 0.05, chi-square test) or age at death (P > 0.05, t-test).