Alzheimer’s disease-associated (hydroxy)methylomic changes in the brain and blood

1The neuropathological cascade of the world’s leading cause of dementia, late-onset Alzheimer’s disease (AD), is characterized by the progressive accumulation of extracellular amyloid plaques and intracellular neurofibrillary tangles, followed by neuronal cell death. The susceptibility to AD is determined by the complex interaction of genetic, environmental, and life-style factors, as well as epigenetic factors. Genetic research has been successful in identifying genetic variants modulating susceptibility to AD, including the first and strongest genetic risk factor for AD in the APOE gene. In addition to APOE, large-scale genome-wide association studies looking at AD have identified a number of independent common variants with a small-to-modest effect size. Besides genetics, recent studies have suggested an important role for epigenetic mechanisms in the etiology of AD, with reports of both global and gene-specific alterations in epigenetic modifications.

2Several types of epigenetic DNA modifications have been described, including DNA methylation (5-methylcytosine, 5mC) and DNA hydroxymethylation (5-hydroxymethylcytosine, 5hmC). While the best studied epigenetic DNA modification, 5mC, plays an important gene regulatory role in most tissues, 5hmC seems to have a different impact on gene expression and is particularly enriched in the brain, where it may play an important role in learning and memory. Unfortunately, conventional bisulfite (BS) conversion, a widely used procedure when quantifying DNA methylation, does not distinguish between 5mC and 5hmC. However, combining measurements from BS- and oxidative BS (oxBS)-converted DNA now allows for the quantification of both 5mC and 5hmC levels (Fig. 1).

3Where genetic factors can identify persons at risk for developing AD from birth, epigenetic markers may offer more dynamic views on trajectories of biological change and may therefore be able to offer an improved, chronological insight into the sequence of events at different stages of AD. As brain tissue cannot be readily sampled in living humans, blood may offer an alternative. Available research on the blood DNA methylome in relation to AD is limited and mainly focuses on the direct comparison of AD cases and healthy controls. Identifying disease-predicting biological profiles at pre-dementia stages of AD may provide improved precision in predicting onset of dementia and give potential treatments a better timeframe to successfully impede, or even halt disease progression.

4In the present study, we explored the association between AD and epigenetic dysregulation by quantifying 5mC and 5hmC, as well as unmodified cytosine (UC) proportions, at a single-site resolution in middle temporal gyrus (MTG) tissue obtained from AD patients (n = 45) and elderly, non-demented controls (n = 35; see Table 1 and the “Materials and methods” section for detailed demographics). This brain region was selected as the MTG is known as a site of early AD pathology, and differences in global levels of DNA methylation and hydroxymethylation have previously been reported in this brain region in AD. While informative on its own, the inclusion of UC measurements also allows us to better compare our findings with previous studies using conventional BS conversion, since UC is determined by subtracting the BS signal (5mC + 5hmC) from 1 (Fig. 1). Even though the effects will be opposite from directly using the BS signal, incorporating UC in our study represents a crucial legacy analysis that enables the comparison with previous studies solely relying on the BS signal. Moreover, mechanistically, as an example, the affinity of a transcription factor may be different in the presence of UC, 5mC, or 5hmC, implicating that differential levels of UC (in the absence of significantly different 5mC or 5hmC levels) may have direct functional implications on gene expression. We followed up the brain analysis exploring DNA methylation in whole blood in an independent cohort, including samples from AD-converters and non-converters at two time points, before (54 converters, 42 controls) and after (41 converters, 42 controls) conversion to clinical AD (see Table 2 and the “Materials and methods” section for detailed demographics). Blood DNA methylomic markers were measured using only BS-converted DNA, as 5hmC has a very low prevalence in blood.

5Site-specific 5mC, 5hmC, and UC levels were determined for the MTG using Illumina’s Infinium HumanMethylation450K microarray (HM 450K array) with BS and oxBS-converted DNA (Fig. 1; see Tables 1 and 2 for cohort demographics). An epigenome-wide association study (EWAS) was performed for each DNA modification to identify the association with AD. The adjusted linear models showed no signs of inflation (all lambda values were between 0.95 and 1.05; see Additional file 2: Figure S1 for QQ plots). None of the AD-associated CpG sites in the MTG passed false discovery rate (FDR) correction (Additional file 1: Tables S1–S3).

6A structural and functional genomic annotation enrichment analysis on the 1000 highest ranked sites indicated a significant enrichment of several CpG island features, gene features, and alternative transcription events. This included an enrichment of mainly gene body sites for the 5mC (fold enrichment = 1.42, p = 1.17E−10) and 5hmC (fold enrichment = 1.17, p = 3.64E−03) results and mainly intergenic sites for the UC (fold enrichment = 1.59, p = 1.67E−09) results (Additional file 2: Figure S8; Additional file 1: Table S7).

7A regional analysis, looking at the spatial correlation of adjacent modified positions, detected 1 differentially methylated region (DMR), 1 differentially hydroxymethylated region (DHR), and 11 differentially unmodified regions (DURs) that were associated with AD in the MTG (Table 3; Additional file 2: Figure S3). Analysis of MTG expression data of genes annotated to DMRs, DHRs, and DURs showed a significant negative correlation between a DUR associated with RHBDF2 and RHBDF2 RNA expression (ρ = -0.39, pFDR = 4.37E−03) (Additional file 1: Table S10). Of note, although the DHR residing in the transcription start site (TSS) of CHRNB1, of which all probes show hyperhydroxymethylation in the AD cases, did not correlate with CHRNB1 mRNA expression (ρ = − 0.09, pFDR > 0.05), a linear regression analysis of regressed MTG expression data of CHRNB1 showed a significant elevation of CHRNB1 mRNA levels in AD cases (estimate = 0.13, p = 1.37E−04) (Additional file 2: Figure S4). For a full transcriptomic investigation of the MTG cohort used in the present study, see the recent publication of Piras et al..

8Next, a gene regulatory network (GRN) analysis was performed with the unique genes annotated to the 1000 highest ranked probes. Because of different numbers of associated genes from each dataset, we obtained contextualized networks with varying number of interactions. The number of interactions in the contextualized GRNs representing the differential 5mC, 5hmC, and UC MTG states were 325, 398 and 244, respectively. Differential GRN analysis identified several candidate genes highly influential in the simulated transition from a diseased towards a healthy phenotype. Based on a score indicating for each gene, when changed, the number of other genes in the network that were predicted to show altered expression, IL6 (score = 55), SIAH1 (score = 78), and EGF (score = 55) were found to be the most influential in the 5mC, 5hmC, and UC networks, respectively (Additional file 1: Table S9).

9Since 5hmC is not enriched in the blood, only BS conversion was used to measure site-specific 5mC levels, also with the HM 450K array. A blood EWAS investigating the association between DNA methylation and conversion to AD was performed at baseline and at follow-up, leading to the identification of 3 differentially methylated positions at baseline and 266 at follow-up (Additional file 1: Tables S4–S6). No significant inflation was detected (Additional file 2: Figure S2; see the “Materials and methods” section for details).

10Genomic annotation enrichment analysis of the top sites in blood showed enrichment of mainly intergenic sites (fold enrichment = 1.32, p = 5.80E−04) at baseline and proximal promoters (fold enrichment = 0.79, p = 1.60E−04) at follow-up (Additional file 2: Figure S9; Additional file 1: Table S8).

11The regional analysis found 15 and 21 DMRs associated with conversion to AD at baseline and follow-up, respectively (Table 4; Additional file 2: Figure S5).

12GRNs representing the blood baseline and follow-up states contained 475 and 277 interactions, respectively. Differential GRN analysis identified WNT3A (score = 50) as the most influential gene in the baseline network, and SHH (score = 33) in the follow-up network (Additional file 1: Table S9).

13Only 1 blood DMR, close to GLIPR1L2, showed hypermethylation in relation to AD conversion at both the baseline (+2.72%, pŠidák = 1.40E−04) and follow-up (+ 1.34%, pŠidák = 6.94E−06) time points. Extracting the probes located in this blood GLIPR1L2 DMR from the MTG EWAS for comparison showed, in AD cases, lower UC levels (9/10 probes with negative log2 fold change [logFC]), mixed changes for 5mC (6/10 probes with positive logFC), and lower 5hmC levels for the probes that passed the detection threshold (2/2 probes with negative logFC). Even though the UC observations in the MTG are in line with the blood findings, only for one UC probe (cg07311024) the change was nominally significant (logFC = − 0.01, p = 3.88E−02). A targeted linear regression analysis of the regressed MTG expression data of GLIPR1L2 showed a significant decrease in AD cases (estimate = -0.10, p = 3.12E−04) (Additional file 2: Figure S6).

14Interestingly, close to the TSS of OXT, we observed a DMR which was detected both in the MTG (− 3.76%, pŠidák = 1.07E−06), as well as in the blood dataset (at baseline, + 3.43%, pŠidák = 7.14E−04) (see Additional file 2: Figure S7 for the probe positions of both OXT DMRs). MTG OXT methylation across Braak stages, as a proxy indicator of disease progression, is displayed in Fig. 2 and suggest OXT hypermethylation towards Braak 3-4 stages and OXT hypomethylation during later stages. Moreover, in the differential GRN analysis, OXT came forward as an influential gene. In case of the 5mC and 5hmC MTG states, a change in OXT was predicted to alter the expression of 39 and 54 other genes in the networks, respectively, and in the blood baseline state, OXT was predicted to alter 41 genes in the network (Additional file 2: Figures S10–S12; Additional file 1: Table S9).

15For the current study, we aimed to identify AD-related changes in epigenetic DNA modifications, comparing brain tissue from AD patients and age-matched controls. In addition, we explored DNA methylation in blood samples from AD-converters and non-converters, both at a preclinical stage and after conversion, identifying an AD-associated DMR in OXT in both the brain and blood datasets.

16The DHR identified in the MTG resided in the promotor of CHRNB1, which encodes acetylcholine receptor subunit beta and is important for cholinergic neurotransmission. In combination with the observed increased levels of CHRNB1 mRNA in the MTG, this potentially reflects a compensatory mechanism to maintain acetylcholine signaling in AD. Indeed, the acetylcholine-related pathway is known to be altered in AD and, as such, remains an important target in the development of novel treatment options. Previous epigenomic studies of AD using standard BS-conversion have found associations between AD and RHBDF2 methylation in multiple cortical regions. We replicated these findings; observing an AD-associated DUR in RHBDF2, which included the previously detected CpG sites (cg13076843, cg05810363, and cg12163800) and showed the same direction of effect as previously reported. For instance, using conventional bisulfite (BS) conversion, a 3.36% increase in DNA methylation level of cg05810363 has been observed across cortical regions in association with AD neuropathology. Interestingly, a negative correlation between UC levels within the RHBDF2 DUR and RHBDF2 mRNA expression was observed in the MTG. RHBDF2 is thought to be important for the release of tumor necrosis factor, a major inflammatory cytokine associated with neuroinflammation observed in AD. C3, another gene with an AD-associated DUR, encodes a central component of the complement system and mediates developmental synapse elimination by phagocytic microglia. C3 has previously been implicated in mediating synaptic loss in the early stages of AD.

17The top DMR from the baseline blood analysis, showing hypermethylation in AD, is close to the LDLRAD4 gene. This gene has previously been associated with schizophrenia and blood pressure and is thought to suppress transforming growth factor (TGF)-β signaling. TGF-β is an inflammatory cytokine playing a role in cell survival and synaptic transmission, and various isoforms have been associated with AD. Additional baseline blood DMRs were close to TENM3, involved in neurite growth, SYMPK, involved in polyadenylation regulation of gene expression and which showed increased expression in AD, SLC44A4, associated with type 1 diabetes mellitus and human aging, ZMAT2, which had decreased expression in AD, ULK1, which may play a role in the autophagic degradation of amyloid beta (Aβ), and RUNX2, which links bone health and cognitive function and anxiety-like behavior. The DMR that was found both at baseline and follow-up is associated with GLIPR1L2. GLIPR1L2 also showed decreased expression in the MTG. The function of this gene is not well known, but it may play a role in tumor suppression and immune function. The top AD-associated blood DMR at follow-up, showing hypomethylation, is located in GSDMD, which encodes a critical factor in pyroptosis; a form a cell death that may be triggered by Aβ. Other genes with a nearby AD-associated blood DMR at follow-up include KHDRBS2, previously identified in a genome-wide association interaction analysis in relation to AD, RARRES2, encoding an adipokine that has been linked to inflammation, obesity, diabetes, and cardiovascular diseases, and GNG7, for which Braak stage-associated differential methylation has been reported in cortical glial cells of AD patients.

18Taken together, the observation of epigenetic modifications in several inflammation-associated genes in both brain and blood aligns with the amyloid cascade-inflammatory hypothesis of AD. These findings could reflect either downstream effects resulting from the inflammatory activation seen in AD, or, particularly in the brain, reflect mediating effects of DNA modifications on inflammation as a causative factor. Exploring the exact nature of the AD-associated epigenetic modifications in inflammation-associated genes and the potential for blood biomarkers is thus a pivotal aim for future studies.

19Strikingly, our methylomic profiling in MTG and whole blood both led to the identification of a common DMR associated with AD, close to the transcription start site of OXT. Our design allowed for the disentanglement of specific 5mC and 5hmC signals in the MTG, which, in the case of OXT, suggests they change in opposite directions in relation to AD. The detection of a DMR near OXT is in line with a recent report of a nearly identical AD-associated OXT DMR (containing 1 additional probe) in the superior temporal gyrus (STG). This area is located directly above the MTG. Furthermore, using GRN analysis addressing the overlap between the top influential genes in the networks and genes with significant differentially modified regions, we observed OXT to consistently appear as one of the most influential genes in both brain and blood GRNs. OXT encodes oxytocin, a neuropeptide involved in the neuromodulation of social behavior, stress regulation, and associative learning. Interestingly, the functional impact of OXT promoter methylation at the same genomic locus has been recently shown. It was linked to several measures of sociability, superior temporal sulcus activity during social cognition tasks, as well as fusiform gyrus gray matter volume, a brain region closely related to the MTG.

20The paraventricular nucleus and supraoptic nucleus are thought to be the main sites of central oxytocin production, areas which reportedly undergo cell loss during AD. The remaining neurons are thought to undergo a, potentially compensatory, hypertrophy. One might hypothesize that this activation could initially lead to higher than normal oxytocin levels, before synthesis collapses during the final stages of AD. Interestingly, enhanced levels of hippocampal oxytocin have been associated with memory impairment, and AD-associated elevations of oxytocin have been reported in the hippocampus and temporal cortex. There is also limited evidence oxytocin is reduced in cerebrospinal fluid of manifest AD. Additionally, it has been suggested that co-damage to the locus coeruleus and hypothalamic nuclei could happen early during AD pathogenesis, substantiating the hypothesis that oxytocin could serve as an early diagnostic biomarker for AD.

21In line with an AD-related increase in temporal cortex oxytocin levels, all ten CpG sites within the MTG OXT DMR showed decreased levels of methylation in AD cases. Conversely, we observed OXT hypermethylation in the DNA from peripheral blood of participants who progressed to dementia. The OXT blood DMR was not observed after conversion. Research has shown that independent mechanisms may be involved in peripheral and central regulation of OXT expression, supporting this apparent discrepancy observed in blood and brain. Alternatively, these observations suggest there may be a temporal change in OXT methylation during AD progression. Looking at MTG OXT methylation across Braak stages appears to support the observation of OXT hypermethylation at earlier stages, as also seen in the blood, and OXT hypomethylation at more advanced stages. Notably, it has recently been reported that oxytocin administration was able to improve social cognition and behavior in frontotemporal dementia patients, illustrating the complex modulatory function of oxytocin in different brain regions and its potential use in the treatment of certain manifestations of dementia. Whether oxytocin represents a suitable therapeutic agent for AD remains to be elucidated.

22Even though we detect several targets relevant in light of AD, a general lack of overlap between the different analyses presented here might be noted, an observation which is true for EWAS and epigenetics studies in AD in general. Others have discussed a myriad of possible reasons for discrepancies between studies, such as methodological differences, differences in tissue type and processing, study designs, and samples sizes. In view of this, the detection of a common OXT DMR in two completely independent cohorts and two different types of tissue, further supported by a recent similar EWAS on the STG, makes it an even more promising target for future studies. However, the differences in direction of change and the OXT methylation pattern observed over Braak stages indicates these epigenetic changes should be further studied in a longitudinal fashion to establish a clear relationship with AD neuropathology, as well as clinical manifestations of AD.

23Given the detection of several regions of interest, it should be noted that the lack of positions significantly associated with AD in the MTG after FDR correction may be the result of a limited sample size. Genome-wide site-specific AD-related epigenetic changes should thus be further investigated using studies with larger sample sizes or meta-analyses. Alternatively, future studies may focus on candidate genes identified in the present work, such as OXT.

24Our novel approach confirms some previous epigenetic findings identified in the central nervous system, including RHBDF2, as well as revealed novel targets, such as in CHRNB1, involving dysregulated DNA hydroxymethylation. Furthermore, the nearly identical OXT DMRs found in both the blood and brain suggest a systemic epigenetic dysregulation in AD involving OXT. The detection of the OXT DMR at pre-dementia stages suggests its potential relevance as a novel biomarker and may offer new treatment strategies to be explored in future studies.

25Informed consent was obtained from all human participants. This includes donors of the Banner Sun Health Research Institute (BSHRI) Brain and Body Donation Program (BBDP), who signed an Institutional Review Board-approved informed consent form, including specific consent to the use of donated tissue for future research. The German Study on Ageing, Cognition and Dementia in Primary Care Patients (AgeCoDe) study protocol was approved by the local ethics committees at the University of Bonn (Bonn, Germany), the University of Hamburg (Hamburg, Germany), the University of Duesseldorf (Duesseldorf, Germany), the University of Heidelberg/Mannheim (Mannheim, Germany), the University of Leipzig (Leipzig, Germany), and the Technical University of Munich (Munich, Germany).

26DNA from the MTG was obtained from 82 AD patients and neurologically normal control BBDP donors stored at the Brain and Tissue Bank of the BSHRI (Sun City, AZ, USA) (Table 1). The organization of the BBDP allows for fast tissue recovery after death, resulting in an average post-mortem interval of only 2.8 h for the included samples. Braak staging was carried out for AD neurofibrillary pathology. A consensus diagnosis of AD or non-demented control was reached by following National Institutes of Health AD Center criteria. Comorbidity with any other type of dementia, cerebrovascular disorders, mild cognitive impairment (MCI), and presence of non-microscopic infarcts was applied as exclusion criteria. Although this may limit the generalizability of the current study, these strict exclusion criteria were applied to enhance the detection of AD-specific dysregulation, not confounded by common comorbidities. Detailed information about the BBDP has been reported elsewhere.

27AgeCoDe is a prospective longitudinal study including 3327 non-demented individuals at baseline, initiated to investigate the early detection of MCI and dementia in primary care. Participants were randomly selected from the general practice registry in six German cities and cognition was assessed at approximately 18-month intervals and 10-month intervals after visit 7, for up to 11 years after baseline. For this study, whole blood DNA was obtained from a subsample of 99 individuals aged above 75 years from this AgeCoDe cohort (Table 2). Of these, 42 were converters: they had no dementia at baseline, had DNA samples available at baseline and follow-up (after ~ 4.5 years), and had sufficient information available for a diagnosis of AD dementia to be made at the 4.5-year follow-up. There were 44 control subjects, who had to adhere to the same criteria, except that they should have no signs of dementia at neither baseline, nor the 4.5-year follow-up, and all subsequent cognitive assessments up to 11 years after baseline. The remaining 13 participants had not yet converted at the 4.5-year follow-up (when blood was drawn), but were diagnosed during a later follow-up, up to a maximum of 11 years after baseline. These samples were grouped together with the other converters.

28The groups were matched for age, gender, and APOE genotype. The presence of dementia was assessed in all subjects with the Structured Interview for Diagnosis of Dementia of Alzheimer Type, Multi-infarct Dementia, and Dementia of Other Etiology based on the DSM-IV criteria. The dementia diagnosis in subjects who were not personally interviewed was based on the Global Deterioration Scale (≥ 4) and the Blessed Dementia Rating subscales. The etiological diagnosis of AD was based on the criteria of the National Institute of Neurological and Communicative Disorders and Stroke and the Alzheimer’s Disease and Related Disorders Association for probable AD and was only assigned in case of sufficient information provided. All final diagnoses were a consensus between the interviewer and an experienced geriatrician or geriatric psychiatrist. More detailed information about the AgeCoDe cohort has been published previously.

29For the BBDP samples, the TrueMethylTM 24 Kit version 2.0 by CEGXTM (Cambridge Epigenetix Limited, Cambridge, UK) was used for BS and oxBS conversion of genomic DNA (gDNA) extracted from frozen MTG tissue. All laboratory procedures were performed at GenomeScan (GenomeScan B.V., Leiden, the Netherlands), without knowledge of the phenotypic characteristics of the samples and according to the manufacturer’s instructions. Prior to conversion, high molecular weight (HMW) gDNA was quantified using a PicoGreen assay (Invitrogen, Carlsbad, CA, USA), and gel-electrophoresis was performed to assess gDNA quality. All samples were of sufficient quantity and quality. A volume of 1 μg HMW gDNA was used per sample, which, after purification and denaturation, was split into two samples that underwent either DNA oxidation (oxBS samples) or mock DNA oxidation (BS samples). Subsequently, all samples were BS-treated, and the yield of the samples was assessed by a Qubit ssDNA assay (Invitrogen). An additional quality control, using a restriction enzyme only able to cut unconverted cytosines, was performed for a qualitative assessment of 5hmC oxidation and BS conversion. From each BS/oxBS-treated DNA sample, 8 μL was amplified and hybridized on HM 450K arrays (Illumina, Inc., San Diego, CA, USA), and the Illumina iScan was used for imaging of the array. Sample preparation, hybridization, and washing steps for the Illumina Infinium Methylation Assay of the BeadChip arrays were performed according to the manufacturer’s protocol.

30For the AgeCoDe samples, gDNA was isolated from whole blood and DNA concentration and purity was determined using the NanoDrop ND1000 spectrophotometer (Thermo Fisher Scientific). All samples were of sufficient quantity and quality. Five hundred nanograms of gDNA was used for BS conversion, using a Qiagen EpiTect 96 Bisulfite Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. A total of 200 ng of BS converted DNA was analyzed using HM 450K arrays according to the manufacturer’s instructions. The Illumina iScan was used for imaging of the array.

31Total RNA extracted from frozen MTG, from matched samples as used for the epigenetic MTG analyses, was isolated with the RNeasy Mini Kit (Qiagen) starting with at least 60 mg of tissue. Raw expression data was obtained at the BSHRI, using the HumanHT-12 v4 BeadChip (Illumina).

32All computational and statistical analyses were performed using the statistical programming language R (version 3.3.2) and RStudio (version 1.0.136), unless otherwise specified. Raw IDAT files from the Illumina iScan were loaded into R using the minfi package (version 1.20.2). To confirm that the longitudinal samples were from the same donor a genetic fingerprinting test was performed based on the 65 SNP probes included on the HM 450K chip, as implemented in the ewastools package. Based on this test, 2 donors with mismatching samples were detected and excluded from the blood data. Next, the gender of the samples was predicted based on X chromosome methylation using the DNAmArray package (version 0.0.2), compared with the assumed gender, and mismatches were excluded (1 mismatched sample was excluded from the blood data). Cross-hybridizing probes and probes containing a common SNP in the sequence or within 10 bp of the sequence were removed. The “pfilter” function of the wateRmelon package (version 1.18.0) was used for probe filtering (6 969 and 1 437 probes were removed from the MTG and blood data, respectively). The remaining probe data was normalized using the dasen method, as implemented in the wateRmelon package. Probes on the X and Y chromosomes were excluded from further analyses.

33Following normalization, two sets of beta values, from the standard BS arrays (5mC + 5hmC) and from the oxBS arrays (5mC), were generated in case of the MTG. By subtracting oxBS beta values from the BS beta values (ΔβBS-oxBS) for each probe in each sample, 5hmC levels were calculated (Fig. 1). UC values were determined as 1-BS (1-βBS). It should be noted that other DNA demethylation intermediates, such as 5-formylcytosine (5fC) and 5-carboxylcytosine may be represented in the BS or UC levels, as it is currently unclear how these intermediates respond to oxBS conversion. However, these intermediates are present at very low levels and are not enriched in brain tissue like 5hmC is. In order to reduce noise and filter out non-hydroxymethylated sites, outliers deviating more than ± 2SD from the probe mean in the 5hmC dataset were determined and set to the mean ± 2SD first, and subsequently, a threshold of zero was applied to the mean of individual probes (218,009 5hmC values were excluded). Boxplots and density plots of raw and normalized beta values per sample were inspected for clear outliers (2 MTG samples were excluded due to clear deviation from the other samples; data not shown). After data processing, 80 MTG and 96 blood samples remained, with 396,600 remaining probes for MTG 5mC and UC, 178,591 5hmC MTG probes, and 402,480 remaining probes in the blood datasets. The case-control analysis of the blood baseline data included all 96 samples (54 converters, 42 controls), while follow-up data included 83 samples, including the 41 converters that had already converted to AD at the 4.5-year follow-up and excluding those that had converted later. All individuals in the follow-up analysis were also included in the baseline analysis.

34An initial model with beta values as outcome, AD diagnosis/conversion as predictor, and age and gender as covariates was used for a surrogate variable (SV) analysis with the sva package (version 3.22.0). The first 5 SVs of this analysis were added to the model to adjust for unobserved confounders, including potential batch effects and differences in cell type composition. As the addition of SVs still resulted in inflation of the regression statistics (lambda = 1.43) of the blood follow-up analysis, and none of the SVs strongly correlated with the HM 450K chip IDs (which was the case for the other analyses), the chip IDs were also added to the model for this analysis. This successfully eliminated the inflation (lambda = 1.00).

35Linear regression was performed using the limma package (version 3.30.11) to test the association between the beta values and AD diagnosis/conversion. Test statistics were adjusted for bias and inflation with the bacon package (version 1.2.0). An FDR correction for multiple testing was applied to the p values to identify differentially (hydroxy)methylated and unmodified positions (probes with pFDR < 0.05). Individual probes were annotated using Illumina UCSC annotation.

36To examine the distribution of 5mC, 5hmC, and UC levels across genomic regions, we annotated the 1000 highest ranking probes (Additional file 1: Tables S2–S7) using ENCODE annotation data, as described by Slieker et al.. Fisher’s exact test was used to assess enrichment in specific genomic regions.

37To identify differentially (hydroxy)methylated and unmodified regions (DHRs/DMRs/DURs), spatial correlations between p values of the association analysis were determined using comb-p with a seeding p value of 0.01 and a window size of 1000 bp. Obtained p values were Stouffer-Liptak-Kechris corrected for adjacent p values and were subsequently corrected for multiple testing using the Šidák correction. Of the regions detected by comb-p, only those containing at least 3 CpGs and having a pŠidák < 0.05 were accepted as differentially modified regions.

38GRNs have been extensively used to achieve deeper understanding of disease related mechanisms. Different topological characteristics of these networks, such as connectivity of nodes or gene-gene interaction tendency in cell/tissue specific contexts, have been used to predict disease-related genes. Here, we have employed an in-house developed differential GRNs inference approach, which relies on gene expression data to infer GRNs specific to a given gene expression program. The initial set of interactions among the genes of interest was compiled from literature-based database ARIADNE and consists of interactions belonging to the categories of “Direct Regulation,” “Expression,” and “Promoter Binding.” The obtained set of interactions is not context-specific as they are reported to happen in different cell/tissue types and organisms. To obtain context-specific networks from the literature interaction maps, the pruning of interactions incompatible with the gene expression state was carried out, which resulted in contextualized networks compatible with the given gene expression state of the system. As a differential expression setting was used here, we obtained two contextualized GRNs for each state, representing the different network topology of diseased and healthy phenotype. The differential network topology helps us in identifying the set of genes that are regulated by different transcription factors in both networks. These genes formulate an ideal set of candidate perturbagens, as to change their expression state we have to perturb them individually. The obtained contextualized networks were used to identify genes in the common elementary circuits (positive and negative circuits) that can also serve as a set of candidate genes for perturbation. Genes in elementary circuits have been reported to play a crucial role in maintaining network stability and are considered as a necessary condition for a network to have an attractive cycle. In this regard, genes present in the common elementary circuits are considered to be the backbone of the network and any perturbations in the expression levels of these genes might lead the system to deviate from the normal steady state of the system, which can be described as a transition from healthy to a diseased state. Once we obtained a set of optimal perturbation candidates, we performed single-gene perturbation simulations to see the effect of change in expression of a single gene on all the other genes in the GRN. This measure tells us about the influential capability of the selected gene in the network; the higher the number of downstream genes being affected by perturbing a candidate gene, the more crucial is its role in the regulation of other genes in the GRN.

39Positions from the AD association analyses were ranked based on a combined p value and log2 fold change ranking score. The GRN analysis was then conducted separately for the genes annotated to the 1000 highest ranked sites in the MTG (5mC, 5hmC, and UC separately) and blood (baseline and follow-up separately) (Additional file 1: Tables S2–S7). Closest UCSC TSS annotation was used to obtain unique genes. After applying the differential GRN analysis on the contextualized networks, we ranked the key candidate genes based on their scores. This score represents the number of genes whose gene expression is changed (shifted from diseased towards the healthy phenotype) upon perturbation of the candidate gene.

40Raw RNA expression data was exported from Illumina’s GenomeStudio (version 2011.1) with the Expression Module (v1.9.0) for further analysis in R. Of the 80 subjects used for the epigenetic analyses, 1 case was not included on the expression array, and 3 additional cases were excluded after quality control of the data, due to extreme outlying values or failed reads, leaving 76 subjects for further analyses. Data was quantile-quantile normalized. Using the same model as for the regression analysis, the sva package was used to determine SVs for the epigenetic and expression datasets. The effects of age, gender, and 5 SVs were regressed out of the epigenetic and expression data using limma (i.e., “regressed data” refers to the residuals of a model fitted with the covariates, excluding the predictor of interest, being AD diagnosis or conversion in this case). Spearman correlations were determined for the expression data and the average of the regressed beta values of the probes in the DMRs, DHRs, and DURs, as well as correlations between the different epigenetic markers (5mC, 5hmC, and UC) for these probes, using the Hmisc package (version 4.0-2).