Alzheimer's disease: early alterations in brain DNA methylation at ANK1, BIN1, RHBDF2 and other loci.

We used a collection of 708 prospectively collected autopsied brains to assess the methylation state of the brain's DNA in relation to Alzheimer's disease (AD). We found that the level of methylation at 71 of the 415,848 interrogated CpGs was significantly associated with the burden of AD pathology, including CpGs in the ABCA7 and BIN1 regions, which harbor known AD susceptibility variants. We validated 11 of the differentially methylated regions in an independent set of 117 subjects. Furthermore, we functionally validated these CpG associations and identified the nearby genes whose RNA expression was altered in AD: ANK1, CDH23, DIP2A, RHBDF2, RPL13, SERPINF1 and SERPINF2. Our analyses suggest that these DNA methylation changes may have a role in the onset of AD given that we observed them in presymptomatic subjects and that six of the validated genes connect to a known AD susceptibility gene network.

Introduction

Evidence is emerging that DNA methylation levels at certain CpG dinucleotides can be both highly variable across individuals and stable over time within an individual[1,2]. This suggests that differences in DNA methylation in certain loci could be correlated with life experiences of a given individual, such as a disease risk factor or a diagnosis[3]. The most compelling evidence that the epigenome may influence AD comes from the manipulation of histone deacetylases (HDAC) in model systems of AD and in off-label treatment of AD patients with HDAC inhibitors[4,5]. To date, results of epigenomic studies of AD are sometimes conflicting and have not yet returned robust associations [6-9]. Further, there is a gradual increase in methylation at many sites throughout the genome with increasing age that has to be carefully considered when studying AD.[8,10-12] In this study, we present the results of a statistically rigorous gene discovery effort to identify regions of the genome that are differentially methylated in relation to the burden of AD neuropathology. We show that 71 discrete CpG dinucleotides in the human genome exhibit altered DNA methylation levels in relation to AD, that these changes are an early feature of AD, that the transcription of genes found in these differentially methylated regions is also independently associated with AD pathology, and that these differentially expressed genes connect to a previously reported genetically defined AD susceptibility network[13].

Description of subjects and data

Our data set consists of methylation measures at 415,848 discrete CpG dinucleotides in 708 subjects. These methylation profiles were generated using the Illumina HumanMethylation450 beadset and a sample of dorsolateral prefrontal cortex obtained from each individual. Since we dissected out the gray matter from each sample, we have profiled a piece of tissue composed primarily of different neuronal populations and other parenchymal cells such as glia. These subjects are part of the Religious Order Study (ROS) or the Memory and Aging Project (MAP), two prospective cohort studies of aging that include brain donation at the time of death. Since the subjects are cognitively non-impaired at study entry, we have studied a random selection of the older population. Over the course of the study, some subjects decline cognitively and display a range of amyloid pathology burden at the time of death, with 60.8% of subjects meeting a pathologic diagnosis of AD[14] (. To technically validate the nature of our data, we compared our Illumina-derived data to genome-wide DNA methylation sequence data generated from the same brain DNA samples in four of the subjects (two non-impaired and two AD subjects): in these four subjects, we see a very strong correlation (mean r= 0.97) between the estimated levels of methylation generated by the two technologies, consistent with prior reports[15].

Interestingly, when examining the nature of human cortical methylation profiles across our subject population, we note that the mean Pearson correlation of methylation levels for all possible subject pairs is 0.98 (), suggesting that that majority of CpG sites do not show significant interindividual variation in methylation levels despite the very different life course of each of these older subjects. As expected, we see many more differences in DNA methylation profiles between our cortical samples and lymphoblastic cell lines from HapMap individuals that were profiled for assessments of data quality in our experiment ().

Discovery study to identify differentially methylated chromosomal regions

Our analytic strategy involves three stages, which are illustrated in : Stage 1 is our DNA methylation screen for chromosomal regions in which methylation levels correlate with AD pathology. Details of the analytic model are presented in the Online Methods section. It is followed by Stage 2 in which we replicate the significantly associated CpGs from Stage 1 in an independent set of subjects. In Stage 3, we attempt to functionally validate the role of the differentially methylated regions that are replicated in Stage 2 using mRNA obtained from AD and non-AD subjects. This strategy accomplishes 2 goals: (a) further confirms the role of a given differentially methylated region by showing that a meaningful biological effect (transcriptional change) relates to the disease and (b) helps to narrow down which of the genes near the differentially methylated CpGs are differentially expressed and may therefore be the target gene(s) in a given region.

In the primary analysis of our cortical methylation profiles (Stage 1), we identified autosomal CpGs whose level of methylation correlates with the burden of neuritic amyloid plaques (NP), a key quantitative measure of Alzheimer's disease neuropathology. NP burden better captures the state of the brain of a deceased subject since cognitively intact individuals display a range of NP pathology, some of which meet neuropathologic criteria for a diagnosis of AD[16,17]. and contain the results of the primary genome-wide analysis: 137 CpGs are associated with the burden of NP pathology at a p<1.20×10−7. This threshold of significance accounts for the testing of all 415,848 tested CpGs by imposing a Bonferroni correction on a standard p<0.05. Since the exact number of functionally independent units of methylation in the genome is currently unknown, we have chosen this simple but conservative strategy to account for the testing of multiple hypotheses and correct for the testing of each CpG that was measured. Since the proportion of neurons found in each sample was not related to AD (p=0.08), we did not include this as a term in the primary analysis. However, to focus only on the most conservatively associated CpGs, we performed a secondary analysis that includes the variable that captures the proportion of neurons as well as surrogate variables that capture structure in the methylation data that do not correlate with known confounders and may capture cryptic technical or other artifacts. Of the 137 CpGs discovered in the primary analysis, 71 CpGs remain significant in the more conservative secondary analysis (). Some of these 71 CpGs are found in the same chromosomal segment and are highly correlated in their level of methylation. Altogether, the 71 CpGs are found in 60 discrete differentially methylated regions distributed throughout the genome (): in 8 of these 60 regions, up to three neighboring CpGs with correlated levels of methylation emerge as significant in our analysis and probably capture the same effect.

Individually, any one of the significantly associated CpGs (, ) has a modest effect on the brain's NP burden: on average, each the 71 CpGs explains 5.0% (range 3.7-9.7%) of the variance in NP burden. However, this is greater than the proportion of variance explained by genetic variants associated with AD, with the exception of APOE. For example, in our subjects, the well-validated CR1 susceptibility allele explains just 1% of variance in NP burden[18], and all known AD variants and APOE ε4 account for 13.9% of the variance in NP burden. If we consider all 71 CpGs in one comprehensive model, they explain 28.7% of the variance in NP burden, suggesting that methylation levels of certain genomic regions is correlated and that cortical DNA methylation of a large number of discrete regions is strongly correlated with a key measure of AD neuropathology.

Notably, two of the 71 significantly associated CpGs () are found in loci that harbor known AD susceptibility alleles: cg22883290 in the BIN1 locus (beta=4.44, p=9.00×10−8) and cg02308560 in the ABCA7 locus (beta=3.62, p=2.45×10−12)[19-22]. cg22883290, is located 5 kb from the 5’ end of the BIN1 gene and 92 kb from the index SNP, rs744373, that best captures the genetic association to AD in this region ()[19]. The susceptibility variant rs744373 is moderately associated with the level of methylation at cg22883290 (p=0.0003). However, the CpG association with AD pathology is not driven by the variant: adjusting for rs744373 does not meaningfully change the effect size of the CpG association to NP burden (model with rs744373 as a covariate: beta=4.37, p=4.91×10−7). Within our dataset of modest size, rs744373 is not associated with AD susceptibility, and we therefore cannot formally test for mediation of the SNP's association to disease by CpG methylation. In the case of ABCA7, the index SNP (rs3764650) is associated with NP burden[23] but has no association (p=0.07) with the level of methylation at cg02308560 which is 25 kb away, so, in both of these regions, SNPs and CpGs appear to have independent effects on AD susceptibility. Overall, risk of AD may therefore be affected by different sources of genomic variation (genetic and epigenetic) that have independent effects on the disease process.

To facilitate the interpretation of our results, we performed a secondary analysis correlating the level of methylation at these 71 CpGs with a post-mortem, neuropathologic diagnosis of AD. 22 of the NP-associated CpGs are also associated with a diagnosis of AD at a genome-wide level of significance ( and ), and all of the CpGs associated with NP burden display at least some evidence of association (p<0.001) with AD. This is not surprising since NP burden is one criterion for a neuropathologic AD diagnosis. We note an interesting polarization in the direction of these associations: 82% of the differentially methylated regions are more methylated in subjects with a diagnosis of AD. As noted above, the increased level of methylation in relation to AD at any one associated probe is modest ().

Validation of the associated CpG in an independent sample set

To further assess the robustness of our results, we evaluated the 71 significantly associated CpGs in an independent collection of 117 subjects with a different quantitative measure of AD pathology (Braak staging)[24]. These subjects were profiled in a sample of frontal cortex using the same Illumina Humanmethylation 450 platform; demographic details are presented in . We imposed a Bonferroni correction in this analysis and find that 12 CpGs are significant in this analysis (); since two of these CpGs are found in the same differentially methylated region near RHBDF2 (), eleven of the differentially methylated regions from the discovery study are validated (). Thus, despite the use of a different but related measure of AD pathology and a much smaller sample size, we see robust replication of our discovery screen's results. In addition to the significant CpGs, many other CpGs display suggestive evidence of association: when evaluating the entire set of 71 CpG, we see that the effect size of most of these CpGs is consistent across the two datasets (), suggesting that most of these CpGs will be validated as larger sample sizes are profiled. Both the BIN1 (p=0.0067) and the ABCA7 (p=0.011) CpGs display suggestive evidence of replication ().

Cognitively non-impaired subjects display the same alterations in methylation

To begin to explore the question of whether the increased level of DNA methylation in the associated regions is a cause or an effect of the neurodegenerative process of AD, we limited the NP analysis to those subjects who were deemed to be cognitively non-impaired at the time of death (no AD and no mild cognitive impairment). As has been well documented in neuropathological and imaging studies[25,26], a large fraction of non-impaired, older individuals demonstrate accumulation of amyloid pathology that is asymptomatic. Within the subset of non-impaired subjects, the p value for the CpG associations is diminished given the reduced sample size (n=237), but the beta values, which capture the magnitude of the association's effect, are not significantly different from the beta values calculated from the entire sample collection (). This suggests that the altered DNA methylation that we have identified in our discovery study is an early feature of AD pathology and occurs in the presymptomatic stage of the disease. These DNA methylation changes are therefore not secondary to the later stages of the dementing process. The question of whether altered DNA methylation contributes to the pathologic process or is an early epiphenomenon of the neurodegenerative process remains open.

Distribution of associated CpGs among different chromatin states

To better understand the functional consequences of the associated CpGs, we interpreted our results in relation to a chromatin map of the dorsolateral prefrontal cortex, generated in collaboration with the Epigenomics Roadmap team (http://www.roadmapepigenomics.org). It is derived from two MAP subjects who were cognitively non-impaired at the time of death and had minimal AD-associated pathology on post-mortem examination. Using histone modification profiles and established methods[27], each 200 bp segment of the genome is annotated as being in one of 11 chromatin states () that capture the transcriptional states and putative regulatory elements found in this tissue. Using this reference map, we see that at least some of the 71 associated CpGs are found in every chromatin state but that there is an enrichment of associated CpGs in regions predicted to be weak enhancers (p=0.0098) or to be in a weakly transcribed chromatin state (p=0.028), (). Further, we see a strong under-representation in regions displaying a strong promoter profile in the reference chromatin map (p=8×10−4). These data suggest that the chromatin architecture of strong promoters that drive fundamental cellular processes of neurons and glia in the healthy brain may not be strongly altered by AD. Rather, methylation changes appear to primarily affect genomic regions that are weakly transcribed or inactive in the healthy older brain. There are no enrichments noted in different genic features or in different structures defined in relation to CpG islands ().

Functional validation of the CpG associations

Focusing on the 12 CpGs that have been validated in the replication stage, we evaluated their role in AD by assessing the level of expression of genes found in the vicinity (±50 kb) of these DMRs in an independent set of 202 AD and 197 non-AD individuals assembled by the Mayo Clinic (see for demographic characteristics) that have RNA data from the temporal cortex. We find that the level of expression of 8 of the 21 selected genes have significant associations (p<0.0023) with AD in these data (): ANK1, CDH23, DIP2A, RHBDF2, RPL13, RNF34, SERPINF1 and SERPINF2.

Integrating our results with known AD genes

To further evaluate the role of these eight genes in relation to well-validated AD genes, we used the DAPPLE algorithm to evaluate the connectivity of these genes with the network of known AD susceptibility genes. We have previously used this method that requires co-expression of interacting protein pairs and adjusts for gene size, and we reported the existence of an AD susceptibility network derived from protein:protein interaction data[13]. Here, we use an updated model that includes the latest results from genome-wide association studies and the studies of rare variation. First, we find that the network of susceptibility genes from genome-wide association studies and mendelian AD genes is significant both in terms of direct connectivity (p=0.0072) and indirect connectivity (proportion of susceptibility genes sharing a common interactor, p=0.037)(). We then repeated the analysis after adding the eight genes found in the validated differentially methylated regions that also display altered RNA expression in AD. As seen in , several of the differentially expressed genes found in the differentially methylated regions - ANK1, DIP2A, RHBDF2, RPL13, SERPINF1 and SERPINF2 - connect to the AD susceptibility network derived from genetic studies. The direct (p=0.0072) and indirect (p=0.042) network connectivity remain significant in the iteration of the network analysis that includes the eight genes with altered RNA expression levels.

Discussion

Overall, while our study has certain limitations due to the Illumina platform - such as surveying only a fraction of the human genome's CpGs and the array's inability to distinguish two closely related chromatin marks (DNA methylation and DNA hydroxymethylation) we nonetheless report several replicated, functionally validated associations between altered DNA methylation and the presymptomatic accumulation of AD pathology. These changes do not appear to be part of a generalized, genome-wide process: specific differentially methylated regions are targeted in AD and are unlikely to be found in genes that are actively transcribed in the healthy older brain. Instead, the associated regions are more likely to be in a poorly transcribed chromatin conformation in the older healthy brain (). Altered DNA methylation () and enhanced mRNA expression (p=1.09×10−4) of the BIN1 gene in AD in the Mayo clinic dataset links our study to a well-validated AD susceptibility locus and a recent report of enhanced BIN1 expression in AD.[28] The BIN1 cg22883290 association is significant in our discovery study, and, while the results are suggestive (p=0.0067) in the small replication analysis, the direction of the effect is consistent (), suggesting that the association is likely to validate with additional subjects. Our results therefore refine our understanding of the BIN1 locus and suggest that different types of genomic variation (SNP and CpG) can have independent effects that integrate on the expression of BIN1 and influence AD susceptibility. Similarly, cg02308560 in the ABCA7 locus is associated to AD pathology, and this association is independent of the susceptibility allele found in its vicinity. These two loci illustrate the point that, while genetic variation can drive differences in DNA methylation for certain CpGs[3], the associations that we report are not driven by genetic associations with AD pathology: we have also recently completed a genome-wide SNP association study with the same trait in the same subjects and found no significant genetic associations[23]. Thus, our CpG associations are not driven by SNP associations.

The colocation of genetic susceptibility with CpG associations to AD pathology, along with the presence of the CpG associations in cognitively non-impaired subjects () and the connection of six of the differentially methylated genes (Figure 4) to an existing AD susceptibility network (), suggest that DNA methylation changes play a role in the onset of AD. However, as with genetic studies, our epigenome-wide scan only reports associations with a trait. Thus, we cannot state that the observed changes in methylation are causal: given the plasticity of the epigenome, it is possible that these changes are an early consequence of AD pathology. Further experimental work needs to be conducted to resolve this important question.

Figure 4

Genes identified in our DNA methylation screen connect to a network of known AD susceptibility genes

Using protein:protein interaction data, the DAPPLE algorithm evaluated the extent of connectivity among known AD genes (susceptibility and Mendelian genes) and the eight genes found in DMRs that are also differentially expressed relative to AD. The figure displays the results of an analysis allowing for one common interactor protein that is not known to be associated with AD. For example, RHBDF2 is displayed at the top of the figure in green and connects to PTK2B, a protein tyrosine kinase genetically associated with AD susceptibility which has a central role in this network. Interestingly, SERPINF1 and SERPINF2 connect to different elements of the amyloid component of the network (bottom left). Further, DIP2A connects the recently described PLD3 gene that has a rare AD susceptibility allele to SORL1, a gene with a common AD susceptibility allele, that connects to the amyloid precursor protein (APP). These interconnections are consistent with the reported effects of both PLD3 and SORL1 on amyloid biology and implicate DIP2A in the same process. Alternative figures presenting the network with all interacting proteins listed and the result of the network analysis with only the genetically associated loci are found in . The colored nodes are the proteins encoded by genes implicated in AD (genetic and epigenomic associations); the colors have no meaning. The connecting proteins not known to be associated AD are shown in gray.

Looking at the AD network map, DIP2A connects directly to the known SORL1 susceptibility gene[29] and indirectly to PLD3, a recently reported AD susceptibility gene that is otherwise not connected to the AD susceptibility network. DIP2A may function as a cell surface receptor protein[30] and, given the putative role of PLD3[31] and SORL1 in amyloid processing, its relation to the burden of NP pathology may well be related to a direct effect on amyloid processing. SERPINF1 and SERPINF2 also connect to elements of the amyloid machinery. Interestingly, SERPINF1 mRNA expression is reduced in AD (), and its knockdown in an in vitro system leads to reduced neurite outgrowth[32], suggesting one potential effect of this gene.

On the other hand, ANK1 and RHBDF2 connect to PTK2B, an AD gene that is a key element of the signaling cascade involved in modulating the activation of microglia and infiltrating macrophages. Several other AD genes, such as CD33 and EPHA1, connect to this molecule as well. While little is known on the potential role of ANK1, the connection of RHBDF2 with PTK2B is consistent with the known role of this molecule in myeloid cells: it is necessary for the transport of TNFα converting enzyme (TACE, also called ADAM17), which releases TNFα from the cell surface[33]. Absence of RHBDF2 in mice impacts the normal release of TNFα from the cell surface[34] and impairs systemic immune responses to pathogens[35]. In vitro work also suggests that RHBDF2 may function in regulating the substrate specificity of the TACE/ADAM17 protease which functions in the release of TNFα but also of other proteins such as epidermal growth factor (EGF)[36] Its exact role in AD is not clear at this point, but our data suggests that RHBDF2 expression is increased in the context of AD (). Its connection to PTK2B further suggests that it may be involved in the role of microglia and infiltrating macrophages in the AD pathophysiological process. Consistent with this, adjusting for an estimate of the number of microglial cells using an RNA-based model seems to account for the AD-associated differences in RHBDF2 mRNA expression ().

Such cell type adjustment using surrogate markers for different cell types[37-39] are crude analyses but are helpful to begin to assign some of the transcriptional alterations to certain cell types. For example, using GFAP expression as surrogate marker for astrocytes, we see that adding a term for GFAP expression in our assessment of CDH23 RNA expression largely abrogates the association of CDH23 RNA expression with AD (). Since GFAP expression is enhanced with astrocyte activation[40,41], we cannot distinguish whether the alteration of CDH23 RNA expression (and presumably its altered DNA methylation) is caused by an increased number of astrocytes near neuritic plaques, the activation of these astrocytes, or a combination of both effects. Regardless of the exact mechanism, in the case of CDH23, our DNA methylation screen has uncovered a robust alteration in methylation that can now be dissected mechanistically.

The human cortex has a complex architecture that includes many different types of neurons, glia, and other cells such as microglia, peripheral immune cells, and endothelial cells from cortical capillaries. The changes in DNA methylation that we report most likely represent the altered methylation state of a subset of cells within our cortical sample since AD pathology accumulates over several decades and only a small number of cells will be affected at a given time. It is too early to confidently differentiate between three possibilities that could explain these modest but robust changes in methylation that occur in relation to AD pathology: (1) a fraction of the constituent cortical cells are changing, such as activated astrocytes in the vicinity of neuritic plaques that overexpress CDH23, (2) the relative proportion of the constituent cell populations of the cortex is changing as some populations such as neurons are lost, or (3) there is a modest influx of immune cells from the systemic circulation that alters the relative abundance of the different cortical cell populations. It is also likely that more than one of these or other, unsuspected processes may be at work.

Overall, the replication of our study's results by an independent study[24] and its functional validation with RNA data makes our strategy for investigating the brain's epigenome more broadly relevant to other epigenomic epidemiology studies. With clear effect sizes in hand, our results can be used to calibrate the design of future human studies in the brain or other organs. We have also made the important observation that these epigenomic changes are happening early in the pathologic process: while subjects display no cognitive impairment but have accumulated amyloid pathology. Going forward in the aging brain, we clearly need to more precisely map the alterations of chromatin structure that contribute to AD pathophysiology and to assess, using model systems, whether remodeling the epigenome is a fruitful goal for the development of AD therapies.

Online Methods

Subjects and genotypes

The analyses in this manuscript included deceased subjects from two large, prospectively followed cohorts maintained by investigators at Rush University medical Center in Chicago, IL: the Religious Orders Study (ROS) and the Memory and Aging Project (MAP). The ROS cohort, established in 1994, consists of more than 1,100 Catholic priests, nuns, and brothers from 40 groups in 12 states who were at least 55 years of age and free of known dementia at the time of enrollment. The MAP cohort, established in 1997, consists of more than 1600 men and women primarily from retirement facilities in the Chicago area who were at least 53 years of age and free of known dementia at the time of enrollment. All participants in ROS and MAP sign an informed consent agreeing to annual detailed clinical evaluations and cognitive tests, and the rate of follow-up exceeds 90%. Similarly, participants in both cohorts signed an Anatomical Gift Act donating their brains at the time of death. The overall autopsy rate exceeds 85%. As in previous manuscripts[17], we analyzed the ROS and MAP cohorts jointly since they were designed to be combined, are maintained by a single investigative team, and a large set of phenotypes collected are identical in both studies. All aspects of these studies were approved by the Institutional Review Boards of Rush University Medical Center and Partners Healthcare. More detailed information regarding the two cohorts can be found in previously published literature[42,43]. Genotypes were available from prior studies and were derived from Affymetrix GeneChip 6.0 or Illumina Omni1-Quad genotypes and imputation using the HapMap reference, as previously described.[23]

The replication DNA methylation analysis uses samples of prefrontal cortex (PFC) obtained from 117 individuals archived in the MRC London Neurodegenerative Disease Brain Bank (http://www.kcl.ac.uk/iop/depts/cn/research/MRC-London-Neurodegenerative-Diseases-Brain-Bank/MRC-London-Neurodegenerative-Diseases-Brain-Bank.aspx). All samples were dissected by trained specialists, snap-frozen and stored at -80°C. Genomic DNA was isolated from ~100mg of each dissected brain region using a standard phenol-chloroform extraction method and tested for degradation and purity prior to analysis.

Temporal cortex expression levels for the autopsied Mayo Clinic subjects were obtained as part of a recently published brain expression GWAS (eGWAS), where the methodology is described in detail (eGWAS)[44,45]. Briefly, expression levels of 24,526 transcripts were measured from the temporal cortex of autopsied brains from subjects with pathologic AD (n=202) and those with other brain pathologies (non-AD, n=197). Total RNA extraction and QC were done using the Ambion RNAqueous kit and Agilent 2100 Bioanalyzer, respectively, according to published methods. Whole genome DASL expression microarrays (Illumina, San Diego, CA) were used for the transcriptome measurements of RNA samples that were randomized across chips and plates using a stratified approach to ensure balance with respect to diagnosis, age, sex, and RNA Integrity Numbers (RIN). Raw probe-level expression data exported from Genome Studio software (Illumina Inc.) were preprocessed with background correction, variance stabilizing transformation, quantile normalization and probe filtering using the lumi package of BioConductor[46,47]. Preprocessed probe transcript levels were used in the downstream analysis.

Phenotypes

Our primary phenotype of interest in this manuscript was the burden of neuritic plaques, a quantitative measure of the amount of AD neuropathology in the brain at the time of death. Brain autopsies in ROS and MAP were performed across the US, as described in detail elsewhere[42,43] Bielschowsky silver stain was used to visualize neuritic plaques within tissue sections from five brain regions: the midfrontal, middle temporal, inferior parietal, and entorhinal cortices, and the hippocampal CA1 sector. As in prior publications[17] a quantitative composite score of neuritic plaque burden was then computed for each individual by dividing the subject's raw count in each of the 5 regions by the population standard deviation in that same region, and then taking the average of the standardized counts across the 5 regions. Because the distribution of these average standardized counts is skewed, we used the square-root transformed values in our statistical analyses.

To put the associations with neuritic plaque burden in the context of an AD diagnosis, we also assessed for associations with a neuropathologic diagnosis of AD, which is determined on post-mortem examination. Specifically, subjects were classified as having a pathologic diagnosis AD if they had intermediate or high likelihood of AD based on the National Institute on Aging (NIA)-Reagan criteria. The NIA-Reagan criteria, which integrates both the Consortium to Establish a Registry for Alzheimer's Disease (CERAD) estimates of neuritic plaque density and Braak staging of neurofibrillary tangle pathology, was implemented as reported[48]. These diagnoses are made by board-certified neuropathologists without access to the clinical data collected during the study. All neuropathologic data is collected in a blinded fashion (relative to clinical diagnosis) by the neuropathology staff.

Experimental Protocol for DNA Extraction from Post-mortem Brain

100 mg sections of frozen dorsolateral prefrontal cortex were obtained from each of 761 deceased subjects from the ROS and MAP studies based at the Rush Alzheimer's Disease Center. These sections were thawed on ice, and the gray matter was carefully dissected from the white matter. DNA extraction was performed using the Qiagen (cat: 51306) QIAamp DNA mini protocol. The Qubit 2.0 Fluorometer was used to quantitate the DNA. 16uL of DNA at a concentration of 50ng/uL as measured by PicoGreen, was used by the Broad Institute's Genomics Platform for data generation by the Illumina InfiniumHumanMethylation450 bead chip assay. The platform produces a data file by implementing the recommended procedures of the proprietary Illumina GenomeStudio software, which includes color channel normalization and background removal. All data generation was conducted by laboratory personnel who were blinded as to the clinical and neuropathological phenotypes of each subject

Subject and Probe Quality Control

For the initial quality check of the data, we used the detection p-value criteria recommended by Illumina. These p-values represent the quality of the probes compared to background noise. We selected probes, which have detection p-value < 0.01 for all samples to ensure the use of good quality probes. Using these criteria, we selected 470,913 out of a total of 485,577 tested probes for further analysis. Of these 470,913 probes, 460,045 are found in one of the autosomes.

However, not all probes are unique: some probes are predicted to cross-hybridize with the sex chromosomes[49] based on sequence alignment. Many of these probes showed strong association to gender in a recent report that recommends discarding probes in which 47/50 nucleotides match the sex chromosome sequence during sequence alignment using BLAT[49]. Specifically, the authors recommend discarding 29,233 probes from the Illumina InfiniumHumanMethylation450 bead chip assay that meet this sequence match criterion. We implemented this recommendation and removed the 29,233 probes from our probe list for downstream analysis. Further, as noted below, we adjusted for gender in our analytic model given the well-described influence of sex on methylation levels.

In addition to cross-reactive probes, a substantial fraction of CpG probes also overlap with known Single Nucleotide Polymorphic sites (SNPs) based on the 1000 genome database[49]. Methylation level of these CpGs could be affected by a subject's genotype. However, SNPs with very low minor allele frequencies should not have major effect on the methylation data given our sample size. Therefore, we removed CpGs in which a SNP with a minor allele frequency (MAF) ≥0.01 exists within 10 base pairs upstream or downstream of the CpG site, where single base extension occurs. Based on the analysis reported by Chen and colleagues[49], we found a total of 14,964 autosomal, polymorphic CpGs with 0.01 ≤ MAF (EU) ≤ 0.99 in subjects of European ancestry, such as our subjects. These CpGs were removed from consideration, leaving 415,848 autosomal CpGs for downstream analysis.

We also evaluated the “per subject” quality of the data after the initial CpG data cleaning. The first step in this component of our quality control pipeline was to remove subjects with poor quality data from further consideration. Specifically, we used Principal Component Analysis (PCA) to select subjects that are within ± 3 SD from mean of a PC for PC1, PC2 and PC3. PCA was performed using a random selection of 50,000 autosome probes for all of the samples. Using these criteria, we removed 13 subjects from a total of 761 subjects for which the Illumina assay was attempted. We excluded an additional 8 samples having for poor bisulfite conversion (BC) efficiency, which is defined as having at least 2 of the 10 BC control probes that fail to reach a value of 0.8. Hence, after our “per subject” quality check, a total 708 samples were selected for downstream analyses.

Normalization of the data

We observed a strong batch effect in our data due to the use of two different thermocyclers during data generation. This batch effect was also confirmed by our PCA analysis. We evaluated different approaches to normalize our data, including COMBAT[50] and independent component analysis (ICA). However, when applying these approaches to our data and comparing them with an approach in which we adjust for given confounding variables (such as batch number and age) that we find to be associated with principal components within our data, these approaches (COMBAT and ICA) are overly conservative in normalizing our data. In particular, ICA includes adjustments for many unknown variables that capture structure within a set of data and are not associated with known confounders. Thus, such ICA-derived variables could be related to elements of disease pathophysiology that is not yet understood.

ICA is a matrix factorization technique that separates a matrix into statistically independent, non-Gaussian factors using the R-package “fastICA” developed by Marchini and colleagues[51]. We applied ICA to the matrix of methylation beta values to infer number of statistically independent components k iteratively from 1 to 40. In each instance, all samples are assigned k surrogate variables representing k entries of the mixing matrix. We then performed a refinement step similar to that used in “iSVA”[52] where for each independent component i (1 <= i <=k) we created a sub-matrix of the beta values populated only by methylation probes that were significantly associated with component i. We then performed ICA on this sub-matrix and selected the independent component most correlated with component i, and used those refined components as the surrogate variables. To set an optimal k we explored the average variance explained in the methylation matrix by each value of k, and determined that 7 was a conservative number of components that explained an appreciable percentage of the variance. These 7 surrogate variables were then used in our evaluation of normalization methods.

In our comparison of the results obtained using the two different approaches (including known technical confounders such as batch vs. COMBAT and ICA normalization), we find that the top 70 CpG in the primary analysis remain significant regardless of the approach used and that, while the p-value of the other 67 CpGs that are significant in our primary analysis fluctuates below our threshold of significance, they retain strongly suggestive evidence of association to the neuritic amyloid trait (95% remain at a p<5×10−5). As a result, and given our strategy of validating the CpG analysis with two rounds of RNA analyses to minimize false positives (), we opted to use the batch variable approach to directly address the source of technical variation and avoid over-correcting our data for variables that are not demonstrated to be associated with technical or demographic confounders. Thus, we prefer to minimize the risk of false negative results and minimize false positives with two rounds of validation.

Accounting for differences in the proportion of cell types in our tissue

Since the proportion of neurons and other cell types found in the human cerebral cortex can change with AD, we evaluated a technique to account for the possible difference in the proportion of neuronal cells in our brain samples. As discussed in the main text, we used an R package[53] to quantify the proportion of NeuN+ cells (primarily neurons) in each brain sample using DNA methylation data. We used the data from the NeuN+ nuclei found in this package to create convex combinations of purified profiles using non-linear least squares. This yielded an estimate of the proportion of NeuN+ nuclei in each of our samples. However, the resulting measure is not significantly associated with a pathologic diagnosis of AD (p=0.08), and we therefore have not included it as a covariate in our primary analysis. However, we add it as a covariate in a secondary analysis that leads us to identify the most conservatively associated subset of 71 CpG.

Data analysis and statistical modeling

For our primary and secondary analyses, we used the β-values reported by the Illumina platform for each probe as the methylation level measurement for the targeted CG site in a given sample instead of M-values (logistically transformedβ-values)[52]. While M-values have certain favorable statistical properties relative to β-values, they are less biologically interpretable than β-values[54]. We therefore have opted to use β-values which range from 0 (no methylation) to 1 (100% methylation), and show good correlation to estimates of DNA methylation derived from whole genome bisulfite sequencing approaches (see below for details). Any missing β-value was imputed using a K-nearest neighbor algorithm for k=100.

In order to perform our methylome-wide association study (MWAS) and discover differentially methylated regions associated with neuritic plaque pathology, we used a linear model, adjusting for age at death, sex, study (ROS or MAP), experimental batch and bisulfite conversion efficiency. A logistic regression was used for the AD analysis, using the same covariates. As noted in the section above, we considered but elected not to include terms for (1) the proportion of neurons in each sample or (2) surrogate variables that displayed no correlation with available technical or demographic variables. To account for the testing of multiple hypotheses, we used a Bonferroni correction which yields a p<1.20 ×10−7 as the threshold for genome-wide significance given the 415,848 autosomal CpG probes tested in our analysis. For annotation of the CpG probes, we used the hg19 human reference genome.

The bisulfite conversion efficiency term is calculated using the bisulfite conversion control probes, based on Illumina guidelines. Ten CpG sites designated by Illumina as control sites (6 CpGs targeted by type I probes and 4 CpGs targeted by type-II probes), where we expect each CpG to be 100% methylated, are used to control for non-complete bisulfite conversion. The bisulfite conversion efficiency term used in the primary analysis is the median methylation estimate from the 10 control sites. The bisulfite conversion term is calculated by taking the median value of the probes that Illumina provides to estimate bisulfite conversion efficiency.

To assess whether changes in DNA methylation are an early feature of AD, we compared the results of our association analysis obtained from the full set of subjects (n=708) to those obtained from the subset of subjects who were cognitively non-impaired at the time of death (n=217) by assessing, for a given CpG, whether the beta estimate from the non-impaired subset was different from the beta estimate derived from the entire cohort. We used the fact that the non-impaired subjects are a subset of the entire cohort and tested, using a 1 sample test and the t distribution, whether the results from this subset, specifically the beta estimate taking into account the standard error, were different from those measures obtained from the entire population (i.e. the total cohort of 708 subjects) from which the subset was drawn.

To assess mediation by a third factor (as in the analyses presented in ), we compared models with and without that third factor. For example, to test if a mRNA expression mediates the association between a CpG and NP burden, we compared the beta estimate of the CpG from the model with only the CpG and our standard confounders to the beta estimate of the CpG from the model with CpG, our standard confounders and a term for the mRNA. If the beta estimate from the model with the mRNA term is >10% less that the model without the mRNA term we say that the mRNA mediates the effect between CpG and NP.

Bisulfite-sequencing data generation and analysis

For four individuals (two with AD and two non-impaired, each pair consists of a man and a woman), we used the same DNA sample profiled using the Illumina HumanMethylation Array to perform whole genome bisulfite sequencing. Genomic DNA was fragmented to 100-500 bp using a Covaris S2 sonicator. DNA fragments were end-repaired, A-tailed, and ligated with methylated paired-end adapters (purchased from ATDBio). Whole genome bisulfite sequencing (WGBS) was performed as previously described[55]. In short, Illumina genomic DNA adapters were added to the fragments, and the adapter-ligated fragments were size-selected prior to two rounds of 5-h bisulfite treatments using the Epitect Bisulfite kit (Quiagen). Libraries were then purified and run on the Illumina HiSeq2000 using a standard 36 base protocol[55].

The WGBS libraries were aligned using BSMap 2.7[56] to the hg19/GRCh37 reference assembly. Subsequently, CpG methylation calls were made using custom software[57], excluding duplicate, low quality reads as well as reads with more than 10% mismatches. Only CpGs covered by ≥5x reads were considered for further analysis.

Methylation profile from 4 samples (2 with AD; 2 non-AD) generated from Illumina450k as well as Bisulfite-sequencing, were compared by selecting a random 50,000 autosomal CpG sites that are present in both datasets for each individual. In each subject, the estimated level of methylation at each CpG measured by the two technologies was compared using a student t-test. The results of the comparison of the technical replicates (Illumina vs. WGBS) are as follows: subject # ROS20963578, r=0.975; subject # MAP5797875, r=0.969; subject # MAP50403446, r=0.971; and subject # ROS20214850, r=0.972. Across the 4 subjects, the mean r=0.972.

Chromatin State Map

In collaboration with the Broad Institute's Roadmap Epigenomic Mapping Center, chromatin immunoprecipitation using antibodies targeting six different chromatin marks (H3K36me3, H3K27me3, H3K4me1, H3K4me3, H3K9ac, and H3K9me3) was performed independently for two MAP subjects who were cognitively non-impaired at the time of death and had minimal evidence of AD, vascular and lewy body pathology on neuropathological examination. In both cases, samples were obtained from the dorsolateral prefrontal cortex that is also sampled in our DNA methylation scan. Chromatin was extracted from each of the cortical samples. Library construction and sequencing were performed as previously reported[58]. Chromatin state maps were then learned from the sequence data (available at http://www.roadmapepigenomics.org). The H3K36me3 and H3K27me3 data were from donor id 149 while the other data sets were from donor id 112. The data were dichotomized using the default settings of the “BinarizeBed” command of ChromHMM[59] except and applying the ‘-center’ option to the already signal extended bed files. The models were trained by applying the “LearnModel” command with default settings from ChromHMM. We selected a model with 11 distinct chromatin states as the optimal model for our tissue sample.

To assess whether certain chromatin states are enriched for associations in our analysis, we compared the distribution of the 137 associated CpGs across the 11 chromatin states to that of all 415,848 probes tested using a chi-square test.

RNA data & analysis

A detailed description of the RNA data from the Mayo Clinic samples has been reported elsewhere[44,45] Since the analysis here explored the relation of mRNA levels to an AD diagnosis, a multivariable linear regression model was implemented, adjusting. covariates to correct for technical or biological variables, including age at death, sex, PCR plate, RIN and (RIN-RINmean)2. For those analyses which corrected for the expression levels of genes that are specific for the main five cell types present in the central nervous system (CNS), the following probes were used as covariates: ENO2 for neurons (ILMN_1765796), GFAP for astrocytes (ILMN_1697176), CD68 for microglia (ILMN_2267914), OLIG2 for oligodendrocytes (ILMN_1727567) and CD34 for endothelial cells (ILMN_1732799). Some or all of these 5 expression levels were included to account for neuronal loss, gliosis and/or vascular tissue in the assessed brain regions, where indicated.

Pathway analysis

We constructed a protein-protein interaction (PPI) network using the online tool DAPPLE [60] in order to determine whether the genes identified in our DNA methylation study significantly interact with known AD associated proteins discovered in genetic studies. We compiled a list of genes associated with AD, including 25 late onset AD GWAS genes, 3 early onset AD associated genes, and TREM1 (Replogle et al., submitted). We produced a PPI network with a cutoff of 2 interacting binding degrees. DAPPLE creates direct and indirect networks of connected proteins using evidence of physical interaction from the InWeb database, which contains 169,810 high-confidence pair-wise interactions involving 12,793 proteins [61]. To assess the statistical significance of PPI networks, DAPPLE applies a within-degree node-label permutation strategy to build random networks that mimic the structure of the original network and evaluates four network connectivity parameters on these random networks to generate empirical distributions for comparison to the original network. Our genetically defined AD network was significantly connected based on its Direct Edge Count (p=0.0072) and Seed Common Interactors Degrees Mean (p=0.037). Once the nine genes that emerge from the RNA functional validation analysis are added to the list of genes to be considered in the analysis, the PPI network is expanded and includes six of these genes. In this iteration of the analysis, both the Direct Edeg Count (=0.0072) and Common Interactors Degrees Mean (p=0.042) measures remain significant.

Table 1

CpGs associated with amyloid burden: validated CpGs and CpGs in known AD loci

CpG	Chr	Position(bp)	Discovery Study						Replication Study		Genes within 50 kb of associated CpG
			EstimateNpburdenModel 1	P-valueNpburdenModel 1	EstimateNpburdenModel 3	P-valueNpburdenModel 3	EstimateAD	P-valueAD	EstimateBraakscore	P-valueBraakscore
cg11724984	12	121890864	3.02	4.76E-09	4.48	9.75E-11	16.7	1.27E-07	3.40	5.44E-06	RNF34,KDM2B
cg23968456	10	73521631	4.97	3.97E-10	4.79	8.97E-09	16.8	1.45E-05	1.62	8.27E-06	CDH23, C10orf105,C10orf54
cg15821544	1	43473840	3.52	1.17E-07	5.11	8.09E-11	18.9	1.46E-07	2.47	2.08E-04	SLC2A1,FLJ32224
cg16733298	16	19127132	2.75	5.24E-08	3.96	1.32E-09	12.9	9.76E-06	3.58	2.30E-04	COQ7,ITPRIPL2
cg22962123	7	27153605	1.7	1.12E-07	3.08	8.79E-13	10.8	5.56E-08	5.13	2.79E-04	HOXA1,HOTAIRM1,HOXA2,AK291164,HOXA3,AK311383,BC035889,HOXA4,LOC100133311,HOXA5,HOXA6,DQ655986,HOXA7,HOXA9,HOXA10
cg13076843	17	74475294	2.35	1.68E-09	2.29	5.99E-08	9.32	5.81E-07	3.02	2.99E-04	UBE2O,AANAT,RHBDF2,AX747521,CYGB,PRCD
cg25594100	7	4786943	3.15	2.54E-11	4.22	3.33E-13	14.9	1.83E-08	4.09	4.19E-04	FOXK1,KIAA0415,RADIL
cg00621289	21	47855916	3.5	6.48E-08	5.06	2.18E-11	17.5	3.98E-07	2.02	4.95E-04	PCNT,DIP2A
cg19803550	17	1637391	4.36	1.04E-08	4.48	1.60E-08	19	6.81E-07	1.51	5.61E-04	PRPF8,TLCD2,MIR22HG,AF070569,MIR22,WDR81,SERPINF2,SERPINF1,SMYD4
cg03169557	16	89598950	4.86	3.99E-10	4.88	3.36E-09	18.6	1.67E-06	1.52	6.05E-04	ANKRD11,SPG7,SNORD68,RPL13,CPNE7
cg05066959	8	41519308	2.69	7.13E-14	2.78	7.56E-13	11.2	4.88E-10	3.45	6.48E-04	AGPAT6,NKX6-3,JA429246,ANK1
cg05810363	17	74475270	2.95	3.68E-10	2.86	8.11E-09	11.3	3.36E-07	2.76	7.93E-04	UBE2O,AANAT,RHBDF2,AX747521,CYGB,PRCD
Known AD loci
cg22883290	2	127800646	4.41	3.73E-08	4.44	8.97E-08	15.9	1.83E-05	0.96	0.0067	BIN1
cg02308560	19	1071176	2.19	3.06E-08	3.62	2.45E-12	13.7	1.77E-08	3.60	0.011	CNN2,ABCA7,HMHA1,POLR2E,GPX4,SBNO2

Table legend: The Np analysis reports the results of a linear regression analysis relating DNA methylation level to the burden of neuritic plaque. “Model 1” refers to the primary analysis, and “Model 3” refers to the secondary that includes a variable for the estimate of neuronal cells and surrogate variables. The “AD” analysis reports the result of a logistic regression relating the level of methylation of a given CpG to a pathologic diagnosis of Ad. The threshold of genome-wide significance is p<1.12×10-7 after accounting for the testing of probes genome-wide in both Model 1 and Model 3. We have listed all genes found within a segment 50 kb upstream and downstream of the associated CpG. In the “Known AD loci”, we have listed those CpGs that are significant in our discovery analysis

Glossary: AD – Alzheimer's disease, Chr – chromosome, NP – neuritic plaque, patho AD – pathologic diagnosis of AD

Table 2

Analyses of transcriptional data within loci associated with amyloid burden

Gene	Mayo Clinic subjects AD diagnosis
	Beta	P value
CDH23	0.58	<0.0001
DIP2A	0.07	<0.0001
RPL13	−0.08	<0.0001
RHBDF2	0.22	0.0001
SERPINF1	−0.14	0.0004
SERPINF2	0.2	0.0008
ANK1	−0.17	0.0012
KIAA0145	0.05	0.0013
SMYD4	0.03	0.065
PRPF8	−0.04	0.13
SLC2A1	0.02	0.13
PCNT	0.03	0.15
COQ7	0.02	0.22
ANKRD11	−0.02	0.24
C10ORF54	0.07	0.25
FOXK1	0.02	0.28
UBE2O	−0.02	0.28
RNF34	−0.05	0.30
WDR81	−0.03	0.39
SPG7	0.04	0.54
CYGB	−0.02	0.62
AGPAT6	0	0.9

The “Beta” refers to effect size in the linear regression associating gene expression and AD diagnosis, with non-AD being the reference subject group. The significance threshold is p<0.0023 given 22 tests. Results are shown for mean expression level across available probes for each gene. Per probe results are presented in Supplementary Table 6.