Integrative epigenomic landscape of Alzheimer’s Disease brains reveals oligodendrocyte molecular perturbations associated with tau

Domain: chemistry
Published: 2026-03-03
Editorial outcome: major_revision

Abstract
Alzheimer’s disease (AD) brains have variable neuropathologic and biochemical changes. Capturing epigenetic factors associated with this variability can reveal novel biological insights into AD pathophysiology. Here, we conduct an epigenome-wide association study of DNA methylation in 472 AD brains with neuropathologic and biochemical brain protein levels core to AD pathogenesis. Using a novel regional methylation (rCpGm) approach, we identify 5478 significant associations, 99.7% of which associate with tau biochemical measures, and 93 concordant associations in external datasets. Transcriptome-methylome integration reveals enrichment in oligodendrocyte genes, including known AD risk gene BIN1, myelination genes MYRF, MBP and MAG previously implicated in AD, and novel genes like LDB3. Further characterization of these perturbations in independent AD and primary tauopathy datasets highlights consistent tau-related associations. In summary, we uncover the integrative epigenomic landscape of AD, demonstrate tau-related oligodendrocyte gene perturbations as a common potential pathomechanism across tauopathies and share findings via our Multiomic Atlas.

Introduction
Amyloid beta (Aβ) plaques and tau neurofibrillary tangles (NFT) are core neuropathological hallmarks of Alzheimer’s Disease (AD) 1,2. The characteristic distribution and levels of these amyloid and tau hallmarks in the brain are often measured by Thal phase 3 and Braak stage 4, respectively. These core measures show high inter- and intra-individual brain variability within AD 5–8 suggesting the presence of distinct mechanisms controlling the levels and biochemical states of their underlying proteins. This is supported by our previous findings that key proteins involved in the pathophysiology of AD including Aβ, tau, and apolipoprotein E (APOE) each have unique genetic architectures that associate with their specific biochemical states and levels in the AD brain 9. Characterization of these endophenotype specific molecular mechanisms in the AD brain are important as each protein and state has potentially unique functions and consequences 10. Insoluble states of Aβ42 and Aβ40 are major components enriched in Aβ plaques in the brain parenchyma and cerebral amyloid angiopathy (CAA) deposits, respectively 6,11. CAA is a co-pathology often seen in AD characterized by amyloid deposits in the cerebral vasculature which can impair blood flow and may contribute to clinical deficits beyond AD 12–14. Insoluble states of tau, particularly hyperphosphorylated tau (pTau), are core components of NFT which are often seen after brain levels of soluble total tau rise and become hyperphosphorylated 15. Recent studies have also shown that membrane bound states of Aβ and tau correlate with clinical severity and that tau increases its propensity to fibrillize when contacting a plasma membrane 16–21. Moreover, the APOE protein which is important in transporting and processing lipids in the brain, is co-deposited in insoluble Aβ plaques and its gene harbors the largest known genetic risk factor for AD, the APOE ε4 tagging variant 22–24.

Although there has been great progress identifying risk for AD including genetic and environmental/lifestyle factors 25,26, the underlying molecular mechanisms that these risk factors act through to impact specific AD endophenotypes remains to be fully characterized. It is likely that these risk factors affect both shared and distinct mechanisms that underpin the hallmarks of AD and investigating these factors will be important in discovering the molecular underpinnings of AD. Epigenetic mechanisms are known to mediate genetic as well as environmental and lifestyle factor effects on gene expression. Previous epigenetic studies of AD found broad dysregulation of epigenetic mechanisms corresponding to transcriptomic dysregulation 27,28. A key epigenetic regulator of gene expression is DNA methylation (DNAm). Previous studies have found differential DNAm sites near known AD risk genes including BIN1, APOE, APP, ADAM10, MAPT, SPI1, TREM2 and others 29–37. DNAm changes have also been linked to AD risk, onset and progression 30–32,38,39. These results suggest that epigenetic regulation, particularly that of DNAm, may be important in modulating expression levels of known AD risk genes and thereby influencing core AD endophenotypes. Nonetheless, large scale, integrative studies have yet to investigate how dysregulation of DNAm impacts specific brain endophenotypes beyond global neuropathology levels. Identification of AD-endophenotype specific molecular perturbations is required to deconstruct the underlying pathophysiology of AD.

We hypothesize that, within the AD brain, neuropathology, levels and biochemical states of core AD-related proteins are influenced by variations in DNAm. To test this, we performed epigenome-wide association studies (EWAS) using an innovate approach to analyze regional levels of DNAm (rCpGm) from the temporal cortex (TCX) and cerebellum (CER) with AD-related neuropathologic measures (Braak, Thal, and CAA) and brain biochemical levels of five proteins (APOE, Aβ40, Aβ42, tau, and pTau) (Fig. 1). Through our deep endophenotype EWAS and integrative omics approach leveraging a large cohort of neuropathologically confirmed AD cases, we identified widespread associations of DNAm with levels of distinct biochemical states of tau connected to significant dysregulation of nearby oligodendrocyte and myelin related gene expression. Interrogation of these tau related, epigenetically dysregulated oligodendrocyte genes and loci revealed concordant tau related associations across multiple independent datasets of AD and other tauopathies. Our findings suggest a key relationship between brain levels and biochemical states of tau and expression of oligodendrocyte genes regulated through DNAm perturbations. This may be a common pathomechanism across AD and other tauopathies. We broadly share our data and findings facilitated through our interactive web application, Multiomic Atlas of AD Brain Endophenotypes.

Results
We analyzed levels of DNA methylation at CpG sites using reduced representation bisulfite sequencing (RRBS) from two brain regions, TCX and CER, in a dataset of neuropathologically confirmed AD donors with multiple AD-related endophenotype measurements. These endophenotypes included neuropathological variables (Braak stage, Thal phase, and CAA scores) and TCX biochemical measures of five AD-related proteins 10 (APOE, Aβ40, Aβ42, total tau, and phospho-Tau (pTau) (Thr231)) from three tissue fractions (soluble (TBS), membrane-bound (TX), and insoluble (FA)) (Fig. S1). Our study dataset consisted of 472 individuals (53% female) with CpGm measured from the TCX and a subset of 200 from the CER (Supplementary Data 1). After stringent quality control and data reduction steps, a total of 455 TCX and 191 CER samples with 1,958,174 TCX and 1,858,631 CER CpGm sites were available for downstream analytics (Supplementary Data 2, Fig. S2).

Grouping CpGm by predicted chromatin state is a biologically relevant method to investigate regional levels of DNA methylation
A commonly used EWAS approach to identify AD-related CpGm loci is to associate individual CpGm sites with a phenotype of interest. However, this individual CpGm based EWAS does not account for methylation variations at neighboring CpGm sites which are highly correlated and likely work together to have a functional effect 40–43. Currently available methods address this by grouping nearby significant individual CpGm sites into regions, however, these approaches are based on post-hoc statistical groupings which lack biological basis 44–47. In this study, we aimed to investigate CpGm on a regional level defined by the biological function of the genomic region to identify functionally relevant CpGm loci. To do this, we developed a novel method that averages CpGm sites into regions (rCpGm) based on biologically based genomic windows defined by the Roadmap Epigenomics Consortium 15-chromatin state model generated from the human temporal lobe 48, a region pertinent to our human brain dataset (Fig. 1).

First, to determine whether our novel method yields biologically consistent results, we investigated the average methylation levels of rCpGm sites grouped by brain chromatin states 48. We found that the average methylation levels for each chromatin state were highly consistent with the predicted biology of these states (Fig. 2). Quiescent, repressed, and heterochromatin annotations had higher methylation, typically associated with lower gene expression levels 49, whereas active transcription start sites (TSS), flanking TSS, and bivalent/poised TSS regions had low methylation levels. Chromatin regions that are expected to have more variability across cell types, like enhancer regions, have around 50% methylation in our bulk brain tissue rCpGm data. Our chromatin state specific rCpGm averages were also highly consistent with published whole-genome bisulfite sequencing (WGBS) DNA methylation data 48 demonstrating the utility of our method for defining biologically relevant CpGm regions.

Next, we compared the average rCpGm levels between the TCX and CER brain regions for each chromatin state annotation. TCX and CER have mostly similar levels of rCpGm for each chromatin state with the largest differences in enhancer and transcription states which may reflect the distinct cell-specific epigenetic regulation of these two brain regions (Fig. 2). Investigating the 190 donors with both TCX and CER rCpGm measures, we found that average global rCpGm levels for each sample were similar between these brain regions (range: 0.27–0.34) and highly correlated in the same individual (r = 0.88–0.94, p < 2.2e-16), although TCX had consistently higher global methylation levels than CER (Fig. 2, Fig. S2).

We next explored the applicability of our novel method into two independent brain DNAm datasets measured via array technologies which are more often used in large-scale AD EWAS. These two datasets are the ROSMAP 30 (n = 686) and Brains for Dementia Research (BDR) 32 cohorts with CpGm measured from the dorsolateral prefrontal cortex (DLPFC) in both datasets and also the occipital cortex (OCC) in BDR (n = 581, n = 576) (Supplementary Data 3). Quality control (QC) processing of the DNAm levels from these datasets through our array-specific pipeline resulted in 429,723 CpGm sites in ROSMAP and 747,577 CpGm sites in BDR. Using our novel rCpGm method, we identified 136,201 rCpGms in ROSMAP and 203,435 rCpGms in BDR (Fig. S3). These independent datasets are similar in the levels and patterns of their average rCpGm for each chromatin state annotation across these datasets (Fig. S3) and in comparison to our dataset (Fig. 2), despite differences in donors, brain regions, collection sites and methylation measurement methods highlighting the replicability of our novel analytic approach. These average rCpGm patterns are also consistent even when examining only those rCpGms common to all datasets (Fig. S3). Notably, the array-based rCpGms have a narrower range of values (Fig. S3) in comparison to our RRBS-based measurements, particularly on the hypomethylated end which may potentially be due to lower precision and dynamic range for the former methodology 50,51. Even within array-based datasets, we observe differences across chromatin states where 450K arrays have collectively higher rCpGm averages than the EPIC v1 array. These differences indicate potential methodology-specific impacts on absolute methylation value measurements, which are unlikely to influence our results since we analyze methylation levels within each dataset. We do not recommend comparisons of absolute rCpGm values across datasets measured with independent technologies. Further, rCpGm levels for each chromatin annotation are more similar between DLPFC and OCC regions for the BDR cohort than for the DLPFC region across BDR and ROSMAP cohorts (Fig. S3). This suggests that inter-individual differences may be greater than intra-individual brain region differences in methylation.

Collectively, these results underscore the robust and replicable rCpGm measurements obtained with our analytic approach which demonstrate the broad utility and biological relevance of our novel method.

DNA methylation associates primarily with tau related brain endophenotypes
To determine DNA methylation (DNAm) associations with AD-related endophenotypes we took two EWAS approaches, namely the more common method of individual CpGm associations and regional methylation associations using our novel method that averages CpGms into regions based on the 15-chromatin state model 48. We also compared results from our novel regional CpGm association method (rCpGm) to the published post-hoc regional aggregation comb-p method 45.

Neuropathology endophenotypes (Braak stage, Thal phase, CAA score) were available for all study donors (i.e., both TCX and CER samples), whereas brain biochemical endophenotypes (protein levels of APOE, Aβ40, Aβ42, Total tau, and pTau measured in soluble (TBS), membrane-bound (TX), and insoluble (FA) brain tissue fractions) were measured only from TCX tissue (Fig. S1). We performed EWAS of individual CpGm levels from both TCX and CER with neuropathology and those from TCX also with biochemical endophenotypes adjusting for age at death, sex, and the first 3 genetic principal components. We identified 2489 FDR significant TCX CpGm associations with 10 endophenotypes and 50 FDR significant CER CpGm associations only for Thal (Fig. S4, Supplementary Data 4, Supplementary Data 5). The majority (97%) of the FDR significant TCX CpGm associations were with tau biochemical measures (n = 2415). The most frequently associated endophenotype was Total tau with 2268 FDR significant TCX CpGms, 21 of which were also associated with pTau. The remaining FDR significant CpGms were uniquely associated with a single endophenotype. Top TCX CpGm associations for each endophenotype were located in or near the genes ACSF3 (Total tau), SLC16A1-AS1 (Total tau), YES1 (AB40), TRAPPC14 (Thal phase), HR (pTau), RGS10 (APOE), PLEKHJ1 (AB42), CIB1/GDPGP1 (APOE), ARPC1B (AB42), and UMODL1 (APOE). Interestingly, CpGms near ACSF3, SLC16A1, HR, RGS10, PLEKHJ1, ARPC1B, and UMODL1 have previously been identified in meta-analyses with Braak Stage, although for different CpGm positions 32,52. ACSF3, ARPC1B, CIB1, HR, PLEKHJ1, RSG10, UMODL1 also had nominal CpGm associations with CERAD score and Thal phase 32. Our top 5 CER CpGm associations with Thal phase were near genes ANKRD33, LINC02314, RAB30, ASIC1, and UBE2Q2. A previous meta-analyses of cortical tissue found CpGms near RAB30 and UBE2Q2 to be nominally associated with Braak stage and Thal phase, respectively 32.

Using our novel regional methylation analysis method, we performed EWAS with 116,014 TCX and 116,811 CER rCpGms with the same endophenotypes. We identified 5478 FDR significant TCX rCpGm associations for 5 endophenotypes (Total tau, pTau, pTau, APOE, Aβ42) (Fig. 3, Supplementary Data 6). There were no FDR significant TCX or CER rCpGm associations with any of the neuropathology endophenotypes (Supplementary Data 6, Supplementary Data 7). Concordant with the individual CpGm findings, the majority of the significant TCX rCpGm associations were with Total tau (n = 4424) followed by pTau (n = 1029). There were 738 rCpGms significantly associated with brain levels of both Total tau and pTau. All 738 rCpGms had opposite directions of association with Total tau and pTau, which is consistent with the negative correlation of these two endophenotypes (ρ = −0.21, p < 0.001) 10. This suggests that Total tau and pTau may have shared but opposing molecular mechanisms that impact their levels in the brain. Each FDR significant rCpGm was on average composed of 23 individual CpGm sites (range: 1–895 CpGms) indicating that most of these associations are not driven by a single CpGm but are a composite methylation measure of the region. The majority (69%) of the FDR significant TCX rCpGm associations reside in more active chromatin regions including active TSS (TssA, N=578), Flanking TssA (N=652), transcription at gene 5’ and 3’ ends (N=47), Strong transcription (N=175), weak transcription (N=707), genic enhancer (N=256) and enhancer (N=1363) annotations (Fig. 3) suggesting that these significant rCpGms may be functionally important in gene expression regulation.

To ensure that these FDR significant rCpGm associations are robust to any biological and technical variations, we repeated these analyses after further adjustments for Braak stage and Thal phase, neuronal cell type proportion estimates, oligodendrocyte cell type proportion estimates, and sequencing batch. Comparisons between the primary and further adjusted models by Pearson correlation revealed high, significant positive correlations of both the coefficient estimates (r = 0.998-0.9997) and p values (r = 0.66–0.989) (Fig. S5). We also investigated if the normal distribution assumptions were violated or if these associations were driven by outliers via a non-parametric Spearman rank correlation test. We found 99.3% (5442/5478) of the FDR significant associations also had a nominally significant (p ≤ 0.05) association in the Spearman rank test (Fig. S5). These highly consistent test statistics between our primary and adjusted models and their persistent significance in a non-parametric test demonstrate the robustness of these associations to potential biological and technical confounders.

To compare our rCpGm approach to the previously published comb-p method 45, we applied comb-p which groups proximal CpGms into regions based a sliding window approach and association significance. We identified 2362 significant associations (Sidak p ≤ 0.05) in 1833 unique regions by comb-p (Fig. S6, Supplementary Data 8). Unlike our rCpGm approach, comb-p does not take into account direction of effect when determining significant regions. Further, since comb-p defines regions based on arbitrary windows and statistical significance, different regions are identified for each endophenotype EWAS. This hinders a direct comparison of comb-p regions across endophenotypes, a feature which is possible by our rCpGm approach. To resolve the comb-p shortcomings on effect direction and regional variability, we defined positive and negative comb-p associations by the effect direction for the majority of the individual CpGs making up that region. We considered regions across different endophenotype EWAS to be common if there was any genomic overlap between these regions. Comb-p based regional associations largely followed similar trends as the individual CpGm and rCpGm EWAS where Total tau had the highest number of significant associations (70.2%) followed by pTau (21.6%) and APOE (6.4%). Additionally, the largest overlap of significant associations between endophenotypes were for Total tau and pTau, all with opposing directions of association, consistent with rCpGm results. There were 1628 regional methylation associations common to comb-p and rCpGm, however rCpGm identified 3763 additional associations whereas comb-p uniquely identified 677 associations.

Although both the individual CpGm followed by comb-p regional grouping and our novel rCpGm analyses found the most associations with brain tau biochemical endophenotypes, there were twice as many FDR significant loci using the latter approach. Of the 2489 FDR significant individual TCX CpGms, 68% were included in at least one of the FDR significant rCpGms. Notably, 1671 individually significant CpGms were part of the 1000 rCpGms associated with Total tau, and 26 such CpGms were within 23 rCpGms associated with pTau. Taken together, our novel rCpGm approach identifies both those loci captured by the more traditional single CpGms followed by comb-p grouping but importantly also many others that are missed by the latter. This suggests that clustering CpGms based on chromatin states is both a more powerful EWAS approach and also provides biological context for these associations. Consequently, we focused our subsequent analyses on the rCpGm results and approach. Importantly, regardless of the CpGm analytic approach, most DNAm across the genome primarily associate with tau protein levels in the AD brain.

Replication of DNAm associations with tau neuropathology and other AD-related endophenotypes in independent datasets
We next sought replication of our findings in the ROSMAP 30 (n = 686) and Brains for Dementia Research (BDR) 32 (n = 581, n = 576) cohorts that have DLPFC DNAm measured with the Illumina 450K Human Methylation and both DLPFC and OCC DNAm measured with the EPICv1 850K arrays, respectively (Supplementary Data 3). Although to our knowledge there are no DNAm datasets that also have the precise brain biochemical measures of AD-related proteins available in our dataset, ROSMAP and BDR are two of the largest cohorts with both brain DNAm and AD-related endophenotypes.

We already demonstrated the broad similarities of rCpGm levels and patterns for each chromatin state annotation across these datasets and ours (Fig. S3, Fig. 2), despite the fact that the individual CpGm sites measured by their DNAm platforms are not identical. To determine whether overlapping rCpGms that are obtained from different datasets correlate with each other, we performed pairwise correlations (Fig. S7). Using DNAm levels averaged across all samples within each dataset for all rCpGms overlapping with those from the other datasets, we conducted Spearman correlations. There are strong positive correlations for DNAm levels of overlapping rCpGm sites for each pairwise dataset analysis. As expected, these correlations are strongest for BDR DLPFC and BDR OCC datasets (Spearman Coefficient = 0.996) obtained from the same brain donors, followed by BDR DLPFC and ROSMAP DLPFC that have measurements from the same brain region both using similar array technologies (Spearman Coefficient = 0.914) and BDR OCC and ROSMAP DLFPC (Spearman Coefficient = 0.918). Nonetheless, overlapping rCpGms from our RRBS-based TCX also had high correlations with those from the other array-based datasets (Spearman Coefficients: 0.74–0.828) (Fig. S7). Thus, while not identical to our DNAm measurements, rCpGms from these other AD brain cohorts could serve as replication datasets for our findings.

For replications, we focused on the 738 TCX rCpGms that are significantly associated with brain levels of both Total tau and pTau. These 738 rCpGms are high confidence both due to their FDR-significant and biologically congruent directions of associations with two brain biochemical AD endophenotypes. Of these 738 rCpGms, 600 were present in at least one of the replication datasets. Although BDR and ROSMAP did not have brain biochemical AD-related protein measurements, they both have Braak stages, which is a measure of tau pathology. To consider an rCpGm as having a biologically congruent association in the replication datasets, the direction of effect would be the same for Braak stage and pTau and opposite for Total tau. This is based on prior work 9 and positive correlations of Braak stage with pTau (Spearman correlation ρ = 0.289, pval = 3.72E-10), but negative with Total tau (Spearman correlation ρ = −0.30, pval = 6.29E-11) in our dataset (Fig. S8).

Using these criteria, 93 out of 600 tested rCpGms had nominally significant (p ≤ 0.05) and biologically congruent associations with Braak stage in BDR and/or ROSMAP (Fig. 4, Supplementary Data 9, Supplementary Data 10). Of these 93, 3 replicated in all 3 datasets (tier 1), 23 in 2 of the datasets (tier 2), and 67 in 1 of the replication datasets (tier 3). We also performed a meta-analysis of Braak stage including the RRBS-TCX, BDR-DLPFC, and ROSMAP-DLPFC datasets and found that 47 of the 93 replicated rCpGms were significant (fixed p ≤ 0.05) including all of tier 1 rCpGms, 16/23 tier 2 rCpGms, and 28/67 tier 3 rCpGms (Fig. 4, Supplementary Data 11). We further tested the association of these rCpGms with Thal phase and AD diagnosis in the replication datasets and also found biological congruence of associations across these additional AD measures (Fig. 4, Supplementary Data 9, Supplementary Data 10). In other words, rCpGms with FDR-significant negative associations with Total tau and positive associations with pTau in our dataset, had positive associations with Braak, Thal and AD diagnosis in one or more of the replication datasets and vice versa.

When we annotated these 93 rCpGms for proximal (±1 Mb) variants identified in AD genome-wide studies 53, we found that a number of them reside near known AD risk loci including JAZF1, SIGLEC11, SCIMP, KLF16, BIN1, SPI1, ABI3, ANKH, WDR12, and SIGLEC11 53 (Supplementary Data 12). Collectively, these results support independent replication of our 93 high confidence TCX rCpGms for their association with brain tau-related and other AD endophenotypes, which may also have implications for AD risk.

Tau related rCpGms associate with expression of oligodendrocyte and myelin related genes
Most of the FDR significant AD endophenotype associated rCpGms in our dataset are located in active chromatin regions (Figs. 3, 4). To further characterize the associations between these rCpGms and brain gene expression levels, we leveraged bulk CER and TCX RNAseq data from the same individuals. After RNAseq QC, we paired rCpGms to genes ±500kb from the start and end of each rCpGm yielding 1,612,499 TCX and 1,834,887 CER cis-rCpGm-gene pairs. We then tested the association between cis-rCpGm-gene expression pairs using linear regression.

We identified 12,973 TCX and 1413 CER FDR significant cis-rCpGm-gene pairs (FDR p ≤ 0.05) (Supplementary Data 13, Supplementary Data 14). Many of the significant cis-rCpGm-gene pair associations had rCpGms annotated to expression regulation terms including active TSS, flanking active TSS, and enhancer sites (51% in TCX and 32% in CER). Genes previously reported to have AD-related methylation changes nearby, including RHBDF2 30,31, IRF8 31, CCND1 31, and SPI1 32 had both FDR significant rCpGm associations with AD endophenotypes and cis-rCpGm-gene expression associations in our dataset (Supplementary Data 13, Supplementary Data 14).

To further annotate the 93 high-confidence tau and other AD-endophenotype associated rCpGms, we tested their corresponding 1464 cis-rCpGm-gene expression pairs (1225 unique genes). Many of these 93 rCpGms reside in more active chromatin regions such as Active TSS, Flasking TSS, or Enhancers, suggesting their functional consequences in gene expression (Fig. 4). Indeed, when tested for their effects on brain gene expression, these rCpGms had FDR significant associations with nearby genes measured from the same donors (Supplementary Data 12–14). There were 535 significant cis-rCpGm-gene expression associations (459 unique genes) which included previously implicated AD genes like BIN1 23,53,54, ANXA2 55,56, CDK18 55,57,58, AKAP13 59, RORA 60,61, KIFC1 62, and KLK6 63,64 (Supplementary Data 12).

To determine if this set of 459 genes, which correspond to the 93 tau-related high confidence rCpGms, implicates common pathways or cell types, we performed gene set enrichment analysis. The top three nominally enriched gene ontology terms for these genes relate to positive regulation of myelination (GO:0031643, p=3.4E-5), gliogenesis (GO:0042063, p=1.4E-4), and microtubule nucleation (GO:0007020, p=1.5E-4) (Supplementary Data 15). Consistent with these GO terms, this gene set was also enriched for oligodendrocyte marker genes. Using the top 1000 BRETIGEA 65 cell type markers genes, there were 45 oligodendrocyte genes amongst these 459 (p=9.8E-05), 11 of which were also in the top 100 BRETIGEA 65 markers (p=4.0E-05). The set top 11 oligodendrocyte marker genes with rCpGm associations are MAG, KLK6, ST18, MYRF, LDB3, MBP, CD22, TGFA, RTKN, SOX8, and CLMN. We next investigated these 11 genes and an additional gene CKD18 further in both our and external datasets. Although part of the 1000, and not the 100 top oligodendrocyte marker genes 65, CKD18 was previously implicated in tau phosphorylation including at the T231 site and tau aggregation 58,66, found contained within NFT 58, differentially expressed in TCX tissue of individuals with AD compared to controls 58, associated with AD DNAm changes 67, and helps prevent DNA damage and promotes genomic stability 68, thereby leading us to include this gene in additional analyses.

All 12 oligodendrocyte genes have negative associations of their brain expression and paired rCpGm methylation levels (Fig. 5). There were 16 rCpGm-gene expression pairs for these 12 oligodendrocyte genes and 15 unique rCpGms. All 15 rCpGms had FDR significant negative associations with Total tau and positive associations with pTau (Fig. 5). These rCpGm~tau endophenotype associations remained significant even after adjusting for oligodendrocyte cell proportion estimates (Supplementary Data 16), indicating that they are not driven by changes in proportions of oligodendrocytes in bulk brain tissue. These rCpGms had additional nominally significant associations with other AD-related endophenotypes, the most consistent of which were negative associations with APOE levels and positive associations with pTau (Fig. S9). Almost all of the individual CpGms that comprise these oligodendrocyte rCpGms have the same direction of association with Total tau or pTau as their corresponding rCpGm, although few individually reach FDR significance (Fig. S10). This again underscores the superior power of our novel and biologically relevant rCpGm in comparison to individual CpGm analysis in detecting DNAm associations.

We further evaluated the top oligodendrocyte genes with tau-related rCpGm associations using our dataset and independent studies. Within the same brain donors, TCX expression levels of these oligodendrocyte genes also had inverse directions of associations with TCX levels of Total tau and pTau. Of these 12 top oligodendrocyte genes, 5 nominally associated with both Total tau and pTau levels (CDK18, LDB3, MBP, CD22, CLMN) and 2 with pTau (KLK6, ST18) in the TCX (Fig. 5A, Supplementary Data 17). All but one of these associations (TGFA-pTau) had directions of effect that are biologically congruent with both rCpGm-tau and rCpGm-gene associations. Notably, TGFA, unlike the other marker genes, is not only expressed in oligodendrocyte cells but is also expressed in oligodendrocyte precursor cells which may underlie this discrepancy. Upon adjustment for oligodendrocyte cell proportion estimates, CDK18, MBP, and CD22 remained nominally associated with both Total tau and pTau levels with consistent effect directions (Supplementary Data 17).

Further investigation of the 15 rCpGms in Illumina 450k DNAm data derived from NeuN+ and NeuN- sorted occipital cortex nuclei 36, revealed 14 rCpGms with positive associations with AD in the glia enriched NeuN- population while only 6 rCpGms had positive associations with AD in the neuronal enriched NeuN+ population (Fig. S11). The regression estimate magnitudes of the glia enriched NeuN- population were also larger than the NeuN+ population suggesting that the associations identified in the bulk DNAm data are likely originating from DNAm changes in the glial cells rather than neuronal cells 32.

Using three independent brain bulk RNAseq datasets of AD and control brains, we next investigated the association of these top oligodendrocyte genes with Braak stage in the Mayo TCX 69, Mount Sinai BM22 70, and ROSMAP 71 datasets from the AMP-AD study 72. Expression of LDB3 was associated with Braak stage in the Mayo TCX (beta = −0.19; FDR p = 2.54E-2) and ROSMAP (beta = −0.058; FDR p = 2.52E-2) datasets with a biologically congruent negative direction of effect (Supplementary Data 18). Using AMP-AD data, MYRF and CDK18 have also been nominated as potential therapeutic targets for AD 73, providing additional support for a role of these oligodendrocyte genes in tau-related outcomes.

To determine whether these brain rCpGm-oligodendrocyte gene-tau associations are unique to AD vs. conserved across other tauopathies, we investigated available data from methylation or gene expression studies of progressive supranuclear palsy (PSP) and Pick’s disease. Further, we analyzed oligodendrocyte cell clusters from published brain single nucleus datasets of tauopathies to establish these associations in the cell type of interest. There was one brain DNAm study conducted from the prefrontal cortex of PSP and control donors 74 to which we applied our QC, rCpGm method and conducted an EWAS of PSP risk and Braak stage (90 PSP and 66 controls) (Supplementary Data 3 and 19, Fig. 5B). We determined biologically congruent associations with Braak stage across the 12 top oligodendrocyte rCpGms that also had consistent PSP risk associations for all but three of the genes (TGFA, ST18, CDK18). This finding supports a model where rCpGm-oligodendrocyte gene-tau associations are not unique to AD but may be a shared disease mechanism across tauopathies.

We next analyzed all available single nucleus RNAseq (snRNAseq) data from sizeable brain datasets of AD and other tauopathies 75–81 to replicate the top oligodendrocyte gene expression associations with disease risk and tau endophenotypes in an oligodendrocyte cell-type specific manner. These independent snRNAseq datasets and brain regions were from Is et al. 75 (24 TCX, AD and controls), Mathys et al. (427 prefrontal cortex, AD and controls) 76, Gamache et al. (24 TCX, AD and controls) 77, Morabito et al. (18 prefrontal cortex, AD and controls) 78, Anderson et al. (15 dorsolateral prefrontal cortex, AD and controls) 79, Rexach et al. (41 insular cortex, AD, PSP, Pick’s disease and controls) 80, and Min et al. (34 TCX, PSP and controls) 81. In addition to these snRNAseq data from a total of 583 independent brain samples from multiple regions and diagnoses, Min et al. 81 also had bulk RNAseq data from 408 PSP and control TCX samples. We tested the top brain rCpGm associated oligodendrocyte genes for their oligodendrocyte-specific expression associations with risk for AD 75–80, PSP 80,81, Pick’s disease 80, tau-related (Fig. 5B) and other available endophenotypes (Fig. S12) in these datasets. Despite differences in the cohorts, brain regions and diagnoses, there was highly consistent association for oligodendrocyte-specific expression for many of these top 12 genes with these outcomes in a direction that would be expected based on the rCpGm-oligodendrocyte gene-tau associations we identified. Most of these genes had lower expression levels in AD and other tauopathies and were negatively associated with tau-endophenotypes such as Braak stage and NFT. Most consistent associations were observed for LDB3, MYRF, MBP, MAG, TGFA, and CD22 loci where the rCpGms positively associate with pTau, negatively with Total tau and brain gene expression in our data (Fig. 5A, Supplementary Data 16–18). These rCpGms associate positively with PSP risk and Braak in Weber et al. 74 DNAm data (Fig. 5B, Supplementary Data 19).

To determine whether genetic variants could account for the tau-related DNAm changes observed for these top oligodendrocyte gene related rCpGms, we performed a cis-methylation quantitative trait locus (mQTL) analysis between these rCpGms and imputed genetic variant dosages within a distance of ±500kb. There were 160 FDR significant SNP-rCpGm associations (Supplementary Data 20). The rCpGms implicated by the mQTL associated with ST18 (14 SNP-rCpGm pairs), MBP (122), CLMN (2), MAG & CD22 (2), KLK6 (1), and RTKN (19) gene expression. Of these 160 unique SNPs, 13 near MBP nominally (p ≤ 0.05) associated with Total tau levels, the top 2 of which were also trending towards association (p < 0.052) with pTau levels (rs13381407, rs12457073: R2=1, D’=1) 9 (Supplementary Data 20), but none associated with AD risk 53. These results suggest that genetic variation is unlikely to account for the rCpGm-gene-tau relationships at these loci, with the possible exception of the MBP locus.

Collectively, these findings support a model for tauopathies such that with increasing brain levels of tau neuropathology, pTau, and decreasing Total tau, DNA hyper-methylation occurs in oligodendrocyte/myelin gene regions leading to reduced expression of these genes and potentially downstream oligodendrocyte dysfunction and dysmyelination (Fig. 5C). However, it is equally plausible that other untested genetic or environmental factors lead to DNA hyper-methylation, subsequently leading to downregulated myelin genes and oligodendrocyte dysfunction, culminating in brain accumulation of toxic tau forms 82.

We have compiled our findings and data in our interactive web application, Multiomic Atlas of AD Brain Endophenotypes. Our web application has the results from the following analyses: GWAS: Association between genetic variant dosages and TCX brain biochemical measures. rCpGm EWAS: Association between regional CpGm (rCpGm) levels and neuropathology or TCX brain biochemical measures. Individual CpGm EWAS: Association between individual CpGm levels and neuropathology or TCX brain biochemical measures. mQTL: Association between regional CpGm (rCpGm) levels and nearby genetic variant dosages (+/- 500kb). rCpGm to GE: Association between regional CpGm (rCpGm) levels and nearby gene expression measures, where rCpGms were tested for association with genes +/- 500kb from start and end of region. These data can be searched by gene names or chromosomal locations, and results are presented as both tables and interactive plots.

Methods
This study was approved by the appropriate Mayo Clinic Institutional Review Board.

Brain samples
Post-mortem superior temporal gyrus cortex (TCX) and cerebellum (CER) tissue samples were collected by the Department of Neuroscience, Mayo Clinic Brain Bank (Jacksonville, Florida, USA). These samples are a part of the Mayo Clinic AD-CAA (MC-CAA) study, previously described 9,84 and also detailed on the AD-knowledge portal (https://adknowledgeportal.synapse.org, see Data Sharing section). All samples had a confirmed AD neuropathological diagnosis with a Braak stage ≥ 4 and Thal phase ≥ 2, were scored for CAA pathology, had an age at death greater than 55 years, and due to availability, were all non-Hispanic white (Supplementary Data 1). Donors with co-pathologies were not excluded to allow for inclusion of the full spectrum of neuropathologies in AD. Braak stage, Thal phase, and CAA scores were measured using previously established protocols 3,4,10,146,147. Intermediate Braak stages were grouped with the next lowest stage as follows, stage 3.5 is 3, 4.5 is 4, and 5.5 is 5 as detailed previously 10,84 (Supplementary Data 1, Fig. S1).

Biochemical measures
Biochemical measures from a subset of superior temporal cortex brain samples included the MC-CAA study described previously 9,10 were utilized for this study. Biochemical measures include five AD-related proteins (APOE, Aβ40, Aβ42, tau, and phospho-tau (Thr231)) from three tissue fractions (soluble, membrane-bound, insoluble). Briefly, supernatant fractions were collected after three sequential buffer treatments of tissue homogenate and resulting pellets: first with tris-buffered saline buffer (TBS), second with detergent (TBS/1% Triton X) buffer (TX), and lastly with formic acid buffer (FA). Protein quantification was performed via ELISA for each fraction and normalized against total protein quantities. All biochemical measures were transformed by either the natural log or square root to approximate a normal distribution including the Aβ40/42 ratio. Transformed measures were then standardized with a mean of zero and standard deviation of one using the z-score formula to allow for comparisons of effect sizes between biochemical measure analyses. Of the 469 MC-CAA samples with biochemical measures, 452 also had DNA methylation data (Fig. S1).

Genotype data
Genome-wide genotypes (GWG) were generated in previous studies and are described in detail elsewhere 9,84. Briefly, DNA was isolated for GWG from brain tissue using the AutoGen 245T instrument according to manufacturer’s instructions including an RNase A (Qiagen) digestion. GWG were collected using the Infinum Omni2.5 Exome 8 v1.3 genotyping array, formatted in PLINK (v1.9) 148,149 and underwent rigorous quality control (QC). Variants passing QC were imputed to the haplotype reference consortium (HRC) keeping only variants with an imputation quality R2 ≥ 0.7 and minor allele frequency (MAF) ≥ 2%. Variant base pair (BP) positions were lifted from hg19 to hg38 using the liftOver package and corresponding chain files provided by UCSC Genome Browser. After mapping coordinates in hg19 to hg38, positional information in the plink files were updated using the ‘--update-map’ flag. Of the N=6,724,981 imputed and quality controlled genetic variants aligned to the hg19 reference genome, N=6,721,677 successfully lifted over to the hg38 reference and could be analyzed.

DNA methylation
DNA from 40-70 mg of TCX tissue was isolated using the AutoGen 245T instrument according to manufacturer’s instructions including an RNase A (Qiagen) digestion prior to loading. The CER samples failed DNA isolation using the AutoGen instrument, so DNA was isolated from 10 to 15 mg of CER tissue using the same AutoGen Reagent kit and instrument steps manually. Pilot experiments were performed comparing TCX DNA isolated in triplicate from 3 samples using the AutoGen instrument vs manual methods as well as 2 CER samples in triplicate manually to determine if extraction method significantly impacted overall DNA methylation (DNAm) levels (Fig. S13). Following DNA extraction of all samples, Reduced Representation Bisulfite Sequencing (RRBS) was performed on 25 ng of genomic DNA for each sample by the Mayo Clinic Genome Analysis Core (GAC) in Rochester, MN following manufacturer’s recommendations for the NuGen Ovation RRBS Methyl-Seq System 1–16 (Tecan Genomics). The EZ DNA Methylation Kit (Zymo Research) was used for bisulfite conversion. RRBS libraries were sequenced on the HiSeq 4000 (Illumina) with 8 samples per lane and 50bp paired end reads. The Streamlined Analysis and Annotation Pipeline for Reduced Representational Bisulfite Sequencing (SAAP-RRBS) 150 was used to align reads to the hg38 reference genome and to get CpGm ratio information. Briefly, FASTQs were trimmed to remove diversity and adaptor sequences using cutadapt, and any reads with less than 15bp were discarded. Trimmed FASTQs were then aligned against the hg38 reference genome using BSMAP 151. Samtools 152 was used to get mpileup and custom PERL scripts were used to determine CpG and non-CpG methylation levels along with bisulfite conversion ratios. All samples had a conversion ratio > 98% and we kept only CpG sites with at least 10x coverage. One triplicate TCX sample (TCX_Autogen_1c) had low read coverage and so was not included in the downstream comparisons (Fig. S13a). In a Principal Component Analysis (PCA) of the top 50% most variable CpG sites with 10x coverage (N=535,165 CpGs), we found that sample replicates and brain tissue types were tightly clustered independent of extraction method (Fig. S13b). We also found CpGm levels were highly correlated within sample replicates regardless of extraction method and to a much greater degree than between tissue types of the same samples (Fig. S13c). These results indicated that extraction method (instrument vs manual) does not significantly contribute to differences in DNA methylation levels. Consequently, we performed automated DNA extraction for TCX samples and manual for CER samples.

All study samples were randomized for brain tissue (TCX, CER), sex, age at death, APOE ε4 dose and CAA score for sample prep and sequencing. Library preparation and sequencing was performed in the same manner as that for the pilot study. One TCX sample failed sequencing resulting in N=200 CER and N=471 TCX samples with available DNAm data. Raw sequencing reads were analyzed using the SAAP-RRBS pipeline described above. Rigorous QC of the CpG sites was then performed independently for each brain tissue type by (1) removing samples with a bisulfite conversion ratio < 98%, (2) removing sites with less than 10x total read coverage and those that mapped to the mitochondrial chromosome, (3) removing samples that failed genome-wide association study QC performed in a previous study 84, (4) removing sites with a total read coverage in the 99.9th percentile likely indicative of PCR bias, (5) removing samples that failed sex check by calculating the average X and Y chromosome read coverage normalized to the average read coverage of chromosome 1 defining failed samples as those that are beyond 4 standard deviations of the mean of their respective female and male clusters, (6) removing CpG sites that mapped to the sex chromosomes, (7) masking individual sites for each sample separately that are C to T nucleotide transitions on the forward strand or G–A transitions on the reverse strand as determined from the imputed genotype data 84, (8) removing outlier samples determined from a principle component analysis (PCA) of the top 50% most variable CpGm sites based on standard deviation using the R function prcomp(.scale=TRUE, center=TRUE) with outliers defined as those 4 times the standard deviation from the mean of principal components (PC) 1-3, (9) removing sites not present in over 50% of the samples, and (10) removing samples with less than 50% of the CpGs that passed QC (Supplementary Data 2).

For analyses investigating regional CpGm (rCpGm), we leveraged the ChromHMM core 15-state chromatin model developed by the Roadmap Epigenomics Consortium 48 for the adult brain Inferior Temporal Lobe (E072). The hg38 lift predicted chromatin bed file was downloaded from https://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/coreMarks/jointModel/final/ on June 29, 2022. Although there were regions in this file that were not annotated to a specific chromatin state as a result of the lift over from hg19 to hg38, these individual regions were kept and annotated as “16_NA” or “NotDefined” for downstream analyses. We grouped all CpGm sites that passed QC into regions based on the 15-state chromatin model annotations and averaged each group to determine rCpGm levels (Fig. 1). rCpGm start and stop sites were defined as the first and last CpGm present in each rCpGm group. As a predicted chromatin state model for the CER did not exist at the time of this study and the fact that we want to analyze rCpGms across tissue types, we used the Inferior Temporal Lobe model to define rCpGms for all brain tissue types. In this study, we had two tiers of DNA methylation measures, one at an individual CpGm site level (CpGm) and one at a regional level according to our novel method defined by the chromatin state model (rCpGm). For each tier, we then performed data reduction steps to remove non-variable sites/regions that are not statistically informative and that would unnecessarily contribute to increasing the multiple testing correction burden. These data reduction steps include, (11) removing CpGm/rCpGm with a methylation ratio range less than 1% as these are likely not biologically informative, and (12) removing sites/regions where 95% or more of the samples are fully methylated (100%) or fully unmethylated (0%) as these provide unreliable statistical estimates (Supplementary Data 2).

Gene expression measures
As part of the MC-CAA study, bulk RNAseq data was also generated for 477 TCX and a subset of 200 CER samples. Sample Preparation and Sequencing: Total RNA was extracted for each tissue sample using Trizol® reagent and cleaned using Qiagen RNeasy columns with DNase treatment. RNA integrity (RIN) was measured using an Agilent Technologies 2100 Bioanalyzer and only samples with a RIN > 5.5 were sent to the Mayo Clinic Genome Analysis Core (GAC) for library preparation and sequencing. The TruSeq RNA Library Prep Kit v2 (Illumina, San Diego, CA) was used for library preparation (non-stranded). Library concentration and size distribution were determined on an Agilent Bioanalyzer DNA 1000 chip. Samples were randomized across sequencing flow cells with respect to average CAA score, sex, age at death, and APOE ε4 genotype. The Illumina HiSeq 4000 was used for sequencing with 100 bp paired-end reads, multiplexing six samples per flow cell lane. Quality Control: The Mayo Clinic MAP-Rseq pipeline v.3.0.2 153 was used for read alignment to hg38 implementing STAR v2.5.2b 154, and featureCounts from the Subread package v1.5.1 for counting 155. Raw read counts were log2-transformed and normalized using conditional quantile normalization (CQN) via the Bioconductor package; accounting for sequencing depth, gene length, and GC content 156. FastQC was used for QC of raw sequence reads, and RSeQC 157 was used for QC of mapped reads. We excluded samples with < 50 million reads mapped to genes, < 50% of reads mapped to genes, and outliers (1.5 inter-quartile range) based on distribution of % junction reads. Sex check was performed by comparing expected expression of Y chromosome genes to recorded sex. Samples that were outliers (> 4 standard deviations from the mean) on PCs 1 and 2 in a PCA (prcomp, CPM, R statistical software) and multiple other QC variables (RSeQC) were also excluded. We also excluded samples that failed previous genetic based QC steps including relatedness and heterozygosity checks as described previously 84. There were 456 TCX and 187 CER samples that passed QC and of these, 449 TCX and 184 CER samples had matching RRBS data. We included in our downstream analysis only autosomal genes expressed above background which we defined based on the distribution of CQN values filtering out genes with a CQN value less than 1 for both the CER and TCX datasets.

AMP-AD datasets
The Mayo RNAseq study 69, The Mount Sinai Brain Bank (MSBB) study 70, and The Religious Orders Study and Memory and Aging Project (ROSMAP) Study 71 were obtained from the AD-knowledge portal (https://adknowledgeportal.synapse.org). The RNA-seq data previously underwent consensus reprocessing by Accelerating Medicines Partnership-AD (AMP-AD, RNAseq Harmonization Study) 72. Additional QC and diagnosis harmonization of these datasets based on neuropathological measures are described in detail elsewhere 84.

Statistical analyses
Power calculation. Power calculations were performed with the pwr.r.test function from the pwr package v1.3-0 in R (v4.0.2). For a sample size of 472 in a two-sided test, we have 80% power to detect a correlation coefficient of 0.239 and 0.129 at an alpha of 1E-05 and 0.05, respectively. For a sample size of 200 in a two-sided test, we have 80% power to detect a correlation coefficient of 0.359 and 0.197 at an alpha of 1E-05 and 0.05, respectively.

Phenotypes. Spearman rank correlation between all biochemical measures and CAA have been reported previously 10.

Cell type proportion estimates. Neuronal and oligodendrocyte cell type proportions were estimated from the bulk tissue RNAseq data using the Digital Sorting Algorithm (DSA) 158 via the DSA R package (v1.0). The top 100 genes from BRETIGEA 65 were used as the cell type specific marker genes for these estimations. Oligodendrocyte precursor cell (OPC) types were not investigated as we previously did not observe robust correlations among the defined OPC marker genes as described in Wang et al. 93.

Epigenome wide association studies (EWAS). Primary association models were run using the linear regression lm() function in R v4.1.2 adjusting for age at death, sex, and the first three genetic PCs. The CpGm variable (CpGm or rCpGm) was included in each association model as the dependent variable. There were three samples with a Thal phase measure of 2 (2 TCX and 1 with both TCX and CER) that passed DNAm QC. Given the low number of this measure in a single group and the possibility that small variations in this group may considerably and unreliably impact our statistical estimates, we set these Thal phase 2 measures to NA in analyses that included Thal phase as a variable. False Discovery Rate (FDR) corrected p values were calculated with the R p.adjust function using the Benjamini & Hochberg method.

Comb-p regional aggregation. To compare our rCpGm method with the published regional aggregation method comb-p 45, test statistics from each individual CpGm EWAS were independently used an inputs to the python based comb-p software following the default “pipeline” steps 45. Comb-p identifies and groups spatially correlated p values through a sliding window approach. Input parameters included regions with ≥ 3 CpGs, a seed of < 1E-04, a sliding window size of 500bp, and significant regions defined as those with a Sidak corrected p ≤ 0.05. While comb-p does not take into account the CpG regression estimates, significant regions were defined as positively or negatively associated with a phenotype by the effect directions of the majority (> 50%) of CpGs included in the defined region. When comparing across EWAS, overlapping regions were identified and concatenated using the GenomicRanges::reduce function.

Association with gene expression measures. In the MC-CAA dataset, there were N=449 TCX and N=184 CER samples that had both RRBS and bulk RNA seq data that could be investigated. Cis-rCpGm-Gene pairs were identified by connecting rCpGm with all genes ±500kb away from the start and end of each rCpGm for the TCX and CER separately. Residuals of rCpGm were extracted by fitting the linear regression model: rCpGm~AgeAtDeath+Sex+genetic PC 1–3. Bulk gene expression (GE) residuals were extracted by fitting the mixed effects linear regression model: GE~AgeAtDeath+Sex+RIN+(1|batch). The association between rCpGm and bulk gene expression residuals were tested by fitting the linear regression model: Gene_Expression_residuals~rCpGm_residuals for each cis-rCpGm-Gene pair independently for each tissue type.

To identify gene expression associations with the AD endophenotypes, we tested the association of normalized gene expression data with the endophenotype of interest. For continuous phenotypes (biochemical measures and CAA), we ran the following linear regression model: cqn(gene expression)~ADphenoytpe+Age at Death+sex+RIN+(1|flowcell) using the lmer function in the lme4 package 159. For ordinal variables (Braak stage and Thal phase), we first extracted residuals from the gene expression values using the following model: cqn(gene expression)~Age at Death+sex+RIN+(1|flowcell), and then associated the gene expression residuals with each phenotype in an ordinal regression model using the clm function in the ordinal R package. As with the DNAm data, because there were only a small number of samples with a Thal phase of 2, these were set to NA in analyses that included Thal as a variable. FDR p values were calculated with the R p.adjust package using the Benjamini & Hochberg method. The AMP-AD datasets were used to investigate the association between CQN normalized gene expression levels of implicated genes with Braak stage in non-Hispanic white, neuropathologically defined AD cases and controls using linear regression in R with the lm() function. Mayo TCX models adjusted for age at death, sex, RNA integrity number (RIN), and sequencing flow cell rounding down intermediate Braak stage numbers to the next lowest stage. MSBB BM22 models adjusted for age at death, sex, RIN, and batch. ROSMAP models adjusted for age at death, sex, RIN, sequencing batch, and study (ROS/MAP). Individuals with an age at death above 90 were set to 90, in compliance with HIPAA rules. Oligodendrocyte proportions were further adjusted for in each model by including the OLIG2 (ENSG00000205927) normalized expression levels as a covariate.

mQTL analysis. Cis-methylation Quantitative Trait Locus (cis-mQTL) analysis was performed using the R package MatrixEQTL 160 for all genetic variant dosage and rCpGm pairs ±500kb apart for the TCX and CER separately adjusting for age at death, sex, and the first 3 genetic PCs. In total, we investigated 289,078,059 TCX rCpGm-SNP pairs and 291,400,204 CER rCpGm-SNP pairs. False Discovery Rate (FDR) corrected p values were calculated with the R p.adjust package using the Benjamini & Hochberg method.

Sensitivity analyses. We performed post-hoc sensitivity analyses to test the robustness of the top associations to effects of batch, neuropathological changes, neuronal and oligodendrocyte cell type proportions. To determine if top associations were influenced by batch, we included batch in a linear regression model as a random effect variable using the lmer() function of the lme4 R package 159: rCpGm~Endophenotype+Age at death+Sex+Genetic PCs 1-3+(1|FlowCell). To determine if top biochemical associations were influenced by neuropathological changes, we included Braak stage and Thal phase as fixed-effect covariates in a linear regression model: rCpGm~Endophenotype+Age at death+Sex+Genetic PCs 1-3+Braak stage+Thal phase. To determine if top associations were influenced by neuronal cell type proportions, we included DSA estimated neuronal cell type proportions as a fixed-effect covariate in the linear regression model: rCpGm~Endophenotype+Age at death+Sex+Genetic PCs 1-3+Neuronal Cell Proportion. To determine if top associations were influenced by oligodendrocyte cell type proportions, we included DSA estimated oligodendrocyte cell type proportions as a fixed-effect covariate in the linear regression model: rCpGm~Endophenotype+Age at death+Sex+Genetic PCs 1-3+Oligodendrocyte Cell Proportion. Lastly, to determine if distribution of a phenotype influenced the top associations, we performed Spearman rank correlation analysis in R using the cor.test() function of the stats package.

Pathway analysis. Pathway enrichment was tested using the Gene Ontology Biological Processes database in the anRichment R package with p values computed via a one-sided hypergeometric test.

Annotations. CpGm sites were annotated according to their regions (Island, shore, shelf, inter-genic) and nearby genes (±50kb) using the R package annotatr() v1.21.0 161 querying AnnotationHub_3.2.2 in R (v4.1.2) or with built-in hg38 genes, respectively.

DNAm replication datasets. Data from The Religious Orders Study and Memory and Aging Project (ROSMAP) Study was downloaded from the AD knowledge portal for 707 Non-Hispanic White, dorsolateral prefrontal cortex (DLPFC) samples with DNAm measured on the Illumina HumanMethylation450 BeadChip array 30 including IDAT and metadata files on May 17th, 2023 (Supplementary Data 3). Additional and updated metadata for these ROSMAP samples was also obtained from the Rush Alzheimer’s Disease Center (RADC) on May 22nd, 2023. Data from the Brains for Dementia Research (BDR) Dataset was downloaded for 630 individuals with DLPFC and occipital cortex (OCC) DNAm measured on the Illumina EPICv1 DNA methylation array 32 including IDAT files from the Gene Expression Omnibus (GEO, accession number GSE197305) on February 23rd, 2023 and metadata from the UK Brain Bank Network (UKBBN) database after project approval (https://brainbanknetwork.ac.uk/) on May 26th, 2023 (Supplementary Data 3). PSP data from Weber et al. 74 was downloaded from GEO (GSE75704) on April 12th, 2023 for 94 PSP and 72 control individuals with prefrontal lobe DNAm measured using the Illumina HumanMethylation450 BeadChip. Sorted nuclei data from Gasparoni et al. 36 was downloaded from GEO (GSE66351) on January 31st, 2023 for 31 occipital cortex samples using the Illumina HumanMethylation450 BeadChip. The ROSMAP, BDR, Weber PSP, and Gasparoni sorted datasets were consensus reprocessed through a unified QC pipeline for each brain tissue types separately using R (v4.1.2). For consistency with the Mayo Clinic discovery cohort, individuals with an indicated ethnicity other than non-Hispanic white or white were removed, individuals with unknown ethnicity were kept. Raw IDAT files were imported as a RGChannelSet object using read.metharray.exp in minfi 162. Sample and probe QC was performed using the following steps: (1) excluding probes that mapped to non-autosomal chromosomes; (2) excluding samples with a low (< 80%) bisulfite conversion rate estimated by bscon in wateRmelon 163; (3) using detection p values obtained by detectionP in minfi to exclude probes where < 99% of the samples had a detection p > 0.05; (4) removing samples where < 99% of the probes had a detection p > 0.05; (5) removing samples that failed sex check by clustering samples on the first PCs of previously identified robust sex-related probes on the X and Y chromosomes as implemented by estimateSex 164 in wateRmelon; (6) removing probes that have previously been identified to have cross-hybridizing potential or polymorphic targets (with European MAF > 0.01) 165, as well as the SNP and control probes; (7) removing outlier samples defined by the outlyx function in wateRmelon. Subsequently, to minimize intra- and inter-array variation, as well as probe design (I/II) bias, quantile normalization was performed using dasen 163 in wateRmelon, and outlier check was performed again on the normalized data. Individuals sex and age at death were compared when available from different sources and when discrepancies were identified, individuals were removed from the analysis. PCA was performed with the normalized β-values using prcomp(.scale=TRUE, center=TRUE) on the top 50% most variable CpGm sites based on standard deviation. Non-CpGm sites were removed and CpGm sites were lifted to the hg38 genome using the liftOver package and corresponding chain files provided by UCSC Genome Browser. rCpGms were calculated and data reduction steps were performed as described for the RRBS data. Intermediate Braak stages were grouped with the next lowest stage. AD cases were neuropathologically defined as individuals with a Braak stage ≥ 4 and a Thal phase ≥ 2. Control samples were defined as individuals with a Braak stage < 4 and a Thal phase < 2 (Supplementary Data 3). Any donor who was not neuropathologically defined as an AD case or control was designated as “Other”. For the Gasparoni sorted dataset, ADs were defined as having a Braak stage ≥ 4 and defined as AD by the original study 36 while controls were defined as those with a Braak stage ≤ 3 and defined as a control by the original study 36. For the Weber PSP dataset, PSP and control sample definitions were used as described 36,74. For ROSMAP and BDR, EWAS were performed in each dataset and tissue types separately for Braak stage, Thal phase, and AD diagnosis looking at AD vs Controls using the lmer function in the lmerTest R package. For the ROSMAP linear mixed models, we adjusted for age at death, sex, study (ROS/MAP), and batch as fixed effects as well as slide as a random effect variable. For the BDR linear mixed models, we adjusted for age at death, sex, and CpGm PC’s 1-3 as fixed effects as well as slide as a random effect variable. We included CpGm PC’s 1-3 in the BDR models because batch, although not available, has previously been found to be a source of variation in this dataset and captured by PC’s 32. PC’s 1-3 were chosen because these explained more than 5% of variation (PC1-23%, PC2-13%, PC3-6%) cumulating to a total of 41%. For the Gasparoni sorted dataset, EWAS were run using lmer and linear mixed models for AD diagnosis adjusting for age at death and sex as fixed effects as well as slide as a random effect. For Weber PSP, EWAS were run using lmer and linear mixed models for PSP diagnosis and Braak stage adjusting for age at death, sex, and CpGm PC’s 1-3 as fixed effects as well as slide as a random effect variable. Similarly, because we did not have information on batch, we included PC’s 1-3 explaining 25%, 12%, and 5% of the variation, respectively. Next, due to consistently high level of genomic inflation, we used the bacon R package 166 to adjust for inflation in each of the EWAS models (Supplementary Data 21).

Meta-analysis. Meta-analysis of Braak stage was performed in R with the meta::metagen function including rCpGm EWAS estimate and standard error statistics from the RRBS-TCX, BDR- DLPFC, and ROSMAP-DLPFC datasets. Test statistics from both fixed and random effect models are reported (Supplementary Data 11).

Previously published single nuclei/cell RNAseq studies. Oligodendrocyte specific gene expression association results from previously published single nuclei/cell RNAseq studies depicted in Fig. 5B and Fig. S12 were downloaded in October, 2024 from the following publications: Is 2024 75, Mathys 2023 76, Gamache 2023 77, Morabito 2021 78, Anderson 2023 79, Rexach 2024 80, Min 2023 81.

Cell type marker gene enrichment. One sided hypergeometric tests using the phyper(lower.tail=FALSE) function from the R stats package were used to test for enrichment of the top 100 and top 1000 cell type marker genes identified by BRETIGEA 65 in our dataset. We used the total number of genes expressed above background tested in our gene expression analysis (N=18,383 genes) as the total universe number.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The study presents a comprehensive bioinformatics approach to investigating DNA methylation (DNAm) in Alzheimer’s disease (AD) brains, with a particular focus on regional CpG methylation (rCpGm) and its association with tau-related pathology. The authors introduce a novel rCpGm method to improve biological relevance by grouping CpGm sites based on chromatin state annotations. The integration of epigenomic and transcriptomic data provides significant insights into the role of oligodendrocyte dysfunction in AD.\n\nThe bioinformatics analysis presented in this study is innovative, methodologically sound, and supported by robust statistical validation. The novel rCpGm method offers improved biological interpretability, and the integration of multi-omic data provides valuable insights into the role of DNA methylation in AD. However, addressing the concerns outlined below—particularly regarding cell-type specificity, causal inference, batch effects, and methodological reproducibility—would further strengthen the study’s impact and applicability.\n\nMajor Concerns\n1. Cell-Type Specificity of Methylation Signals\nGiven that bulk tissue methylation data were used, the primary conclusions regarding oligodendrocyte involvement may be influenced by cell composition effects. While adjustments for neuronal proportions were performed, additional validation using reference-free deconvolution methods (e.g., CIBERSORT, Epigenetic Clock models) or single-cell methylation data would provide greater confidence that the observed rCpGm associations are not confounded by cellular heterogeneity. The authors should clarify whether such methods were considered or could be integrated into future analyses.\n2. Causal Relationship Between Methylation and Gene Expression\nThe study demonstrates strong correlations between rCpGm and gene expression, particularly in oligodendrocyte-related genes. However, correlation does not imply causation. While the authors acknowledge this limitation, causal inference methods such as Mendelian Randomization, Bayesian networks, or mediation analyses could help determine whether DNAm changes drive gene expression alterations or vice versa. Implementing these approaches would further substantiate the functional role of methylation in regulating AD-related gene expression changes.\n\nMinor Concerns\n1. Validation of Regional Methylation Approach\nThe study utilizes a novel rCpGm method to cluster CpGm sites within biologically relevant chromatin regions. While this approach offers advantages over traditional single CpGm analyses, a direct comparison to other existing regional methylation aggregation methods (e.g., weighted co-methylation networks, tiling windows) is not provided. It remains unclear how the proposed method performs in terms of sensitivity and specificity relative to these alternatives.\n2. Technical and Batch Effect Corrections\nThe study integrates data from multiple platforms, including Reduced Representation Bisulfite Sequencing (RRBS), Illumina EPIC arrays, bulk RNA-seq, and single-nucleus RNA-seq. While batch effect corrections were performed, it is unclear whether cross-platform normalization techniques (e.g., ComBat, limma-based adjustments) were applied to ensure comparability across datasets. Given the potential impact of batch effects on methylation and expression levels, a more detailed discussion on this aspect would be beneficial.\n3. Longitudinal Validation of Methylation Changes\nThe study relies on cross-sectional data, which limits the ability to infer temporal relationships between methylation changes and disease progression. Future studies incorporating longitudinal methylation profiling would provide stronger evidence for the dynamic role of DNAm in AD pathogenesis.\n4. Reproducibility and Data Availability\nThe study claims high reproducibility across independent datasets (e.g., ROSMAP, BDR), which strengthens the validity of its findings. However, reproducibility also depends on the availability of full bioinformatics pipelines, processing scripts, and custom analytical methods. The authors should clarify whether their rCpGm approach has been made publicly available (e.g., via GitHub or Bioconductor) to allow other researchers to implement and validate the method in additional datasets.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Integrative Epigenomic Landscape of Alzheimer’s Disease Brains Reveals Oligodendrocyte Molecular Perturbations Associated with Tau\n\nThis manuscript presents a new method for clustering methylation cpg data in the brain, provides robust evidence of methylation alterations across the brain in Alzheimer’s disease and in relation to the major neuropathological features of Alzheimer’s disease, replicated the most robust signals in independent datasets, and provides evidence of gene expression consequences and some proposed biological pathways through the white matter that may underlie the observed associations. The paper is impressive in its scope and rigor and provides exciting new biological insights into the pathogenesis of Alzheimer’s disease and related dementias. There are some notable limitations that should be addressed, and I lay out my concerns by category (Major or Minor) below:\n\nMajor Concerns\n(1) The replication of individual signals appears to include nominal associations in the replication dataset. I would recommend relying fully on the correction for multiple comparisons for the number of targets that were brought to replication. FDR are reported in the figure, just lean on the ones that correct for replication. Perhaps the 93 are the FDR corrected, but it was challenging to follow.\n(2) It is unclear from the methods why only the top 50% were used for PCs during quality control and to what degree lane, participant, relatedness, or technical artifacts were accounted for in QC procedures. It seems like there is a substantial mean level difference across brain regions, I’m curious if some type of regression based correction like combatmet was needed or applied?\n(3) It seems that the cis-rCpG to gene analysis was performed across all samples, which would confound any association in which disease was driving both gene expression and methylation. I would recommend focusing the cis-QTL analysis on the normal cognition/path negative participants to remove this potential confound.\n\nMinor Concerns\n(4) I wonder if it would be worth using Hi-C and plaq-seq data from brain to try and infer some potential gene targets beyond the cis-QTL analysis provided?\n(5) The rationale for including CDK18 in the final list is weak and should be strengthened.\n(6) The analysis of the oligodendrocyte genes seems circular. The CpGs were identified in relation to tau, then the authors focus on genes that correlate with those same CpGs in the same dataset, and then argue that correlation of those same targets with tau demonstrates the relevance oligodendrocyte expression to tau.\n(7) The final step of replication of gene expression associations with tau pathology should also include a correction for multiple comparisons.\n(8) The causal order that is proposed could be better evaluated with either an MR approach (perhaps challenging if the methylation changes aren’t under genetic control) or a mediation approach. It seems equally likely that methylation is driving oligo expression resulting in a feed forward exacerbation of tau deposition/propagation, as it is that the tangles drive methylation difference resulting in white matter damage. Without evidence of white matter damage in these participants, the causal argument is challenging to follow or justify in the current manuscript. I’m not sure what my recommendation is, I just wasn’t very convinced by the argument around order of events.", "reviewer_4": "What are the noteworthy results?\nThe authors, using temporal cortex and cerebellum samples from over 400 donors, provide a comprehensive exploration of RRBS analysis of individual CpG sites and chromatin-based regional analysis for combining CpG sites. These two methods are then associated with several AD-related endophenotypes, including amyloid and tau. Both these methods of analysis identify several genomic locations with altered methylation, with the majority of FDR significant results being in tau measures. These results provide overlap with previous results on brain pathological array-based methods. Using ROSMAP and BDR, which have array-based DNA methylation measures pre-frontal cortex and occipital cortex, the authors show consistency in identified sites and chromatin regional analyses, although these associations use Braak, Thal and add cerebral amyloid angiopathy score.\nInvestigating significant findings, the authors show oligodendrocyte markers and myelination-related genes using gene-set enrichment. They also show a correlation between regional and single CpG methylation and gene expression and identify mQTLs from tau-associated sites.\n\nWill the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not original, please provide relevant references.\nThe work is of significance to the field of AD research and the epigenetics of AD. It shows good overlap with previous analyses into similar methods and tissues, but also provides an alternate regional analysis for DNA methylation sites with potentially more biological relevance. The work is original in that it uses a new regional DNA methylation method and applies it to previously published data, showing its applicability and validity. This method also allows for the identification of new regions and genes of interest for integrative genomic analyses.\n\nDoes the work support the conclusions and claims, or is additional evidence needed?\nI believe the work supports the conclusions given.\n\nAre there any flaws in the data analysis, interpretation and conclusions? Do these prohibit publication or require revision?\nThe data analyses seem appropriate in method, and appropriate covariates, multiple testing adjustments, and post-hoc analyses have been used. My only revision is the ability to find the full script. The authors state that scripts are available via the AD Knowledge Portal, however, the portal is sometimes difficult to navigate, and I was unable to find them. This may be due to my lack of familiarity with the portal, but a direct link or instruction is needed for readers who face a similar problem.\n\nIs the methodology sound? Does the work meet the expected standards in your field?\nSamples included are well described and appropriate to the analyses. Methodology is sound and meets the standards in the field. Work has been subjected to power calculations for identifying of appropriate effect sizes and using the sample sizes given. Appropriate use of Spearman’s and Pearson’s correlation for normal and non-normal/categorial measures is used. Multiple testing using FDR is used. CpG ~ AgeAtDeath + Sex + genetic PC 1-3 is used for EWAS analyses, which seems appropriate. Similarly, the methodology for gene expression association and mQTLs seems appropriate. Post-hoc analyses are carried out as per standards.\n\nIs there enough detail provided in the methods for the work to be reproduced?\nOnce access to the code is made easier and can be checked, yes.\nMost recommendations for adjustments are minor.\n\nMajor\nThe website doesn’t allow access as it requires a login. Please make it accessible for the exploration of data.\nLine 890 – Please provide a direct link or details on where to find the code on the knowledge portal.\nMinor\nA lot of these are writing or consistency errors which may be present in other points.\nLine 14 (and others) – 6 to six and 2 to two\nFigure 1 – Define TCX in the acronym list in the legend\nLine 75 – remove “CpG sites (CpGm)” as already defined in intro\nLine 223-229 and Fig S6 – Although I agree that the correlation is there, the distribution of beta values mean the data is often skewed due to the ends of the association, especially in very large datasets. I would suggest using a random sampling or subsetting of data (especially in the RRBS vs. array correlations).\nLine 281-293 – Please include a table for GSEA results in supplementary and refer in the text\nLine 596 – four to 4, two to 2\nTable S1 (and others) – Define acronyms (TCX, CER) in legend\nFigure S1 - stdev to SD for consistency with Table 1\nLine 619 – Genome wide to Genome-wide for consistency\nLine 622-3 – Details or reference on the QC pipeline should be provided\nLine 622-623 – QC is used for quality control and quality controlled. Suggest changing line 622-3 to be “and underwent rigorous quality control (QC)”.\nLine 624 – R2 to R^2 (meaning superscript, limit of the text box)\nLine 664 – quality control to QC\nLine 829 – non-Hispanic\nLine 851-2 – “defined as individuals with a Braak stage ≥ 4 and a Thal phase ≥ 2. Control samples were defined as individuals with a Braak stage ≤ 4 and a Thal phase ≤ 2” from reading table S3, this should be Braak stage < 4 and Thal phase < 2.\nFigure S2 and others - Acronyms for CER, TCX, CpGm and rCpGm need defining\nFigure S3 – For a more direct comparison between datasets, it would be good to include the current study data with BDR and ROSMAP to more easily compare the distributions. RRBS is going to cover more sites, but it would be interesting to see the array selection of sites and how they compare. As stated, arrays have a narrower range of values, so it would be interesting to see this clearly. ROSMAP shows higher average methylation in rCpG analyses (FigS3). Is this for a reason, or a concern with comparison, or proposed to be an artifact if the differences in the arrays (in probe number) and site selection?\nTable S4 and S5 - Add gene annotations to the table for easier interpretation and alignment with the text (line 156-163).\nFigure S5 - Several plots appear to only have one or very few data points, so correlation is not assessable. These should be removed and/or explained.\nFigure 4 – Very busy and a lot of information, especially with the multiple outcomes being associated and is difficult to parse all the information, especially since some measures don’t overlap. This is especially difficult since association in the associations in some of the cohorts defining tiers is from one measure.\nFigure 5 – Text at A4 size is too small to see clearly and image looks blurry until very zoomed in. This figure needs to have text increased or possibly removing some information to make it viewable for a main figure. Advice may be to keep part C or A and B and move the others to supplementary if figure limit has been reached."}, "author_response": "We thank the Editors and the Reviewers for their careful review of our manuscript, their comments, input, and suggestions. We are grateful for all the positive comments from the Reviewers. Below, we provide a point-by-point response to all comments. We believe that the modifications to address the comments of our reviewers have further enhanced the impact of our study.\n\nSpecific points provided by the reviewers are highlighted in italic font and any text edits to the manuscript are provided in blue Times New Roman font. In our submission, we provide two versions of our Main and Supplementary files; one in which all changes are tracked and another one in which all changes are accepted. Page and line numbers align with the tracked changes manuscript version. Table and figure numbers have been updated throughout to reflect new additions.\n\nReviewer 1\nThe first reviewer summarized the study describing “The study presents a comprehensive bioinformatics approach to investigating DNA methylation (DNAm) in Alzheimer’s disease (AD) brains, with a particular focus on regional CpG methylation (rCpGm) and its association with tau-related pathology. The authors introduce a novel rCpGm method to improve biological relevance by grouping CpGm sites based on chromatin state annotations. The integration of epigenomic and transcriptomic data provides significant insights into the role of oligodendrocyte dysfunction in AD.\n\nThe bioinformatics analysis presented in this study is innovative, methodologically sound, and supported by robust statistical validation. The novel rCpGm method offers improved biological interpretability, and the integration of multi-omic data provides valuable insights into the role of DNA methylation in AD. However, addressing the concerns outlined below—particularly regarding cell-type specificity, causal inference, batch effects, and methodological reproducibility—would further strengthen the study’s impact and applicability.”\nWe thank the reviewer for their evaluation of our study, positive comments and suggestions for further improvement.\nThe reviewer has noted two major concerns (1.1 and 1.2) and four minor concerns (1.3-1.6) to which we have responded as detailed below.\n\n1.1. The reviewer comments on the “Cell-Type Specificity of Methylation Signals” stating “Given that bulk tissue methylation data were used, the primary conclusions regarding oligodendrocyte involvement may be influenced by cell composition effects. While adjustments for neuronal proportions were performed, additional validation using reference-free deconvolution methods (e.g., CIBERSORT, Epigenetic Clock models) or single-cell methylation data would provide greater confidence that the observed rCpGm associations are not confounded by cellular heterogeneity. The authors should clarify whether such methods were considered or could be integrated into future analyses.”\nWe agree that the methylation signals could be impacted by cell type composition, which is why, as the reviewer notes, we performed sensitivity analyses of our top hits adjusting for estimated neuronal proportions. As we identified a potential oligodendrocyte specific signal, we have also now adjusted for estimated oligodendrocyte cell proportions in our top hits. Using the same methods described in our Methods section, we estimated oligodendrocyte proportions by leveraging the bulk gene expression data and DSA deconvolution1 using the top 100 BRETIGEA2 cell type markers for oligodendrocytes. Similar to our neuronal proportion adjustments, we included oligodendrocyte proportion estimates in our main model (rCpGm ~ Phenotype + Age at Death + Sex + genetic PC1-3 + Oligodendrocyte proportion estimates), finding highly consistent test statistics suggesting that our findings are unlikely to be due to differences in oligodendrocyte proportions (Rebuttal Figure 1).\n\nRebuttal Figure 1: Comparison of Primary EWAS model results with EWAS model results adjusting for estimated oligodendrocyte proportions for all FDR significant TCX rCpGm associations. Comparison of test statistics (estimates and p-values) for all FDR significant rCpGm associations from the primary rCpGm EWAS models and models adjusting for oligodendrocyte proportion estimates. To test the correlation between estimates and p-values, we combined top hits from all phenotypes and performed Spearman’s rank correlation tests which are reported at the top of each panel. Primary model statistics are the x-axes, sensitivity model statistics are the y- axes. Green dashed line is the –log10(p-value = 0.05) value.\nWe have included Rebuttal Figure 1 as a new panel (D) in Figure S5 along with the other sensitivity model comparisons.\nWe have also edited the text as follows:\nResults text on page 14 was changed from: “To ensure that these FDR significant rCpGm associations are robust to any biological and technical variations, we repeated these analyses after further adjustments for Braak stage and Thal phase, neuronal cell type proportion estimates, and sequencing batch.” To:\n“To ensure that these FDR significant rCpGm associations are robust to any biological and technical variations, we repeated these analyses after further adjustments for Braak stage and Thal phase, neuronal cell type proportion estimates, oligodendrocyte cell type proportion estimates, and sequencing batch.”\nMethods text on page 41 from: “Neuronal cell type proportions were estimated from the bulk tissue RNAseq data using the Digital Sorting Algorithm (DSA)1 via the DSA R package (v1.0).” To:\n“Neuronal and oligodendrocyte cell type proportions were estimated from the bulk tissue RNAseq data using the Digital Sorting Algorithm (DSA)1 via the DSA R package (v1.0).”\nMethods text on page 44 from: “We performed post-hoc sensitivity analyses to test the robustness of the top associations to effects of batch, neuropathological changes, and neuronal cell type proportions.” To:\n“We performed post-hoc sensitivity analyses to test the robustness of the top associations to effects of batch, neuropathological changes, neuronal and oligodendrocyte cell type proportions.”\nWe also added the following sentence to the Methods on page 44:\n“To determine if top associations were influenced by oligodendrocyte cell type proportions, we included DSA estimated oligodendrocyte cell type proportions as a fixed-effect covariate in the linear regression model: rCpGm ~ Endophenotype + Age at death + Sex + Genetic PCs 1-3 + Oligodendrocyte Cell Proportion.”\nAdditionally, we investigated the association between the 15 oligodendrocyte-tau-related rCpGms and total tau and pTau adjusting for the estimated oligodendrocyte proportions (Rebuttal Table 1) finding largely the same test statistics.\n\nRebuttal Table 1: Association of 15 oligodendrocyte- tau-related rCpGms with pTau and total tau levels from our primary model and a model adjusting for estimated oligodendrocyte proportions.\n\nAligning with the reviewer’s suggestion, we also investigated the 15 rCpGms of interest in cell type sorted DNAm data generated by Gasparoni et al 20183. Gasparoni and colleagues sorted 31 occipital cortex tissue samples provided by the Neurobiobank Munich using NeuN immunolabeling followed by fluorescence associated cell sorting (FACS) obtaining a neuronal nuclei population (NeuN+) and a non-neuronal/glial population (NeuN-). We downloaded the available raw Illumina 450K IDAT files from GEO (GSE66351) and processed these independently for each cell type population as described in our Methods for the Brains for Dementia Research (BDR), ROSMAP, and Weber et al 2018 datasets. We obtained a final quality-controlled dataset that included 26 glial samples (13 AD, 11 controls, 2 other) with 137,186 rCpGms and 27 neuronal samples (13 AD, 12 controls, 2 other) with 137,169 rCpGms (Rebuttal Table 2).\n\nRebuttal Table 2: Demographics of Replication Dataset. Summary demographics of the Gasparoni et al. 2018 sorted cell-type dataset post-quality control.\n\nOf our 15 oligodendrocyte-tau-related rCpGms of interest, 14 were measured in this dataset. Association models between ADs vs Controls were run using mixed-effect linear regression models adjusting for age at death and sex as fixed effects and slide as a random effect variable in the glial and neuronal cell-type populations. We found that all but one of our rCpGms of interest in the glial cell population positively associated with AD matching what we expect based on our own data while only 6 out of the 14 rCpGms in the neuronal cell population positively associated with AD (Rebuttal Figure 2). While none reached a nominal statistical significance likely owing to the very small sample size of this sorted cell-type methylation dataset, there are three rCpGms with a trending association level (p-value < 0.07) in the glial enriched population with consistent direction of effect as our study (SOX8, KLK6 rCpGms). In addition, the regression estimates of the glial enriched population are larger in magnitude than the neuronal enriched populations which parallel previous evidence that AD-related DNAm associations are likely a result of changes in non-neuronal rather than neuronal populations4. Together, these results suggest that the rCpGm associations we have identified in our dataset likely originates from the glial cells which largely consists of oligodendrocytes rather than the neuronal cells.\n\nRebuttal Figure 2: Association of AD diagnosis with 15 oligodendrocyte-related rCpGms of interest in Illumina 450k DNAm data derived from glial (NeuN-) and neuronal (NeuN+) enriched occipital cortex nuclei populations from Gasparoni et al. 2018. After quality control, there were 13 AD and 11 Control NeuN- samples and 13 AD and 12 Control NeuN+ samples. Within each tile is the regression estimate (top) and p-value (bottom) labeled with ‘^’ if the p-value is < 0.1. Tiles are colored by regression estimate direction where red is a positive association with AD diagnosis and blue is a negative association.\nWe included Rebuttal Figure 2 in the supplemental figures as Figure S11 and Rebuttal Table 2 in the supplemental tables within Table S3. We have also included the following text in the Results section on page 21-22:\n“Further investigation of the 15 rCpGms in Illumina 450k DNAm data derived from NeuN+ and NeuN- sorted occipital cortex nuclei3, revealed 14 rCpGms with positive associations with AD in the glia enriched NeuN- population while only 6 rCpGms had positive associations with AD in the neuronal enriched NeuN+ population (Figure S11). The regression estimate magnitudes of the glia enriched NeuN- population were also larger than the NeuN+ population suggesting that the associations identified in the bulk DNAm data are likely originating from DNAm changes in the glial cells rather than neuronal cells4.”\nWe also included Methods text describing the quality control and analyses for these new results as follows:\nPage 45: “Sorted nuclei data from Gasparoni et al. 20183 was downloaded from GEO (GSE66351) on January 31st, 2023 for 31 occipital cortex samples using the Illumina HumanMethylation450 BeadChip.”\nPage 46: “For the Gasparoni sorted dataset, ADs were defined as having a Braak stage ≥ 4 and defined as AD by the original study while controls were defined as those with a Braak stage ≤ 3 and defined as a control by the original study3.”\nPage 47: “For the Gasparoni sorted dataset, EWAS were run using lmer and linear mixed models for AD diagnosis adjusting for age at death and sex as fixed effects as well as slide as a random effect.”\nLastly, although the sorted nuclei population generated by Shireby et al 20224 is not publicly available, to the best of our knowledge, the authors did make available the estimated NeuN+, SOX10+, and NeuN-/SOX- cell type proportions for their bulk DNAm data. Shireby and colleagues used their sorted DNAm data as a reference and the Houseman5 method incorporated in the R package minfi6,7 to deconvolve cell proportions from the bulk dorsolateral prefrontal cortex (DLPFC) and occipital cortex (OCC) data. Using these cell type proportion estimates, we investigated how further adjusting for Neuronal (NeuN+) and Oligodendrocyte (SOX10+) estimated proportions would impact the results of our BDR association models. We compared test statistics from our primary model (rCpGm ~ phenotype + age at death + sex + PC1-3 +(1|slide)) with and without bacon correction to a model that also included SOX10+ and NeuN+ cell proportion estimates for the 15 oligodendrocyte-tau related rCpGms. We found highly similar test statistics across all 3 models for each phenotype (AD diagnosis, Braak stage, and Thal phase) in the DLPFC and OCC (Rebuttal Figure 3). Using our definition for replication, even after adjustment for these estimated cell proportions, we would have still included all 15 rCpGms in downstream analyses. This is yet another piece of evidence that our findings are replicable in other datasets and unlikely to be a consequence of changes in cell proportions. Rather, they reflect methylation changes predominantly occurring in oligodendrocytes and associated with tau phenotypes in AD brains.\n\nRebuttal Figure 3: Association of AD traits with the 15 oligodendrocyte-tau-related rCpGms of interest in the BDR dataset with and without adjustment for estimated NeuN+ and SOX10+ cell type proportions. Within each tile is the regression estimate and (p-value) colored by p-value where a dark color indicates a p-value < 0.05 and light color p-value > 0.05. Red indicates a positive association and blue a negative association. The ‘primary.bacon’ model aligns with the results analyzed in our manuscript after applying bacon while the ‘primary’ model reflects results prior to bacon adjustment, both are shown to allow for comparisons with the cell type proportion models where bacon was not applied.\nTo the best of our knowledge, the data generated by Gasparoni et al. is the only cell type specific dataset currently available to us from AD donors which we could investigate our rCpGms of interest. While we believe that cell type heterogeneity does not largely impact the findings we present in our manuscript, we do agree with the reviewer that future investigation of our results in cell-type specific DNAm data should be performed to fully determine the potential contributions of cell type heterogeneity. We have edited the following sentence in the Discussion on page 33 from:\n“Although the associations for the top hits we identified remained even after adjusting for neuronal proportions, future single cell type specific DNAm studies, particularly those from glial cell types4, will be important in parsing out potential cell type specific rCpGm effects in these associations.” To:\n“Although the associations for the top hits we identified remained even after adjusting for neuronal or oligodendrocyte proportions and showed consistent and larger magnitude associations in glial enriched nuclei, future single cell type specific DNAm studies, particularly those from glial cell types4, will be important in parsing out potential cell type specific rCpGm effects or the role of cellular heterogeneity in these associations.”\n\n1.2. The reviewer further comments on “Causal Relationship Between Methylation and Gene Expression” that “The study demonstrates strong correlations between rCpGm and gene expression, particularly in oligodendrocyte-related genes. However, correlation does not imply causation. While the authors acknowledge this limitation, causal inference methods such as Mendelian Randomization, Bayesian networks, or mediation analyses could help determine whether DNAm changes drive gene expression alterations or vice versa. Implementing these approaches would further substantiate the functional role of methylation in regulating AD-related gene expression changes.”\nWe thank the reviewer for their suggestion to more formally investigate causation. First, we applied the Mendelian Randomization (MR) approach using SNPs identified through our mQTL as instrumental variables (IV). Of the 15 tau-oligodendrocyte-related rCpGms we identified, only 4 rCpGms significantly associated (FDR p-value ≤ 0.05) with nearby genetic variants and so could be leveraged for MR. These SNPs included 92 associated with the MBP-related rCpGm, 3 with RTKN, 3 with ST18, and 1 with KLK6. As IVs need to be independent of each other to meet MR assumptions, we used PLINK28 to clump these 99 SNPs defined by linkage disequilibrium (R2 < 0.2) resulting in 8 independent loci, 4 with the MBP-related rCpGm, 2 with RTKN, 1 with ST18, and 1 with KLK6 (Rebuttal Table 3). We next tested the prediction ability of each IV SNP set for their respective rCpGm by including the rCpGm IV SNP set in a linear regression model, finding all F-statistics were greater than the suggested threshold of ten9 (Rebuttal Table 3). We also investigated if the 8 independent SNPs associated with expression of their related gene, finding no significant associations (p-value > 0.05). Together, the F-statistic and lack of association to gene expression suggest that these IV sets are valid for MR testing.\n\nRebuttal Table 3: SNPs used for Mendelian Randomization as Instrumental Variables.\n\nNext, we performed an MR analysis using the TwoSampleMR10 R package to investigate if there is evidence that these four rCpGms may cause changes in expression of their related genes. For each of the 4 rCpGms, we did not find significant evidence of causation (p-value > 0.05) (Rebuttal Table 4).\n\nRebuttal Table 4: Mendelian Randomization results for rCpGms causing changes in related gene expression.\n\nWhile these MR results do not support causation between the 4 rCpGms and related gene expression levels which were testable, it is important to note that our ability to detect a causal relationship through MR is dependent upon the prediction strength of the IV to the exposure (rCpGm), validity of MR assumptions, and statistical power. While all IV sets had an F-statistic > 10, this does not definitively confirm an IV is strong but rather indicates that the IV may be of sufficient strength as bias has been reduced to an acceptable level11. Additionally, the power of an MR test depends on sample size and the proportion of variance in the exposure that the IV explains12. The IV sets we were able to investigate overall explained a small proportion of variance in the rCpGms (3% - 9.4%, Rebuttal Table 3). This suggests that while the SNPs have predictive ability, there may not be enough power to detect a causal association with these IVs as MR studies typically require thousands of samples to detect evidence of causation. This is supported by a power analysis which revealed that we likely lack power to detect a significant MR signal in our dataset of 455 samples because if the true magnitude of association between MBP and related rCpGm was at the highest end of the 95% confidence interval (CI) (β = -0.98), we would need 694 samples to have 80% power to identify a causal association (Rebuttal Table 5). More than likely, the true magnitude of association at the MBP locus is somewhere in between the 95% CI and if it was closer to the magnitude estimated by our model (β = -0.72), a sample size > 1,200 would be required. This is similarly true for the other three testable loci as well where we would need a sample size >1,000 to have 80% power to identify causal evidence (Rebuttal Table 5).\n\nRebuttal Table 5: Power calculation for Mendelian Randomization of rCpGms causing gene expression changes.\n\nFurthermore, a key assumption for MR is that the IVs are independent of confounders related to the exposure and outcome as well as the outcome itself, in this case gene expression. While the IV SNPs do not associate, even nominally, with expression levels of our genes of interest, pTau , total tau , nor AD diagnosis14, we cannot rule out the possibility that unknown confounders are influencing this relationship invalidating the IVs.\nAs a result, we next performed mediation analysis to determine if there was evidence of causal mediation between the oligodendrocyte-tau-related rCpGms, oligodendrocyte gene expression and pTau or total tau . We performed two types of mediation analyses using the mediation15 R package, one testing if gene expression mediates the association between rCpGm and tau (Rebuttal Figure 4A), and the second testing if rCpGm mediates the association between tau and gene expression (Rebuttal Figure 4B).\n\nRebuttal Figure 4: Diagram of causal mediation schemes tested. (X) is the exposure variable, (M) is the mediator, and (Y) is the outcome variable. GE = Gene expression\nPrior to performing mediation, we standardized through z-score transformation the rCpGm and gene expression levels to allow for comparisons of the mediation test statistics. Tau measures were already standardized. We adjusted for age at death and sex in all input linear regression results and estimated confidence intervals through bias-corrected and accelerated (bca) bootstrapping with 1,000 simulations. We also performed sensitivity analyses including estimating how robust these associations would be to different types of unmeasured confounders.\nIn our first analysis (Rebuttal Figure 4A), we evaluated 32 mediation schemes for the 16 oligodendrocyte-tau-related rCpGm-gene pairs with pTau and total tau of which, 7 showed statistically significant mediation effects (p-value ≤ 0.05). These included the KLK6, ST18, CDK18, CLMN, MBP, and CD22 rCpGm-gene pairs with pTau and CD22 with total tau with a significant proportion mediated (p-value ≤ 0.05) ranging from 7% to 24% suggesting these may reflect partial mediation (Rebuttal Table 6). In our second mediation scheme (Rebuttal Figure 4B), while all 32 schemes met modern requirements of mediation, only 13 met previous requirements that the total effect of X on Y needs to be significant. Because we have good theoretical basis for this scheme, we used the modern requirements finding 31 schemes showed statistically significant mediation effects with 12 of these having a significant proportion mediated ranging from 15% to 50%. These included the KLK6, CDK18, LDB3, CLMN, MBP, and CD22 rCpGm-gene pairs with pTau and CDK18, LDB3, CLMN, MBP, and CD22 with total tau. Next, we performed sensitivity analyses finding there were no significant interactions between the mediator (M) and exposure (X) on the outcome (Y). We additionally assessed the robustness of the average mediated and direct effects to unmeasured confounders through the medsens function finding that although the direct effects are often more robust to unmeasured confounders, the mediated effects are more likely to be affected by unmeasured confounders16,17 (Rebuttal Tables 8 and 9).\n\nRebuttal Table 6: Mediation Results for scheme A (Rebuttal Figure 4A).\nRebuttal Table 7: Mediation Results for scheme B (Rebuttal Figure 4B).\nRebuttal Table 8: Sensitivity results for mediation analyses performed for scheme A.\nRebuttal Table 9: Sensitivity results for mediation analyses performed for scheme B.\n\nGiven that many of the mediated effects we identified were significant in both scheme A and B, we hesitate to draw conclusions of causation from this analysis. First, while we adjust for common confounders like age at death and sex, we cannot entirely rule out the possibility of unmeasured confounders which our sensitivity analyses suggest the mediated effect estimates would be sensitive to. Second, while casual mediation analysis is helpful in cases where the data is from a randomized-control trial, is longitudinal, or there is definitive knowledge about temporal ordering like genetic variants (for the most part) being fixed at birth, our data is cross-sectional where temporal ordering is not definitive. In cases where temporal ordering is not definitive, causal mediation analysis is not able to reveal whether a certain temporal order is true over another and its application is ill-advised18-20.\nTaken together, we are not able to definitively determine a causal relationship between the identified rCpGms and expression of nearby genes. While previous literature does support this is a likely relationship (i.e. methylation changes lead to gene expression changes), the only definitive way to conclude that these rCpGms are causing changes in expression of nearby genes is to functionally test these through experimentation.\nAs a result, we have included the following text in the Results section on page 24:\n“However, it is equally plausible that other untested genetic or environmental factors lead to DNA hyper-methylation, subsequently leading to downregulated myelin genes and oligodendrocyte dysfunction, culminating in brain accumulation of toxic tau forms21.”\nIn the Discussion on page 29, we included the following sentence:\n“Nonetheless, future functional experiments will need to verify the casual direction of effect between these rCpGms and their related genes.”\nAs well as:\n“However, other untested factors may also be driving CpG hypermethylation, resulting in oligodendrocyte gene expression reductions and dysmyelination, which culminates in toxic tau propagation84.”\nWe edited Discussion text on page 31 from: “Our epigenome-wide association study of deep brain AD protein biochemical levels, coupled with our novel rCpGm approach, and integrated transcriptome and epigenome data from >4,000 collective brain samples demonstrate strong, replicable and biologically consistent associations between brain tau and hypermethylation near oligodendrocyte genes likely leading to their lower expression, dysmyelination and neurodegeneration in AD and other tauopathies (Figure 5C).” To:\n“Our epigenome-wide association study of deep brain AD protein biochemical levels, coupled with our novel rCpGm approach, and integrated transcriptome and epigenome data of donors from >4,000 collective brain samples demonstrate strong, replicable and biologically consistent associations between brain tau and hypermethylation near oligodendrocyte genes possibly leading to their lower expression, dysmyelination and neurodegeneration in AD and other tauopathies (Figure 5C). However, functional experimental studies are needed to decipher the actual sequence of events in this cascade.”\nWe edited Discussion text on page 32-33 from: “Our results uncovered opposite roles for brain total tau (likely beneficial) and ptau (detrimental) and support a model where brain tauopathy influences rCpGm resulting in hypermethylation and reduced expression of oligodendrocyte genes culminating in dysmyelination.” To:\n“Our results uncovered opposite roles for brain total tau (likely beneficial) and pTau (detrimental) and support two hypothetical models, one where brain tauopathy influences rCpGm resulting in hypermethylation and reduced expression of oligodendrocyte genes, and another where reduced expression of oligodendrocyte genes from rCpGm hyper-methylation disrupts tau homeostasis, both of which could culminate in dysmyelination.”\nFinally, in combination with reviewer point 1.5, we have included the following sentence in the Discussion on page 33:\n“In addition, the results of this study are based on cross-sectional data at the time of autopsy and future longitudinal blood or CSF DNAm studies and functional experimentation investigating our top hits will be important in understanding the temporal relationship between DNAm, gene expression changes and tau pathology.”\n\n1.3. A minor concern of the reviewer was the “Validation of Regional Methylation Approach” stating that “The study utilizes a novel rCpGm method to cluster CpGm sites within biologically relevant chromatin regions. While this approach offers advantages over traditional single CpGm analyses, a direct comparison to other existing regional methylation aggregation methods (e.g., weighted co-methylation networks, tiling windows) is not provided. It remains unclear how the proposed method performs in terms of sensitivity and specificity relative to these alternatives.”\nWe thank the reviewer for the opportunity to compare our method with other existing regional aggregation methods.\nWe compared our rCpGm method to results from comb-p22, a post-hoc sliding window approach often used in DNA methylation studies23-31. Comb-p groups CpGs into regions based on their proximity within a defined window and association significance. After applying the comb-p method to our individual CpGm EWAS results, we identified 2,362 significant associations (Sidak p-value ≤ 0.05) with 1,833 unique comb-p regions (Rebuttal Figure 5). These comb-p regions mainly followed similar trends as our individual CpGm EWAS and rCpGm EWAS where total tau had the highest number of associations (70.2%) followed by pTau (21.6%) and APOE (6.4%). Looking across measures, we similarly saw the largest overlap of significant associations between total tau and pTau , all with opposing directions of association, also consistent with our findings. While comb-p does not calculate an effect direction of a region, we defined positive and negative comb-p associations by the direction of associations for the majority of the individual CpGs within that region.\n\nRebuttal Figure 5: Upset Plot of significant CpGm association regions identified with Comb-p.\nNext, we compared the comb-p results to our rCpGm results to determine regions either common to both or uniquely identified by one of the methods. We defined regions as being the same if the comb-p regions overlapped at any point with the rCpGm regions. This definition resulted in a total of 5,391 associations in the rCpGm method and 2,305 associations in the comb-p method. Of these, there were 1,628 associations that were common and significant in both methods (Rebuttal Table 10), including 1,336 with total tau , 288 with pTau , and 4 with APOE . All of these associations except a single one with total tau had the same direction of effect detected by both methods. On the other hand, there were 677 associations uniquely identified by comb-p and 3,763 by our rCpGm method.\n\nRebuttal Table 10: Comparison between associations identified using the comb-p method vs. rCpGm method.\nWe further investigated if the 93 rCpGms we found to be significantly associated with both total tau and pTau and that showed replication in at least one of the BDR or ROSMAP datasets were identified by comb-p as these are high confidence, replicated associations. We found 33/93 regions (35%) were identified in both the total tau and pTau comb-p analyses including an additional 7 with pTau alone (40 total) and additional 24 with total tau alone (57 total).\nLastly, we investigated the 15 oligodendrocyte-tau-related rCpGms we identified. Of these 15 rCpGms, eight were present in the comb-p pTau and 13 with total tau analyses, of which only six (40%) were associated with both tau measures. These six included the rCpGms related to expression of the CDK18, TGFA, RTKN, two with SOX8, and MBP genes.\nWhile we identified similar patterns of associations across measures, it is important to note that there are major weaknesses to the comb-p method when compared to our novel rCpGm method. First, comb-p does not take into account direction of effect when determining significant regions which allows individual CpGs with different effect directions to be grouped into the same region. We found that 6.6% of our comb-p associations had individual CpGms with opposing effect directions grouped into the same region. Comp-p’s lack of accounting for effect direction therefore requires the investigator to perform additional interrogation of significant regions. Moreover, because comb-p does not account for effect direction, biologically meaningful comparisons of effect direction, as we have done with our rCpGm data, cannot be conducted for comb-p regions of association.\nFor the purpose of comparing comb-p with our rCpGm results, we defined positive and negative comb-p associations by the effect direction for the majority of the individual CpGs making up that region. Some of the comb-p regions with discordant underlying CpGms had a slim difference in the proportion of CpGs with opposing effect direction. Thus, lack of having a direction of effect can lead to relatively arbitrary assignment of direction, again reducing biological relevance. Additionally, comb-p does not determine the magnitude of effect which limits the interpretations and conclusions that can be drawn using this method. Importantly, comb-p does not incorporate data from every CpGm within a region, but rather it only includes CpGms within a region that collectively pass a significance threshold. Defining regions through this method potentially allows for the overemphasis of some CpGm’s contributions to an association while completely ignoring others. Of the 2,362 significant comb-p regional methylation associations, 255 had between 1 to 125 CpGs (average = 10.9) within the defined regional window that were not incorporated in the comb-p results.\nAnother limitation of comb-p is that its CpGm grouping is not based on any biological scaffold, but rather based on arbitrarily defined windows and statistical thresholds. This poses a significant drawback when comparing associations of the same methylation data with different phenotypes from the same cohort or when conducting comparisons of different methylation datasets. Comb-p neither identifies biologically based regional associations nor enables a head-to-head comparison of the same regional associations. In contrast, our rCpGm method groups CpGms into regions based on biological definitions. Hence, by using our rCpGm method, the exact same region of CpGms is tested across phenotypes and even across different methylation datasets. Thus, our rCpGm method enables testing of biologically meaningful regions yielding results that can be directly compared and interpreted in the same manner across datasets.\nIn summary, results of our comparison suggest that comb-p identified associations detected by our rCpGm method, especially for total tau and pTau associations. However, our rCpGm method has numerous strengths over this method as follows: rCpGm improves the detection of regional signals across EWAS, has many more unique associations including those that replicate in additional, independent datasets, takes into account all CpGms within a region, and conducts biologically meaningful analyses across consistently defined methylation regions thereby enabling a direct head-to-head comparison of results across datasets.\nWe have edited and included the following text on our comb-p comparison in the Results section:\nOn page 13, we added: “We also compared results from our novel regional CpGm association method (rCpGm) to the published post-hoc regional aggregation comb-p method22.”\nAdded on page 15-16: “To compare our rCpGm approach to the previously published comb-p method22, we applied comb-p which groups proximal CpGms into regions based a sliding window approach and association significance. We identified 2,362 significant associations (Sidak p-value ≤ 0.05) in 1,833 unique regions by comb-p (Figure S6, Table S8). Unlike our rCpGm approach, comb-p does not take into account direction of effect when determining significant regions. Further, since comb-p defines regions based on arbitrary windows and statistical significance, different regions are identified for each endophenotype EWAS. This hinders a direct comparison of comb-p regions across endophenotypes, a feature which is possible by our rCpGm approach. To resolve the comb-p shortcomings on effect direction and regional variability, we defined positive and negative comb-p associations by the effect direction for the majority of the individual CpGs making up that region. We considered regions across different endophenotype EWAS to be common if there was any genomic overlap between these regions. Comb-p based regional associations largely followed similar trends as the individual CpGm and rCpGm EWAS where total tau had the highest number of significant associations (70.2%) followed by pTau (21.6%) and APOE (6.4%). Additionally, the largest overlap of significant associations between endophenotypes were for total tau and pTau , all with opposing directions of association, consistent with rCpGm results. There were 1,628 regional methylation associations common to comb-p and rCpGm, however rCpGm identified 3,763 additional associations whereas comb-p uniquely identified 677 associations.”\nOn page 16 edited from: “Although both the individual CpGm and our novel rCpGm analyses found the most associations with brain tau biochemical endophenotypes, there were twice as many FDR significant loci using the latter approach. Of the 2,489 FDR significant individual TCX CpGms, 68% were included in at least one of the FDR significant rCpGms. Notably, 1,671 individually significant CpGms were part of the 1,000 rCpGms associated with Total tau , and 26 such CpGms were within 23 rCpGm associated with pTau . Thus, our novel rCpGm approach identifies both those loci captured by the more traditional single CpGms but importantly also many others that are missed by the latter.” To:\n“Although both the individual CpGm followed by comb-p regional grouping and our novel rCpGm analyses found the most associations with brain tau biochemical endophenotypes, there were twice as many FDR significant loci using the latter approach. Of the 2,489 FDR significant individual TCX CpGms, 68% were included in at least one of the FDR significant rCpGms. Notably, 1,671 individually significant CpGms were part of the 1,000 rCpGms associated with Total tau , and 26 such CpGms were within 23 rCpGm associated with pTau . Taken together, our novel rCpGm approach identifies both those loci captured by the more traditional single CpGms followed by comb-p grouping but importantly also many others that are missed by the latter.”\nWe also edited text in the Discussion on page 32 from: “Additionally, investigating DNAm grouped by chromatin regions32, rather than arbitrary groupings or individual CpGm, has better functional alignment providing biologically relevant results.” To:\n“Additionally, investigating DNAm grouped by chromatin regions32, rather than individual CpGm or regions grouped by arbitrary sliding windows and significance as with comb-p22, improves identification of biologically relevant results.”\nWe have included Rebuttal Figure 5 as a new supplemental figure, Figure S6, and the comb-p results as a new supplemental table, Table S8, a small portion of which we included below given space limitations:\n“Table S8: Significant Regions Identified by Comb-p aggregation of individual CpGm EWAS. Table of results for all CpGms present with significant regions (Sidak p-value ≤ 0.05) identified through comb-p using individual CpGm EWAS results as input. To identify commonly significant regions across phenotypes, regions that overlapped across EWAS were concatenated into one region as specified by the 'Overlap_region_index' and genomic positions. 'Combp-region' columns indicate the original grouping of CpGms comb-p identified, while 'CpG' columns indicate the CpGms grouped into a region by comb-p.\nAB = Amyloid-beta; CER = Cerebellum; FA = Insoluble formic acid fraction; pTau = phospho-tau; TBS = soluble tris buffered saline fraction; tTau = total tau; TCX= Temporal Cortex; TX= membrane associated Triton-X fraction”\n\nFinally, we included the following text in the Methods on page 42:\n“Comb-p regional aggregation:\nTo compare our rCpGm method with the published regional aggregation method comb-p22, test statistics from each individual CpGm EWAS were independently used an inputs to the python based comb-p software following the default ‘pipeline’ steps45. Comb-p identifies and groups spatially correlated p-values through a sliding window approach. Input parameters included regions with ≥ 3 CpGs, a seed of < 1E-04, a sliding window size of 500 bp, and significant regions defined as those with a Sidak corrected p-value ≤ 0.05. While comb-p does not take into account the CpG regression estimates, significant regions were defined as positively or negatively associated with a phenotype by the effect directions of the majority (> 50%) of CpGs included in the defined region. When comparing across EWAS, overlapping regions were identified and concatenated using the GenomicRanges::reduce function.”\n\n1.4. In addition, the review commented on “Technical and Batch Effect Corrections” across data types stating “The study integrates data from multiple platforms, including Reduced Representation Bisulfite Sequencing (RRBS), Illumina EPIC arrays, bulk RNA-seq, and single-nucleus RNA-seq. While batch effect corrections were performed, it is unclear whether cross-platform normalization techniques (e.g., ComBat, limma-based adjustments) were applied to ensure comparability across datasets. Given the potential impact of batch effects on methylation and expression levels, a more detailed discussion on this aspect would be beneficial.”\nAs the reviewer notes, we did correct for batch by including a random effect variable in the linear models for the gene expression data and as a sensitivity analysis for the DNAm data. We chose not to apply cross-platform normalization techniques like ComBat or limma for technical reasons discussed in reviewer point 3.2. Importantly, the only cross-platform data integration we performed in our study was for data from the same individuals, namely between the DNAm and gene expression data and DNAm and genotype data. The other data types we investigated including the Illumina arrays, additional bulk RNAseq, and single nucleus RNAseq data were analyzed individually and then their association test statistics were compared across platforms. These comparisons of association test statistics consist of a meta-analysis of Braak stage which we have now included to address reviewer point 3.1 below and comparisons of effect directions and significance.\nPrior to integrating the aforementioned RRBS DNAm and gene expression data from the same individuals, we regressed out important biological and technical variables including age at death, sex, and the first three genetic principal components from the DNAm data, as well as age at death, sex, RIN, and batch from the gene expression data. We have thus performed an association between the DNAm and gene expression residuals that were obtained following regressing out these variables. Given we did not integrate sample level data across data types or platforms beyond those mentioned, we did not apply cross-platform normalization techniques.\nTo clarify this, we have edited the Discussion on page 31 from: “Our epigenome-wide association study of deep brain AD protein biochemical levels, coupled with our novel rCpGm approach, and integrated transcriptome and epigenome data from >4,000 collective brain samples demonstrate strong, replicable and biologically consistent associations between brain tau and hypermethylation near oligodendrocyte genes likely leading to their lower expression, dysmyelination and neurodegeneration in AD and other tauopathies (Figure 5C).” To:\n“Our epigenome-wide association study of deep brain AD protein biochemical levels, coupled with our novel rCpGm approach, and integrated transcriptome and epigenome data of donors from >4,000 collective brain samples demonstrate strong, replicable and biologically consistent associations between brain tau and hypermethylation near oligodendrocyte genes possibly leading to their lower expression, dysmyelination and neurodegeneration in AD and other tauopathies (Figure 5C).”\n\n1.5. Another minor concern of the reviewer was the potential to include “Longitudinal Validation of Methylation Changes” stating that “The study relies on cross-sectional data, which limits the ability to infer temporal relationships between methylation changes and disease progression. Future studies incorporating longitudinal methylation profiling would provide stronger evidence for the dynamic role of DNAm in AD pathogenesis.”\nWe fully agree with the reviewer that future studies incorporating longitudinal DNAm profiling would provide stronger evidence for the role of our top hits in the pathogenesis of AD. Since longitudinal profiling of brain tissue for neurodegenerative diseases is essentially impossible, such studies will have to rely on either cerebrospinal fluid or blood samples longitudinally collected in well characterized cohorts. The omics signatures from such studies of CSF or blood will need to be compared to omics signatures obtained in brain tissue. As such, brain-based datasets such as ours will continue to serve a critical role to anchor peripheral omics and disease endophenotypes to central (i.e. brain-based) omics and endophenotype data. To acknowledge this point, we have added the following sentence in combination with reviewer point 1.2 to the Discussion on Page 33:\n“In addition, the results of this study are based on cross-sectional data at the time of autopsy and future longitudinal blood or CSF DNAm studies and functional experimentation investigating our top hits will be important in understanding the temporal relationship between DNAm, gene expression changes and tau pathology.”\n\n1.6. Lastly, Reviewer 1 comments on “Reproducibility and Data Availability” describing that “The study claims high reproducibility across independent datasets (e.g., ROSMAP, BDR), which strengthens the validity of its findings. However, reproducibility also depends on the availability of full bioinformatics pipelines, processing scripts, and custom analytical methods. The authors should clarify whether their rCpGm approach has been made publicly available (e.g., via GitHub or Bioconductor) to allow other researchers to implement and validate the method in additional datasets.”\nWe completely agree with Reviewer 1 regarding making the approaches and pipelines utilized in this study publicly available. This is a practice we have been routinely conducting in all our publications for over a decade. We have also deposited the bioinformatic scripts and pipelines utilized for the main findings of this study on the AD knowledge portal. Access to this space can be gained after creation of a Synapse account.\nWe have edited the following text to incorporate a link to the deposited scripts on Page 48:\n“Code used in this manuscript for data analyses can be accessed via the AD Knowledge Portal: https://www.synapse.org/Synapse:syn64465911.”\n\nReviewer 2\nWe thank Reviewer 2 for co-reviewing our manuscript and fully support the initiative to facilitate training in the peer review process for early career researchers.\n\nReviewer 3\nThe third reviewer summarized our manuscript stating “This manuscript presents a new method for clustering methylation cpg data in the brain, provides robust evidence of methylation alterations across the brain in Alzheimer’s disease and in relation to the major neuropathological features of Alzheimer’s disease, replicated the most robust signals in independent datasets, and provides evidence of gene expression consequences and some proposed biological pathways through the white matter that may underlie the observed associations. The paper is impressive in its scope and rigor and provides exciting new biological insights into the pathogenesis of Alzheimer’s disease and related dementias. There are some notable limitations that should be addressed, and I lay out my concerns by category (Major or Minor) below.”\nWe thank the reviewer for their positive comments regarding the scope, rigor, and biological insights of our study. We have addressed each of the comments of our Reviewer 3 and believe that our responses and edits for the major (3.1 - 3.3) and minor concerns (3.4 - 3.8) have further strengthened our manuscript.\n\n3.1. “The replication of individual signals appears to include nominal associations in the replication dataset. I would recommend relying fully on the correction for multiple comparisons for the number of targets that were brought to replication. FDR are reported in the figure, just lean on the ones that correct for replication. Perhaps the 93 are the FDR corrected, but it was challenging to follow.”\nWe appreciate the opportunity to clarify the methods we used for rCpGm replication. We should first note that the rCpGms selected for external replications are already FDR-significant and have a biologically relevant direction of effect for two different brain endophenotypes, namely Total tau and pTau in our dataset. Having selected these rCpGms for replication with these stringent criteria, we chose to apply a more relaxed nominal p-value ≤ 0.05 threshold for external replications. The reason for this was to minimize false negative findings that may result from the limitations of our replication cohorts, BDR and ROSMAP. These limitations include the use of methylation array technologies rather than sequencing, different brain regions, dissimilar cohort compositions, and lack of tissue fraction specific tau levels, the endophenotypes of interest that emerged from our cohort. Instead of the quantitative brain biochemical measures we analyzed in our dataset, BDR and ROSMAP only had semiquantitative neuropathology measures or binary diagnostic information (AD vs. control) available, with limited power compared to quantitative endophenotypes. Due to these limitations, we anticipated that even true biological signal replications may have a lower effect size. To minimize false negative findings while maximizing true positives, we conducted multiple downstream analyses in further replication cohorts to support the functional relevance of these rCpGms and tau associations. These additional analyses include integration of methylation with matched gene expression measures and replications in an additional eight independent datasets with related tau measures. Remarkably, many of the rCpGm loci and genes we identified to be associated with tau measures and related to oligodendrocytes showed parallel associations in these additional independent datasets in directions that are biologically congruent.\nNevertheless, to address our reviewer’s comment, we have conducted additional analyses. Rebuttal Figure 6 below illustrates the number of rCpGms that would pass different multiple testing correction thresholds. If we applied an FDR p-value ≤ 0.05 threshold for the 738 rCpGms taken to replication where 600 were present in at least one of the datasets, only eight rCpGms would survive this correction. Of these 8, two (TGFA, RTKN) are part of the top 15 oligodendrocyte genes highlighted and found to have biologically consistent replications in many of the additional independent datasets analyzed in our study. The remaining 12 tau-related oligodendrocyte genes would not have passed this more stringent statistical threshold.\n\nRebuttal Figure 6: Number of rCpGms that pass various multiple testing correction thresholds in the BDR and ROSMAP datasets.\nTo further enhance the power of associations of rCpGms with Braak stage, the only tau-related endophenotype available across all cohorts, we performed a meta-analysis of the independent data from our RRBS TCX, BDR DLPFC, and ROSMAP DLPFC datasets. For this analysis, we did not include the BDR OCC dataset, given overlap of its donors with BDR DLPFC in order not to violate the independent sample assumption of meta-analysis. We chose BDR DLPFC as this is a region that is closer to TCX in terms of progression of tau pathology, and we already include DLPFC data from ROSMAP. Of the 738 rCpGms selected, 600 were present in two or more of the three datasets and could be tested. We identified 71 nominally significant rCpGm meta-analysis associations with Braak stage (fixed-effect p-value ≤ 0.05) including 14 with a negative effect direction and 57 with positive effect direction (Rebuttal Table 11).\n\nRebuttal Table 11: Summary table of significant rCpGms in the Braak stage meta-analysis of the RRBS TCX, BDR DLPFC, and ROSMAP DLPFC datasets. ‘N significant’ indicates a fixed-effect p-value ≤ 0.05.\nRemarkably, of these 71 associations, 47 are also members of the 93 rCpGms we defined as replicated suggesting that although within individual datasets these 93 rCpGms may not pass an FDR ≤ 0.05 multiple testing correction threshold, at least half (50.5%) significantly associated with Braak stage when evaluated together. This included all three tier 1 rCpGms (100%), 16/23 tier 2 (70%), and 28/67 tier 3 (42%) rCpGms. Moreover, 9/15 tau-oligodendrocyte-related rCpGms significantly associated with Braak stage in our meta-analysis including the rCpGms related to RTKN, TGFA, CLMN, MBP, CDK18, ST18, MYRF, KLK6, and LDB3 gene expression. These results suggest that although we applied a nominal p-value threshold to identify the 93 replicated rCpGms, our meta-analysis results provide support that at least half of these rCpGms associate more broadly with the tau related neuropathology measure of Braak stage. If we applied an FDR p-value ≤ 0.05 threshold, we would exclude many of these associations suggesting that given the limitations of this replication analysis outlined above and our downstream efforts to further replicate top signals in multiple datasets, an FDR p-value ≤ 0.05 threshold is likely too conservative.\nWe have incorporated the Braak stage meta-analysis results into Figure 4 in combination with changes suggested in reviewer point 4.8.8.\nWe edited the Figure 4 legend from: “The y-axis depicts the 93 rCpGms that are FDR significant (FDR pval < 0.05) in the Mayo Clinic RRBS-TCX total tau and pTau EWAS and have opposite directions of associations with these tau endophenotypes. These rCpGms are nominally significant (p-value < 0.05) and have congruent direction of associations with Braak stage in ROSMAP and/or BDR. Tier 1 rCpGms (N= 3) are associated with Braak stage in all three replication datasets (BDR-DLPFC, BDR-OCC, and ROSMAP-DLPFC); Tier 2 (N=23) are associated with Braak stage in two of the replication datasets; Tier 3 (N=67) are associated with Braak stage in one of the replication datasets. Also included are additional associations of the rCpGms with Thal phase and AD diagnosis (Dx) with significant associations defined as those with a p-value < 0.05. Non-colored associations in ROSMAP indicate that this region was not measured on the array. Each rCpGm is additionally annotated with the top FDR significant gene expression association (GE~rCpGm) colored by direction of effect, their nearest known AD risk loci (+/- 1Mb, Bellenguez 2022), and by their chromatin state annotation from the Roadmap Epigenomics 15 chromatin state model. See supplemental data Table S12 to connect rCpGm Index numbers with rCpGm positional data. Red indicates a positive direction of effect, blue indicates a negative direction of effect, and the transparency indicates if the association passed nominal significance threshold.\nBDR = Brains for dementia research; CpGm = DNA methylation at CpG sites; DLPFC = Dorsolateral prefrontal cortex; Dx = Diagnosis; EWAS = Epigenome-wide association study; FDR = false discovery rate; GE = Gene expression; N = Number; OCC = Occipital cortex; pTau = phospho-tau; rCpGm = Regional DNA methylation grouped by chromatin state; RRBS = Reduced representation bisulfite sequencing; TBS = soluble tris buffered saline fraction; TCX = Temporal cortex; TX= membrane tissue fraction” To:\nFigure 4: Replication of FDR significant brain tau-associated rCpGms in independent datasets. The y-axis depicts the 93 rCpGms that are FDR significant (FDR pval < 0.05) in the Mayo Clinic RRBS-TCX total tau and pTau EWAS and have opposite directions of associations with these tau endophenotypes. These rCpGms are nominally significant (p-value < 0.05) and have congruent direction of associations with Braak stage in ROSMAP and/or BDR. Tier 1 rCpGms (N= 3) are associated with Braak stage in all three replication datasets (BDR-DLPFC, BDR-OCC, and ROSMAP-DLPFC); Tier 2 (N=23) are associated with Braak stage in two of the replication datasets; Tier 3 (N=67) are associated with Braak stage in one of the replication datasets. In addition, a meta-analysis of Braak stage was conducted including the RRBS-TCX, BDR-DLPFC, and ROSMAP-DLPFC datasets. Meta-analysis fixed-effect estimates and p-values are reported. Also included are additional associations of the rCpGms with Thal phase and AD diagnosis (Dx) with significant associations defined as those with a nominal p-value < 0.05. Non-colored associations in ROSMAP indicate that this region was not measured on the array. Each rCpGm is additionally annotated with the top FDR significant gene expression association (GE~rCpGm) colored by direction of effect and by their chromatin state annotation from the Roadmap Epigenomics 15 chromatin state model. See supplemental data Table S12 to connect rCpGm Index numbers with rCpGm positional data. Red indicates a positive direction of effect, blue indicates a negative direction of effect, and the transparency indicates if the association passed nominal significance threshold.\nBDR = Brains for dementia research; CpGm = DNA methylation at CpG sites; DLPFC = Dorsolateral prefrontal cortex; Dx = Diagnosis; EWAS = Epigenome-wide association study; FDR = false discovery rate; GE = Gene expression; N = Number; OCC = Occipital cortex; pTau = phospho-tau; rCpGm = Regional DNA methylation grouped by chromatin state; RRBS = Reduced representation bisulfite sequencing; TBS = soluble tris buffered saline fraction; TCX = Temporal cortex; TX= membrane tissue fraction”\nWe have also included the following text in the Results on page 18:\n“We also performed a meta-analysis of Braak stage including the RRBS-TCX, BDR-DLPFC, and ROSMAP-DLPFC datasets and found that 47 of the 93 replicated rCpGms were significant (fixed p-value ≤ 0.05) including all of tier 1 rCpGms, 16/23 tier 2 rCpGms, and 28/67 tier 3 rCpGms (Figure 4, Table S11).”\nWe edited the following text in the Discussion on page 28 from: “Remarkably, we found that of the 600 tau-related rCpGms in our cohort, which were also available from these replication datasets, 93 showed significant and biologically congruent associations in BDR and/or ROSMAP” To:\n“Remarkably, we found that of the 600 tau-related rCpGms in our cohort which were also available from these replication datasets, 93 showed significant and biologically congruent associations in BDR and/or ROSMAP of which half were also significant in a Braak stage meta-analysis.”\nWe added the following Methods text on page 47:\n“Meta-Analysis\nMeta-analysis of Braak stage was performed in R with the meta::metagen function including rCpGm EWAS estimate and standard error statistics from the RRBS-TCX, BDR- DLPFC, and ROSMAP-DLPFC datasets. Test statistics from both fixed and random effect models are reported (Table S11).”\nWe also included the meta-analysis results as a new supplemental table, Table S11, a portion of which we have included below:\nTable S11: Braak Stage Meta-Analysis Results for rCpGm Replication.\n\n3.2. “It is unclear from the methods why only the top 50% were used for PCs during quality control and to what degree lane, participant, relatedness, or technical artifacts were accounted for in QC procedures. It seems like there is a substantial mean level difference across brain regions, I’m curious if some type of regression based correction like combatmet was needed or applied?”\nWe choose to include the top 50% most variable CpGs as this aligns with the default methods used is the R package MethylKit33.\nWhile we agree there is a significant difference in global mean DNAm levels between the CER and TCX tissue types, the magnitude of this difference was 0.0217 in an unpaired t-test across all samples or 0.0204 in a paired t-test with only the 190 paired samples. These magnitudes indicate the TCX has 2.04%-2.17% more DNAm on average than CER tissue in the genomic regions we captured through the reduced representation bisulfite sequencing (RRBS) method we used. While RRBS captures more CpG sites than current array technologies and enriches for regions with CpG islands, the captured regions are still only a fraction of the entire genome. The average percent difference of DNAm between TCX and CER indicates that of the RRBS captured CpGs, there is a small fraction of difference between methylation levels of these two brain regions.\nTo determine whether these DNAm differences were due to technical artifacts, we further investigated the global average DNAm level differences between TCX and CER tissue across each batch of sequencing. We found that regardless of sequencing batch flowcell, we consistently see higher levels of average DNAm levels in TCX tissue compared to CER tissue (Rebuttal Figure 7), which suggests that the differences in global mean averages are not a result of sequencing batch differences.\n\nRebuttal Figure 7: Average Global CpGm levels across tissue type and sequencing batch flowcell.\nWe did not account for participant in our QC procedures because we performed QC separately for the TCX and CER samples such that there were no samples from the same individual in the same QC workflow. We also did not account for participant in downstream analyses because we performed our association tests for the TCX and CER data separately. As such, we did not have any donors with DNAm measures from multiple samples in the same analysis.\nTo address relatedness, we leveraged relatedness estimates from our previous genetic studies34,35 investigating genotype data from the same individuals used in this study. Briefly, we evaluated relatedness via the identity by descent (IBD) function in PLINK (-genome) using all independent autosomal genetic variants pruned on linkage disequilibrium (r2 > 0.8). Related pairs were defined as those with a PI_HAT score > 0.125 (3rd degree relatives or more) and the sample with the highest call rate was retained. This procedure was iteratively performed until no related individuals were detected. In our current study, the third step of our QC workflow removed individuals that previously failed the genome-wide association study QC metrics including the relatedness metric. This identified only two related pairs and therefore led to the exclusion of two samples from those pairs in this study.\nWhile we considered applying regression-based correction for additional technical variables using techniques like ComBat-met, ComBat, or limma during our QC procedure, these methods require no missing values to calculate corrected values. As a result, we would need to either remove these missing values or impute them to apply these methods. While others have imputed missing data by using the average of all other values at that position or using a reference database, we hesitate to introduce unmeasured values into our dataset and were not able to find a reference suitable for our RRBS dataset of AD samples from brain tissue. In addition, methods like ComBat and limma were built for use with gene expression data. Although these methods have previously been applied to other data types like DNAm data, typically β-values are transformed into M-values before applying this method. We, on the other hand, are interested in using rCpGm values which are β-values averaged across a region rather than transformed M-values. Moreover, the authors of the limma package36 advise that the removeBatchEffect function not be used for linear modeling and that it would be better to include batch factors as variables in the linear model which we did as a sensitivity analysis for our top hits.\nLastly, for sample preparation and sequencing, we randomized all samples for brain tissue, sex, age at death, APOE-ɛ4 dose and CAA score to help prevent introduction of batch effects. Nonetheless, to further address the reviewer’s comment, we investigated the correlation of RRBS TCX sequencing batch with the total tau and pTau measures through ANOVA, and found no significant association between RRBS sequencing batch and either total tau (p-value = 0.719) or pTau (p-value = 0.86) measures (Rebuttal Figure 8). In summary, there is no effect of sequencing batch on rCpGm values nor tau endophenotypes.\n\nRebuttal Figure 8: Distribution of Total tau and pTau measures across RRBS sequencing batches.\n\n3.3. “It seems that the cis-rCpG to gene analysis was performed across all samples, which would confound any association in which disease was driving both gene expression and methylation. I would recommend focusing the cis-QTL analysis on the normal cognition/path negative participants to remove this potential confound.”\nWe appreciate the reviewer’s comment regarding the effects of disease in the relationship between cis-rCpGms and genes. We would like to re-iterate that our cis-rCpGm to gene expression associations are conducted in the Mayo Clinic RRBS dataset, in which all samples are neuropathologically diagnosed Alzheimer’s Disease cases (please see Methods, page 35). All donors have a Braak stage ≥ 4 and Thal phase ≥ 2. Hence, we would not expect there to be any confounding based on diagnosis.\n\n3.4. “I wonder if it would be worth using Hi-C and plaq-seq data from brain to try and infer some potential gene targets beyond the cis-QTL analysis provided?”\nWe think this is an interesting suggestion and has the possibility to identify potentially new connections between rCpGms and genes, however, we consider this beyond the scope of our current study. We believe one of the strengths of our study is that the DNAm and gene expression data originate from the same donors improving our ability to identify regulatory networks. In doing this, we identified potentially important networks of regulation between tau, DNAm and nearby gene expression which we then focus on and characterize further. Nonetheless, future studies that generate similar types of data or additional data types from the samples we analyzed in this current study should consider generating chromatin conformation data types to better understand all possible chromatin regulatory networks. To express this, we have included the following sentence in the Discussion on page 33-34:\n“Moreover, our investigation of DNAm regulating gene expression was restricted to nearby genes, so it is possible we may have missed important longer range regulatory networks which future studies may be able to identify by incorporation of chromatin conformation data.”\n\n3.5. “The rationale for including CDK18 in the final list is weak and should be strengthened.”\nTo strengthen the rationale for including CDK18 in the final list, we have changed the following sentence in the Results on page 20 from: “Although part of the 1000, and not the 100 top oligodendrocyte marker genes2, CKD18 was previously implicated with tau37 and AD DNAm changes38, thereby leading us to include this gene in additional analyses.” To:\n“Although part of the 1000, and not the 100 top oligodendrocyte marker genes2, CKD18 was previously implicated in tau phosphorylation including at the T231 site and tau aggregation37,39, found contained within NFT37, differentially expressed in TCX tissue of individuals with AD compared to controls37, associated with AD DNAm changes38, and helps prevent DNA damage and promotes genomic stability40, thereby leading us to include this gene in additional analyses.”\n\n3.6. “The analysis of the oligodendrocyte genes seems circular. The CpGs were identified in relation to tau, then the authors focus on genes that correlate with those same CpGs in the same dataset, and then argue that correlation of those same targets with tau demonstrates the relevance oligodendrocyte expression to tau.”\nWe appreciate the opportunity to clarify this seemingly circular argument. The reviewer is correct that we are investigating a stepwise pattern of associations, however, we believe that doing so strengthens our conclusions. This is because across these separate datatypes we find association patterns that not only make statistical sense but are also biologically congruent. We identified a group of rCpGms that not only associate with total tau and pTau in opposite effect directions but also associate negatively with expression of nearby oligodendrocyte genes. If these associations are biologically relevant rather than, for example, a potential technical artifact in the DNAm data, then we would expect the association of gene expression (a different data modality) to the tau phenotypes to also be significant and have a direction of effect consistent to that of the rCpGm to tau phenotype associations. In testing gene expression to tau phenotype associations, we indeed largely identified this biologically congruent pattern further supporting our conclusion that these findings reflect true biology. These related but most importantly consistent associations across three separate data types suggest that the oligodendrocyte related signal we identified and importantly further replicated across multiple tauopathies is biologically important. While we acknowledge these are only associations, we believe the consistent connection across levels of DNAm, gene expression, and tau within our dataset as well as across multiple independent datasets provide support to our conclusions.\n\n3.7. “The final step of replication of gene expression associations with tau pathology should also include a correction for multiple comparisons.”\nWe have now included FDR corrected p-values for the gene expression associations with tau pathology in Table S18, shown below:\nTable S18: Associations of the top oligodendrocyte marker gene brain levels with Braak stage in three bulk brain RNAseq AMP-AD datasets.\nWe have also edited the Results to report these FDR corrected p-values on page 22 from: “Expression of LDB3 was associated with Braak stage in the Mayo TCX (beta = -0.19; p-value = 8.47E-3) and ROSMAP (beta = -0.058; p-value = 7.22E-3) datasets with a biologically congruent negative direction of effect (Table S15).” To:\n“Expression of LDB3 was associated with Braak stage in the Mayo TCX (beta = -0.19; FDR p-value = 2.54E-2) and ROSMAP (beta = -0.058; FDR p-value = 2.52E-2) datasets with a biologically congruent negative direction of effect (Table S18).”\n\n3.8. “The causal order that is proposed could be better evaluated with either an MR approach (perhaps challenging if the methylation changes aren’t under genetic control) or a mediation approach. It seems equally likely that methylation is driving oligo expression resulting in a feed forward exacerbation of tau deposition/propagation, as it is that the tangles drive methylation difference resulting in white matter damage. Without evidence of white matter damage in these participants, the causal argument is challenging to follow or justify in the current manuscript. I’m not sure what my recommendation is, I just wasn’t very convinced by the argument around order of events.”\nWe agree that given the association results, both scenarios are likely, i.e. changes in rCpGm drive gene expression changes resulting in changes in tau levels vs. changes in tau levels drive rCpGm changes and lead to gene expression perturbations. To examine potential evidence for causation one way or the other, we first employed Mendelian Randomization (MR) analysis. As described above in reviewer point 1.2 where we test for evidence that changes in rCpGm levels cause changes in gene expression through MR, only 4 rCpGms had significant mQTL associations (FDR p-value ≤ 0.05) with nearby SNPs and could be tested with MR. Using the 8 independent genetic variants detailed in Rebuttal Table 3, we tested if rCpGms cause changes in total tau or pTau using the methodology described in point 1.2.\nWe found a significant MR result suggesting that changes in the MBP related rCpGm leads to changes in levels of total tau using the weighted median MR method (p-value = 0.04) and a suggestive association using the inverse weighted mean MR method (p-value = 0.056) (Rebuttal Table 12). We performed multiple sensitivity tests and found no evidence of horizontal pleiotropy (MR Egger intercept = 0.05 & p-value = 0.744) nor heterogeneity (Cochran Q statistic p-value > 0.5). Causal directionality via the Steiger test41 suggests that this is likely in the correct direction (Steiger p-value = 4.14E-08). Investigating the contributions of each SNP to the overall inverse variance weighted effect size also suggests that a single SNP is not driving this outcome (Rebuttal Figure 9).\n\nRebuttal Table 12: Mendelian Randomization results for rCpGms causing changes in tau levels.\nRebuttal Figure 9: Forest Plot of SNPs contributing to the Mendelian Randomization effect size for the MBP rCpGm on Total tau .\nThese results suggest that genetically predicted rCpGm levels at the MBP locus may cause a reduction of total tau levels. However, although the 4 SNPs that were included in the IV had prediction ability (F-statistic > 10) and they explained ~10% of the rCpGm variance (Rebuttal Table 3), we only identified significance using the weighted median method and only suggestive significance using the inverse variance weighted method. Collectively, these results imply uncertainty in the casual signal identified and do not remain significant when accounting for multiple test corrections (0.05/2 = 0.025).\nWe were also not able to identify definitive evidence employing casual mediation as described in reviewer point 1.2 for the rCpGm, gene expression and tau variables. Taken together, despite conducting all suggested analyses, we are not able to confidently determine a causal direction unequivocally. We have therefore edited the text of the manuscript in multiple sections to allow for the multiple possibilities in causal direction and the need for functional experiments. These edits are described above in reviewer point 1.2.\n\nReviewer #4\nIn response to the question “What are the noteworthy results?”, reviewer 4 states “The authors, using temporal cortex and cerebellum samples from over 400 donors, provide a comprehensive exploration of RRBS analysis of individual CpG sites and chromatin-based regional analysis for combining CpG sites. These two methods are then associated with several AD-related endophenotypes, including amyloid and tau. Both these methods of analysis identify several genomic locations with altered methylation, with the majority of FDR significant results being in tau measures. These results provide overlap with previous results on brain pathological array-based methods. Using ROSMAP and BDR, which have array-based DNA methylation measures pre-frontal cortex and occipital cortex, the authors show consistency in identified sites and chromatin regional analyses, although these associations use Braak, Thal and add cerebral amyloid angiopathy score. Investigating significant findings, the authors show oligodendrocyte markers and myelination-related genes using gene-set enrichment. They also show a correlation between regional and single CpG methylation and gene expression and identify mQTLs from tau-associated sites.”\nWe thank reviewer 4 for their in-depth summary of our manuscript.\n\n4.1. In response to the question “Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not original, please provide relevant references.” The reviewer comments “The work is of significance to the field of AD research and the epigenetics of AD. It shows good overlap with previous analyses into similar methods and tissues, but also provides an alternate regional analysis for DNA methylation sites with potentially more biological relevance. The work is original in that it uses a new regional DNA methylation method and applies it to previously published data, showing its applicability and validity. This method also allows for the identification of new regions and genes of interest for integrative genomic analyses.”\nWe greatly appreciate the positive feedback regarding the significance of our study to the AD and epigenetics of AD fields as well as the novelty of our regional DNA methylation method.\n\n4.2. In response to the question “Does the work support the conclusions and claims, or is additional evidence needed?” the reviewer affirms “I believe the work supports the conclusions given.”.\nWe thank the Reviewer for their assessment.\n\n4.3. In response to the question “Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit publication or require revision?” the reviewer notes that “The data analyses seem appropriate in method, and appropriate covariates, multiple testing adjustments, and post-hoc analyses have been used. My only revision is the ability to find the full script. The authors state that scripts are available via the AD Knowledge Portal, however, the portal is sometimes difficult to navigate, and I was unable to find them. This may be due to my lack of familiarity with the portal, but a direct link or instruction is needed for readers who face a similar problem.”\nAs described in point 1.6 above, we have now made the scripts available via the AD knowledge portal at this link: https://www.synapse.org/Synapse:syn64465911\n\n4.4. Reviewer 4 answers the questions “Is the methodology sound? Does the work meet the expected standards in your field?” stating “Samples included are well described and appropriate to the analyses. Methodology is sound and meets the standards in the field. Work has been subjected to power calculations for identifying of appropriate effect sizes and using the sample sizes given. Appropriate use of Spearman’s and Pearson’s correlation for normal and non-normal/categorial measures is used. Multiple testing using FDR is used. CpG ~ AgeAtDeath + Sex + genetic PC 1-3 is used for EWAS analyses, which seems appropriate. Similarly, the methodology for gene expression association and mQTLs seems appropriate. Post-hoc analyses are carried out as per standards.”\nWe thank the reviewer for their careful review of the methodologies we undertook in this study, finding our application of these to be appropriate and our methodology to be sound.\n\n4.5. In response to the question “Is there enough detail provided in the methods for the work to be reproduced?” the reviewer states that “Once access to the code is made easier and can be checked, yes.”\nWe agree and have made these available as described in more detail in points 1.6 and 4.3 above.\n\n4.6. A major concern reviewer had was that “The website doesn’t allow access as it requires a login. Please make it accessible for the exploration of data.”\nWhile the manuscript is under review, we have required a login to access the data available through the shiny app. Once accepted for publication, we will remove this requirement making the app freely accessible to the scientific community.\nWe had already provided the login details in the cover letter to enable access for the reviewers and will provide them again here:\nUsername: Reviewer\nPass: passForReviewer!\n\n4.7. In addition, the reviewer points to Line 890 requesting to “Please provide a direct link or details on where to find the code on the knowledge portal.”\nWe have now included a direct link to the code outlined above in points 1.6 and 4.3.\n\n4.8. The reviewer also had multiple minor concerns regarding “writing or consistency errors which may be present in other points.”\nWe thank the reviewer for noting these writing and consistency points.\n\n4.8.1. For brevity, we have made each of the following simple text changes listed below:\n• “Line 14 (and others) – 6 to six and 2 to two” – now lines 14-15\n• “Figure 1 – Define TCX in the acronym list in the legend” – completed\n• “Line 75 – remove “CpG sites (CpGm)” as already defined in intro” – completed\n• “Line 596 – four to 4, two to 2” – now line 656\n• “Table S1 (and others) – Define acronyms (TCX, CER) in legend” – completed for relevant tables\n• “Line 619 – Genome wide to Genome-wide for consistency” – now line 679\n• “Line 622-623 – QC is used for quality control and quality controlled. Suggest changing line 622-3 to be “and underwent rigorous quality control (QC)”.” – now lines 682-683\n• “Line 624 – R2 to R^2 (meaning superscript, limit of the text box)” – now line 684\n• “Line 664 – quality control to QC” – now line 748\n• “Line 829 – non-Hispanic” – now line 929\n• “Line 851-2 – “defined as individuals with a Braak stage ≥ 4 and a Thal phase ≥ 2. Control samples were defined as individuals with a Braak stage ≤ 4 and a Thal phase ≤ 2” from reading table S3, this should be Braak stage < 4 and Thal phase < 2.” – now lines 953-956\n• “Figure S2 and others - Acronyms for CER, TCX, CpGm and rCpGm need defining” – completed for relevant figures\nThe reviewer also suggested a change in “Figure S1 - stdev to SD for consistency with Table 1”. While Figure S1 does not have a ‘stdev’ acronym used, there is the acronym ‘std’ which we have defined as ‘standardized’ meaning these values were z-score transformed to have a mean of zero and standard deviation of 1.\nThe reviewer also notes additional non-simple text changes including:\n\n4.8.2. “Line 223-229 and Fig S6 – Although I agree that the correlation is there, the distribution of beta values mean the data is often skewed due to the ends of the association, especially in very large datasets. I would suggest using a random sampling or subsetting of data (especially in the RRBS vs. array correlations).”\nWe agree that the skewed distribution of beta values may impact correlation tests which is the reason we choose to use the non-parametric Spearman’s rank correlation test to account for this. We have performed additional Spearman’s rank correlation tests for each of the comparisons using a random sampling of the same 5,000 rCpGms across datasets. We found consistent and even stronger correlation coefficients when using a random sampling of 5,000 rCpGms compared to using all common rCpGms as outlined in Rebuttal Table 13 as well as Rebuttal Figure 10. These are consistent with our previous Figure S6 now Figure S7.\n\nRebuttal Table 13: Spearman’s Rank test results for all common rCpGms and a random sampling of 5,000 common rCpGms across the DNAm datasets investigated.\nRebuttal Figure S10: Correlation of rCpGm levels averaged within a study dataset across samples to overlapping rCpGm sites from other study datasets.\nWe have now updated the previous Figure S6 to include both the original correlations as panel A and Rebuttal Figure 10 as panel B in the new Figure S7.\n\n4.8.3. “Line 281-293 – Please include a table for GSEA results in supplementary and refer in the text”\nWe have now included a new supplemental table (Table S15) that includes the gene set enrichment analysis results for the top 10 enriched pathways originally referenced. A condensed version of this new Table S15 is below:\n“Table S15: Top 10 enriched GO:BP Pathway results for genes significantly associated with 93 rCpGms of interest.\nWe have also edited the Results text on page 19-20 to reference this new table from: “The top three nominally enriched gene ontology terms for these genes relate to positive regulation of myelination (GO:0031643, pval = 3.4E-5), gliogenesis (GO:0042063, pval = 1.4E-4), and microtubule nucleation (GO:0007020, pval = 1.5E-4).” To:\n“The top three nominally enriched gene ontology terms for these genes relate to positive regulation of myelination (GO:0031643, pval = 3.4E-5), gliogenesis (GO:0042063, pval = 1.4E-4), and microtubule nucleation (GO:0007020, pval = 1.5E-4) (Table S15).”\n\n4.8.4. “Line 622-3 – Details or reference on the QC pipeline should be provided”\nWe have now included a direct link to the QC pipeline scripts detailed above in points 1.6 and 4.3.\n\n4.8.5. “Figure S3 – For a more direct comparison between datasets, it would be good to include the current study data with BDR and ROSMAP to more easily compare the distributions. RRBS is going to cover more sites, but it would be interesting to see the array selection of sites and how they compare. As stated, arrays have a narrower range of values, so it would be interesting to see this clearly. ROSMAP shows higher average methylation in rCpG analyses (FigS3). Is this for a reason, or a concern with comparison, or proposed to be an artifact if the differences in the arrays (in probe number) and site selection?”\nTo better see the relationship between the RRBS and array studies, we have now included the RRBS, ROSMAP, BDR, and Weber 2018 datasets in the same figure as shown below in Rebuttal Figure 11 and included this in the manuscript as the new Figure S3.\n\nRebuttal Figure 11: Characterization of the rCpGms across annotated predicted chromatin states for all datasets.\nWe think the reviewer’s suggestion that site selection could be the reason why we see differences in average chromatin states is interesting. We investigated this possibility by analyzing all sites common across the six datasets used in our manuscript (Mayo Clinic RRBS CER & TCX, BDR DLPFC & OCC EPICv1, and ROSMAP DLPFC & Weber PFC 450k). There were a total of 70,008 rCpGms that originated from the same genomic chromatin state windows across all six datasets which we define as common rCpGms (Rebuttal Figure 11C). Using only these 70,008 rCpGms, we calculated the average rCpGm values for each chromatin state as we had done previously finding largely similar trends for each chromatin state and dataset (Rebuttal Figure 11D) as we do when including all rCpGms (Rebuttal Figure 11B). We do, however, observe that the ROSMAP and Weber PSP datasets, both of which were measured on the 450K array, are highly similar suggesting that the higher rCpGm averages we observe in ROSMAP are a technical result of the array rather than dataset specific. Interestingly, comparing across datasets, it seems like the arrays do not measure low levels of methylation well in either the 450K or EPIC versions of arrays compared to the RRBS sequencing data (e.g. see Active TSS, Bivalent enhancer etc.), yielding higher methylation averages. Additionally, the 450K may be able to measure hypermethylated regions better than the EPICv1 (especially see Heterochromatin, Quiescent/low regions). Moreover, we observe an overall upward shift of 450K averages across all of the chromatin states compared to EPIC arrays and across most chromatin states in comparison to RRBS.\nWhile these observations are interesting and crucially important to consider when comparing absolute methylation values across datasets generated with different methodologies, this is not what we do in our study. Rather, we compare association results from tests performed within each dataset to identify rCpGms of interest and rely on direction of association rather than magnitude of association such that these differences in absolute measures should not impact our conclusions.\nWe have now edited Results text on page 12 from: “These independent datasets are strikingly similar in the levels and patterns of their average rCpGm for each chromatin state annotation across these datasets (Figure S3) and in comparison to our dataset (Figure 2), despite differences in donors, brain regions, collection sites and methylation measurement methods highlighting the replicability of our novel analytic approach. Notably, the array-based rCpGms have a narrower range of values (Figure S3) in comparison to our RRBS-based (Figure 2) measurements, likely due to lower precision and dynamic range for the former methodology42,43.” To:\n“These independent datasets are strikingly similar in the levels and patterns of their average rCpGm for each chromatin state annotation across these datasets (Figure S3) and in comparison to our dataset (Figure 2), despite differences in donors, brain regions, collection sites and methylation measurement methods highlighting the replicability of our novel analytic approach. These average rCpGm patterns are also remarkably consistent even when examining only those rCpGms common to all datasets (Figure S3). Notably, the array-based rCpGms have a narrower range of values (Figure S3) in comparison to our RRBS-based measurements, particularly on the hypomethylated end which may potentially be due to lower precision and dynamic range for the former methodology42,43. Even within array-based datasets, we observe differences across chromatin states where 450K arrays have collectively higher rCpGm averages than the EPICv1 array. These differences indicate potential methodology-specific impacts on absolute methylation value measurements, which are unlikely to influence our results since we analyze methylation levels within each dataset. We do not recommend comparisons of absolute rCpGm values across datasets measured with independent technologies.”\n\n4.8.6. “Table S4 and S5 - Add gene annotations to the table for easier interpretation and alignment with the text (line 156-163).”\nWe have now added gene annotations to Table S4 and Table S5 and included a small portion of Table S5 below:\nTable S5: Nominally significant (pvalue <= 0.01) association results of CER individual CpGms for each AD neuropathology endophenotype.\nWe have also edited Methods text on page 44 to reflect this annotation from: “CpGm sites were annotated (Island, shore, shelf, inter-genic) with the R package annotatr() v1.21.044 querying AnnotationHub_3.2.2 in R (v4.1.2).” To:\n“CpGm sites were annotated according to their regions (Island, shore, shelf, inter-genic) and nearby genes (+/- 50kb) using the R package annotatr() v1.21.044 querying AnnotationHub_3.2.2 in R (v4.1.2) or with built-in hg38 genes, respectively.”\n\n4.8.7. “Figure S5 - Several plots appear to only have one or very few data points, so correlation is not assessable. These should be removed and/or explained.”\nThank you for raising this concern for Figure S5, now Figure S6. To assess correlation between primary and sensitivity model test statistics, we combined all phenotypes per sensitivity test and then performed the Spearman’s rank correlation test.\nTo improve clarity of this figure, we have removed the blue [1,1] lines and added correlation statistics to the top of each panel as shown in Rebuttal Figure 12. Panel D was added in response to Reviewer comment 1.1. We have also edited the legend text from:\n“Figure S5: Comparison of Primary EWAS model results with sensitivity EWAS model results for FDR significant TCX rCpGms.\nComparison of test statistics (estimates and p-values) for all FDR significant rCpGms from the primary rCpGm EWAS models and the sensitivity EWAS models. Sensitivity models included additional covariates in the linear regression models including A) sequencing batch as a random effect, B) Braak stage and Thal phase, C) Neuronal proportion estimated from bulk RNAseq data using DSA, and D) a spearman rank test. Primary model statistics are the x-axes, sensitivity model statistics are the y- axes. Blue diagonal line is the [1,1] line. Green line is the –log10(pval = 0.05) value.\nFA= insoluble tissue fraction, ln = natural log transformed, sqrt = square root transformed, std = standardized, rCpGm = Regional DNA methylation grouped by chromatin state, TBS = soluble tissue fraction, TCX = Temporal Cortex, TX= membrane tissue fraction” to that of Rebuttal Figure 12.\n\nRebuttal Figure 12: Comparison of Primary EWAS model results with sensitivity EWAS model results for FDR significant TCX rCpGms.\n\n4.8.8. “Figure 4 – Very busy and a lot of information, especially with the multiple outcomes being associated and is difficult to parse all the information, especially since some measures don’t overlap. This is especially difficult since association in the associations in some of the cohorts defining tiers is from one measure.”\nWe appreciate the reviewers’ comment and have revised Figure 4 to improve clarity, please see reviewer point 3.1 for the revised figure and figure legend. For simplification, we have now removed the nearby AD risk genome-wide significant loci annotations from Figure 4 and added them to Table S12. We edited the Results text on page 18 from: “When we annotated these 93 rCpGms for proximal (+/- 1Mb) variants identified in AD genome-wide studies14, we found that a number of them reside near known AD risk loci including JAZF1, SIGLEC11, SCIMP, KLF16, BIN1, SPI1, ABI3, ANKH, WDR12, and SIGLEC1114 (Figure 4).” To:\n“When we annotated these 93 rCpGms for proximal (+/- 1Mb) variants identified in AD genome-wide studies14, we found that a number of them reside near known AD risk loci including JAZF1, SIGLEC11, SCIMP, KLF16, BIN1, SPI1, ABI3, ANKH, WDR12, and SIGLEC1114 (Table S12).”\n\n4.8.9. “Figure 5 – Text at A4 size is too small to see clearly and image looks blurry until very zoomed in. This figure needs to have text increased or possibly removing some information to make it viewable for a main figure. Advice may be to keep part C or A and B and move the others to supplementary if figure limit has been reached.”\nWe have now increased the resolution and font size of all text within Figure 5.\n\nIntegrative Epigenomic Landscape of Alzheimer’s Disease Brains Reveals Oligodendrocyte Molecular Perturbations Associated with Tau.\nResponses to Reviewers’ Comments\nNCOMMS-25-13777A\nWe thank the Editors and the Reviewers for their careful review of our point-by-point response and revised manuscript. We are grateful for all the positive comments from the Reviewers. Below, we provide a point-by-point response to all comments. We believe that these additional modifications to address the second round of reviewer comments have further enhanced the impact of our study.\n\nSpecific points provided by the reviewers are highlighted in italic font and any text edits to the manuscript are provided in blue Times New Roman font. In our submission, we provide two versions of our Main text file; one in which all changes are tracked and another one in which all changes are accepted. Page numbers align with the tracked changes manuscript version.\n\nReviewer #2\nReviewer 2 comments that “The authors responses to my prior comments were comprehensive and the additional analyses have clearly strengthened their claims.”\nWe thank the reviewer for their evaluation of our revisions and encouraging comments on the updated manuscript.\nIn addition, they state that “A few additional points have been provided below for potential improvement” and that “[w]ith these minor clarifications this manuscript would be greatly strengthened for final acceptance”.\nWe greatly appreciate this Reviewer’s positive comments on the final acceptance of our manuscript.\n\n2.1. First, “Cell-Type Specificity of Methylation Signals: Adjustments implemented using BRETIGEA/DSA demonstrated consistent results as well as further replication using Shireby et al was explored. It may be helpful to explicitly state limitations of relying on marker-based deconvolution in regions with complex pathology due to AD.”\nWe have included this additional limitation by editing the following section in the Discussion on page 32 from “Although the associations for the top hits we identified remained even after adjusting for neuronal or oligodendrocyte proportions and showed consistent and larger magnitude associations in glial enriched nuclei, future single cell type specific DNAm studies, particularly those from glial cell types32, will be important in parsing out potential cell type specific rCpGm effects or the role of cellular heterogeneity in these associations.”\nTo:\n“Although the associations for the top hits we identified remained even after adjusting for neuronal or oligodendrocyte proportions and showed consistent and larger magnitude associations in glial enriched nuclei, the marker-based deconvolution methods we used to estimate cell type proportions have limitations, particularly in tissue impacted by pathology. As a result, future single cell type specific DNAm studies, particularly those from glial cell types32, will be important in parsing out potential cell type specific rCpGm effects or the role of cellular heterogeneity in these associations.”\n\n2.2. Second, “Causal Inference Between Methylation and Gene Expression: Both Mendelian Randomization and Mediation analysis were performed and although the authors acknowledge the power limitations of MR, recommending future study designs like larger mQTL reference panels or integration with existing consortia data would be beneficial.”\nWe have included the following sentence in the Discussion on page 33:\n“Similarly, future large scale mQTL studies would be beneficial to further understand underlying mechanisms of genetic regulation as well as aid in methods of causal inference.”\n\nReviewer #4\nReviewer #4 remarks “After reviewing the resubmission, I am happy that the authors have addressed all concerns. The results add to the field by providing a continuation of the epigenetics of AD, while expanding on the work investigating transcriptome-methylome relationships. I believe the work has significance to the field, and the work supports the conclusions provided. Having now accessed the code and supporting website, I am satisfied with the data analysis and methodology taken, and believe the interpretation and conclusions are sound. The methodology is now detailed enough to reproduce the results.”\nIn addition, they remark regarding code availability that the “README is present explaining all code parts. Code is present for all analyses and well annotated. The code was installable but I didn't fully run it, but was reviewed in detail to check steps and annotations.”\nWe thank reviewer 4 for their affirmative comments after reviewing our revised manuscript and detailed examination of the code and supporting website.\n\nReviewer #5\nReviewer 5 “[o]n behalf of reviewer #3” states that “I am confident the authors have appropriately addressed all concerns raised. This is a significant piece of research both in size and relevance for the field. In particular the method used to investigate regional CpGm (rCpGm) will be of considerable interest for the field.”\nWe thank and appreciate reviewer 5’s positive comments on behalf of reviewer 3."}]

Decision letter