Sex-specific DNA methylation differences in Alzheimer’s disease pathology

Sex is an important factor that contributes to the clinical and biological heterogeneities in Alzheimer’s disease (AD), but the regulatory mechanisms underlying sex disparity in AD are still not well understood. DNA methylation is an important epigenetic modification that regulates gene transcription and is known to be involved in AD. We performed the first large-scale sex-specific meta-analysis of DNA methylation differences in AD neuropathology, by re-analyzing four recent epigenome-wide association studies totaling more than 1000 postmortem prefrontal cortex brain samples using a uniform analytical pipeline. For each cohort, we employed two complementary analytical strategies, a sex-stratified analysis that examined methylation-Braak stage associations in male and female samples separately, and a sex-by-Braak stage interaction analysis that compared the magnitude of these associations between different sexes. Our analysis uncovered 14 novel CpGs, mapped to genes such as TMEM39A and TNXB that are associated with the AD Braak stage in a sex-specific manner. TMEM39A is known to be involved in inflammation, dysregulated type I interferon responses, and other immune processes. TNXB encodes tenascin proteins, which are extracellular matrix glycoproteins demonstrated to modulate synaptic plasticity in the brain. Moreover, for many previously implicated genes in AD neuropathology, such as MBP and AZU1, our analysis provided the new insights that they were predominately driven by effects in only one sex. These sex-specific DNA methylation differences were enriched in divergent biological processes such as integrin activation in females and complement activation in males. Our study implicated multiple new loci and biological processes that affected AD neuropathology in a sex-specific manner.


Introduction
Late-onset Alzheimer's disease (LOAD) is the most common cause of dementia. With the aging population in the U.S., Alzheimer's disease (AD) has become a major public health concern and one of the most financially costly diseases [1]. Almost two-thirds of AD patients in the U.S. are women [2]. After diagnosis, women also progress faster with more rapid cognitive and functional declines [3][4][5][6][7][8]. On the other hand, it has also been reported men with AD have an increased risk for death [9][10][11]. However, the molecular mechanisms underlying these observed disparities in AD are still not well understood. Previous studies have shown that epigenetics is an important contributor to the sex differences in brain functions and vulnerability to diseases [12][13][14][15][16]. Among epigenetic modifications, DNA methylation profiles differ significantly between males and females at many loci in adult brains [17]. Importantly, alterations of DNA methylation levels have also been implicated in multiple neurological disorders including AD [18][19][20][21][22].
However, thus far, a comprehensive characterization for the contribution of sex to DNA methylation differences in AD neuropathology has not been performed. In the identification of sex-specific effects, statistical power is a major challenge [23]. Stratifying by sex reduces the sample size of both groups. Also, comparing methylation to disease associations between the sexes by testing the interaction effect would require a much larger sample size than detecting the main effect with the same magnitude [24]. To address these challenges, we performed a comprehensive meta-analysis of more than 1000 postmortem brain prefrontal cortex (PFC) samples, collected from four recent AD epigenome-wide association studies [18][19][20][21], to identify the most consistent DNA methylation differences affected by AD neuropathology in a sex-specific manner. Within each cohort, to identify sexspecific differences in AD neuropathology, we employed two complementary approaches, a sex-stratified analysis that examined methylation-Braak stage associations in female and male samples separately, and a sex-by-Braak stage interaction analysis that compared the magnitude of these associations between different sexes. As sex is a strong factor in driving inter-personal variabilities in AD, the results of this study are particularly relevant for precision medicine.

Pre-processing of DNA methylation data
All the samples in this study were measured using the Infinium HumanMethylation450 BeadChip. Additional file 2: Table S2 shows the number of CpGs and samples removed at each quality control step. Quality control for CpG probes included removing probes with detection P-value < 0.01 in all samples of a cohort using the minfi R package, removing the 2623 CpGs associated with cigarette smoking identified in Joehanes et al. [25] (at P-value < 1 × 10 -7 threshold), and removing CpGs having a single nucleotide polymorphism (SNP) with minor allele frequency (MAF) ≥ 0.01 present in the last 5 base pairs of the probe using the DMRcate R package (with function rmSNPandCH and option dist = 5, mafcut = 0.01). For the quality control of samples, we removed samples with low bisulfite conversion efficiency (i.e., < 88%) or detected as outliers in principal component analysis (PCA). More specifically, PCA was performed using the 50,000 most variable CpGs for each cohort. Samples that were within ± 3 standard deviations from the means of PC1 and PC2 were selected to be included in the final sample set. The quality-controlled methylation datasets were next subjected to the QN.BMIQ normalization procedure as recommended by a recent systematic study of different normalization methods [26]. More specifically, we first applied quantile normalization as implemented in the lumi R package to remove systematic effects between samples. Next, we applied the β-mixture quantile normalization (BMIQ) procedure [27] as implemented in the wateRmelon R package [28] to normalize beta values of type 1 and type 2 design probes within the Illumina arrays. Finally, to remove batch effects, we applied the linear model methylation M value ~ methylation slide to M values of each CpG. The methylation residuals from these linear models were then used for subsequent analysis.

Single cohort analysis
To identify sex-specific DNA methylation differences in AD neuropathology, we performed both a sex-stratified analysis and a sex-by-Braak stage interaction analysis for each of the four brain sample cohorts. In the sexstratified analysis, we tested methylation-Braak stage associations in female and male samples separately. In sex-by-Braak stage interaction analysis, we analyzed both female and male samples simultaneously and compared slopes for methylation-Braak stage associations in females and males.
More specifically, in the sex-stratified analysis, for each CpG, we applied the model methylation residuals ~ age at death + Braak stage + estimated neuron proportions to female samples and male samples separately, where methylation residuals were obtained in "Pre-processing of DNA methylation data" described above. Here, the neuron proportion for each sample was estimated using the CETS R package [29], an R software that quantifies neuronal proportions from DNA methylation data using cell epigenotype specific (CETS) marks. In sex-by-Braak stage interaction analysis, for each CpG, we applied the model methylation residuals ~ age at death + sex + Braak stage + sex*Braak stage + sex*age at death + estimated neuron proportions to samples including both sexes.
For the analysis of differentially methylated regions (DMRs), we used the coMethDMR R package [30] to analyze 40,010 pre-defined genomic regions on the Illumina 450 k arrays and identify co-methylated DMRs associated with Braak stage. The pre-defined genomic regions are regions on the Illumina array covered with clusters of contiguous CpGs where the maximum separation between any two consecutive probes is 200 base pairs. First, coMethDMR selects co-methylated regions within these pre-defined contiguous genomic regions. Next, we summarized methylation M values within these co-methylated regions using medians and tested them against the AD Braak stage. The same linear models described for the analysis of CpGs were then applied to the median value of each DMR. We considered CpGs (or DMRs) with a false discovery rate (FDR) less than 0.05 in female samples or male samples to be significant.

Inflation assessment and correction
For genome-wide discoveries to be valid, the false positive rate of the study should be properly controlled. Because systematic inflation of test statistics can lead to an increase in the number of false-positive results, traditionally genomic inflation factor [31] is typically used to quantify the amount of inflation in genome-wide association studies (GWAS) of genetic variants. However, as shown by simulation studies [32], real datasets [32], and theory [31], the conventional genomic inflation factor (lambda or used interchangeably below) is dependent on the expected number of true associations. Because in a typical epigenome-wide association study (EWAS), it is expected that small effects from many CpGs might be associated with the phenotype, the genomic inflation factor would overestimate actual test-statistic inflation. To estimate genomic inflations more accurately in EWAS, Iterson et al. [32] developed a Bayesian method that estimates inflation in EWAS based on empirical null distributions, which is implemented in the Bioconductor package bacon.
In this study, to assess inflation of the test statistics, we used quantile-quantile (QQ) plots and estimated genomic inflation factors using both the conventional approach and the bacon method [32] (Additional file 1: Fig. S1). The bacon method was also used to obtain inflation-corrected effect sizes, standard errors, and P-values for each cohort, which were then combined by inversevariance weighted meta-analysis models using R package meta.

Meta-analysis
The evidence for heterogeneity of study effects was tested using Cochran's Q statistic [33]. The inverse-variance weighted fixed effects model was first applied to synthesize statistical significance from individual cohorts. Even though the fixed effects model for meta-analysis does not require the assumption of homogeneity [34], for the CpGs (or genomic regions) with nominal evidence for heterogeneity (nominal P heterogeneity < 0.05), we also applied random effects meta-analysis [35] and assigned final meta-analysis P-value based on the random effects model. For each CpG (and for each genomic region), we used the R package meta to obtain Braak stage effect in female samples and male samples separately in sex-stratified analysis, as well as meta-analysis P-values for sex-by-Braak stage interaction.

Identifying sex-specific differences
In the sex-stratified analysis, we selected significant CpGs (or genomic regions) with FDR < 0.05 in female samples or male samples separately. In sex-by-Braak stage interaction analysis, because the standard error of interaction effect sex × Braak stage is typically much larger than those for main Braak stage effects, the conventional approach for controlling false discovery rate often results in low power for discovering interaction effects [36]. Therefore, we used a stagewise analysis approach, previously proposed by van de Berge et al. (2017) [36], to help improve power in high-throughput experiments where multiple hypotheses are tested for each gene. More specifically, in the screening step, for each CpG (or genomic region), we tested the global null hypothesis that there is methylation-Braak stage association in either male or female samples. Next, in the confirmation step, we considered three individual null hypotheses for each CpG (or DMR): (a) there is no methylation-Braak stage association in male samples, (b) there is no methylation-Braak stage association in female samples, and (c) the methylation-Braak stage associations in male samples and female samples are the same. For the CpGs (or genomic region) selected in the screening step, these three individual hypotheses were then tested while controlling familywise error rate (FWER) as described in van de Berge et al. (2017) [36].
The stagewise analysis described above was implemented using the stageR package to identify CpGs (or genomic regions) with significant differential methylation-Braak stage associations in females and males. In the screening step, we considered meta-analysis P-values for the Braak stage in female samples and male samples (p.meta.female, p.meta.male), and used the minimum of these two meta-analysis P-values to represent each CpG (or genomic region). In the confirmation step, the parameter pConfirmation was defined using three P-values for each CpG (or genomic region): p.meta.female, p.meta. male, and p.meta.interaction (meta-analysis P-value for sex × Braak stage).

Enrichment and pathway analysis
The probes on the Illumina 450k array are annotated according to their locations with respect to genes (TSS1500, TSS200, 5′UTR, 1stExon, gene body, 3′UTR, intergenic) or CpG islands (island, shore, shelf, open sea). To understand the genomic context of sex-specific DNA methylation differences in AD neuropathology, we compared the FDR significant methylation differences from sex-stratified analysis with different types of genomic features. As Braak-associated methylation differences can occur at both significant individual CpGs and significant DMRs, we considered the CpGs located at significant individual CpGs or within significant DMRs jointly in this analysis, by testing their over-and under-representation in different types of genomic features using Fisher's exact test. More specifically, the proportion of significant CpGs mapped to a particular type of genomic feature (e.g., CpG islands) (foreground) was compared to the proportion of CpGs on the array that mapped to the same type of genomic feature (background).
Also, we used Fisher's test to assess enrichment of significant CpGs and DMRs in different chromatin states by comparing with the 15-chromatin state data for bulk PFC tissue samples (E073) from the Roadmap Epigenomics Project [37]. Using combinations of histone modification marks, ChromHMM [38] was previously used to annotate segments of the genome with different chromatin states (repressed, poised, and active promoters, strong and weak enhancers, putative insulators, transcribed regions, and large-scale repressed and inactive domains), which were shown to vary across sex, tissue type, and developmental age [39]. Similarly, we tested enrichment of significant CpGs and DMRs in binding sites of transcription factors and chromatin proteins from the ENCODE project [40] and CODEX database [41] using Locus Overlap Analysis as implemented in the LOLA R package [42].
Finally, we performed pathway analysis by comparing the genes with significant DNA methylation differences in AD neuropathology (identified in sex-stratified analysis) with the canonical pathways and biological process GO terms in MSigDB using Gene Set Enrichment Analysis (GSEA) [43]. First, we linked each CpG and each predefined genomic region tested in DMR analysis to genes by annotating them using the GREAT (Genomic Regions Enrichment of Annotations Tool) software, which associates genomic regions to target genes. With the default "Basal plus method", GREAT links each gene to a regulatory region consisting of a basal domain that extends 5 kb upstream and 1 kb downstream from its transcription start site, and an extension up to the basal regulatory region of the nearest upstream and downstream genes within 1 Mb [44]. Next, we represented each gene by the smallest P-value if there are multiple CpGs or genomic regions associated with them. To remove selection bias due to different numbers of CpGs or genomic regions associated with each gene (i.e., the smallest P-value for a gene with many CpGs or genomic regions linked to it is likely to be smaller than the smallest P-value for a gene with only a few linked CpGs or genomic regions), we next fit a generalized additive model [45] using the R package mgcv: Y i ∼ f (n.links i ) where Y i is negative log (base 10) transformation of the smallest P-value for gene i in the analysis of female samples (or male samples), n.links i is the number of CpGs or genomic regions linked to gene i, and f is a penalized spline function. We assumed gamma distribution for Y i , as under the null hypothesis of no association, Y i follows the chi-square distribution (a special case of gamma distribution). The residuals from this model were estimated and used to generate a ranked gene list, which was then used as input for GSEA (in preranked mode) to identify canonical pathways and gene ontology terms (MsigDB C2:CP and C5:BP collections of gene sets) enriched with significant methylation differences in female samples and male samples separately.

Integrative methylation and gene expression analysis
To systematically evaluate transcriptional differences near the observed sex-specific DNA methylation differences, we next performed integrative methylation-gene expression analysis using data on 333 female samples and 196 male samples from the ROSMAP study with matched DNA methylation and RNA-seq gene expression profiles measured on the prefrontal cortex. To this end, normalized FPKM (Fragments Per Kilobase of transcript per Million mapped reads) gene expression values for the ROSMAP study were downloaded from the AMP-AD Knowledge Portal (Synapse ID: syn3388564).
First, we linked significant CpGs (or DMRs) to nearby genes using GREAT [44], which associates genomic regions to target genes. Next, we removed confounding effects in DNA methylation data by fitting the model methylation M value ~ neuron.proportion + batch + sample.plate + ageAtDeath and extracting residuals from this model; these are the ROSMAP methylation residuals. Similarly, we also removed potential confounding effects in RNA-seq data by fitting model log2(normalized FPKM values + 1) ~ ageAtDeath + markers for cell types. The last term, "markers for cell types, " included multiple covariate variables to adjust for the multiple types of cells in the brain samples. More specifically, we estimated expression levels of genes that are specific for the five main cell types present in the CNS: ENO2 for neurons, GFAP for astrocytes, CD68 for microglia, OLIG2 for oligodendrocytes, and CD34 for endothelial cells, and included these as variables in the above linear regression model, as in previous large study of AD samples [18]. The residuals extracted from this model are the ROSMAP gene expression residuals.
Finally, for each gene expression and CpG (or DMR) pair, we then tested the association between gene expression residuals and methylation residuals using a linear model: ROSAMP gene expression residuals ~ ROSMAP methylation residuals + Braak stage. For significant DMRs, this analysis was repeated, except that methylation M value was replaced with median methylation M value from multiple CpGs in the DMR.

Sex-specific mQTL analysis
To identify methylation quantitative trait loci (mQTLs) for the significant DMRs and CpGs, we tested associations between the methylation levels with nearby SNPs, using the ROSMAP study dataset, which had matched genotype data and DNA methylation data for 434 female samples and 254 male samples. The ROSMAP genotype data was downloaded from AMP-AD (syn3157325) and imputed to the Haplotype Reference Consortium r1.1 reference panel [46]. There were two batches of genotype data, measured by Affymetrix GeneChip 6.0 (Affymetrix, Inc, Santa Clara, CA, USA) and Illumina HumanOmni-Express (Illumina, Inc, San Diego, CA, USA).
The male samples and female samples were analyzed separately. To reduce the number of tests, we focused on identifying cis mQTLs located within 500 kb from the start or end of the DMR (or position of the significant CpG) [47]. We additionally required SNPs to (1) have a minor allele frequency of at least 1%, (2) be imputed with good certainty: information metric (info score) ≥ 0.4, and (3) be associated with AD case-control status (as determined by clinical consensus diagnosis of cognitive status), after adjusting for age, batch, and the first 3 PCs estimated from genotype data, at nominal P-value less than 0.05. Next, for the ROSMAP methylation residuals obtained in section "Integrative methylation-gene expression analysis", we fit the linear model ROSMAP methylation residual ~ SNP dosage + batch + PC1 + PC2 + PC3, where PC1, PC2, and PC3 are the first three PCs estimated from genotype data, to test the association between methylation residuals in CpGs and the imputed allele dosages for SNPs to identify mQTLs. The analysis for DMRs was the same except that we replaced ROSMAP methylation residual with median (ROSMAP methylation residuals) of all CpGs located within the DMR. SNPs with FDR less than 0.05 in the linear model described above were considered to be significant mQTLs.

Drug target analysis
We compared our list of sex-specific DNA methylation differences with targets of drugs prescribed to AD patients or in the development for AD in the ChEMBL database [48] (https:// www. ebi. ac. uk/ chembl/). To this end, we overlapped genes mapped to significant CpGs or DMRs with the genes targeted by compounds annotated to "Alzheimer Disease" in ChEMBL.

Description of EWAS cohorts and data
Among the four cohorts (Additional file 2: Table S1), the mean age at death ranged from 79.3 to 87.2 years in females and from 67.5 to 85.0 years in males. The number of CpGs and samples removed at each quality control step are presented in Additional file 2: Table 2. For females, inflation factor lambdas ( ) by the conventional approach ranged from 1.060 to 1.154, and lambdas based on the bacon approach [32] ( bacon ) ranged from 1.021 to 1.059 (Additional file 1: Fig. S1). Similarly, for males, ranged from 0.906 to 1.265, and bacon ranged from 0.957 to 1.114. These values are comparable to those obtained in other recent large-scale EWAS [49].

Sex-specific DNA methylation differences in AD neuropathology
In the sex-stratified analysis, at 5% FDR, our meta-analysis identified 381 and 76 CpGs, mapped to 245 and 51 genes in female and male samples, respectively ( Fig. 1, Table 1, Additional file 2: Tables S3, S4). Similarly, at 5% FDR, we also identified 72 and 27 DMRs, mapped to 66 and 22 genes, in female and male samples, respectively ( Table 2, Additional file 2: Tables S5, S6). Among them, 3.6% (16 out of 441 unique FDR-significant CpGs) and 12.5% (11 out of 88 unique FDR-significant DMRs) were significant in both females and males with the same direction of change. The average number of CpGs per DMR was 6.5 ± 8.9. Also, the FDR-significant methylation differences at CpGs and DMRs did not completely overlap. Only 89 out of the 381 (23.4%) significant CpGs in females, and 13 out of the 76 (17.1%) significant CpGs in males overlapped with the significant DMRs. Among all CpGs and all DMRs, the effect estimates in males and females correlated only modestly (r CpG = 0.124, r DMR = 0.170, Additional file 1: Fig. S2), and about half (53% of CpGs, 54% of DMRs) were in the same direction of change in males and females, similar to what would be expected by chance.
In sex-by-Braak stage interaction analysis, we identified significant interaction at 14 CpGs, but no significant interactions at DMRs at 5% FDR. There was also little overlap between significant DNA methylation differences identified in sex-stratified and sex-by-Braak stage interaction analyses. Only 4 CpGs were identified by both analyses (Table 3). To understand this discrepancy, note that the sex-stratified analysis detected many differences that are attenuated but might be in the same direction in one sex group compared to the other. More specifically, among the FDR-significant CpGs identified in the sex-stratified analysis, many of them (370 out of 381 significant CpGs in the analysis of female samples, and 65 out of 76 significant CpGs in the analysis of male samples) had the same direction of change for methylation-Braak stage association in both sexes (Additional file 2: Tables S3, S4). In Table 1, among the 10 most significant CpGs from the sex-stratified analysis, 9 female-specific and 6 male-specific CpGs had the same direction of methylation-Braak stage association in both sexes. On the other hand, in sex-by-Braak stage interaction analysis, 13 out of the 14 significant CpGs had the opposite directions of differences for methylation-Braak stage associations in females and males (Table 3). Therefore, the interaction analysis was able to identify CpGs with large differences in sex-specific effect estimates, often in different directions, but these effects might not have reached FDR significance in sex-stratified analysis. For example, in Table 3, the CpG with the most significant interaction (cg13212831) had effect estimates of 0.083 and − 0.139 for females and males, respectively. In sex-stratified analysis, although the methylation-Braak stage associations were highly significant (P-value female = 0.006, P-valuemale = 4.1 × 10 -5 ), they did not reach 5% FDR significance threshold (FDR female = 0.413, FDR male = 0.097). Therefore, the results from sex-stratified analysis and sex-by-Braak stage interaction analysis complemented each other.  Table S7). Significant hypermethylated DMRs and CpGs in males are over-represented in CpG island shores, 5′UTRs, and TSS1500s. In contrast, significant hypomethylated differences in females and males are over-represented in open seas (Additional file 1: Fig. S3b, Additional file 2: Table S7). These observations are consistent with the knowledge that during aging, human brain DNA methylation levels gradually increase (hypermethylation) at genomic loci located at CpG islands and gene promoters, whereas intergenic CpG sites are marked with hypomethylation [50,51].

Enrichment analysis of sex-specific DNA methylation differences across genomic features
Our enrichment analysis for chromatin states showed that significant hyper-methylated differences in females were enriched in bivalent enhancer, flanking active TSS, repressed polycomb, and transcription at gene 5′ and 3′ regions (Additional file 1: Fig. S3c, Additional file 2: Table S8). On the other hand, significant hypermethylated differences in males were enriched in active TSS, flanking active TSS, and repressed polycomb regions. Hypomethylated differences in females were enriched in enhancers, weak repressed polycomb, and weak  Table S8). Notably, the enrichment of hypermethylated differences in repressed polycomb regions in both female and male samples is consistent with our previous sex-combined meta-analysis, which also highlighted the enrichment of hyper-methylated Braak-associated DNA methylation differences in polycomb repressed regions [22].
Similarly, enrichment tests for regulatory elements using the LOLA software also supported the potential functional relevance of these significant differences in DNA methylation. Significant DMRs and CpGs in females and males were both enriched in binding sites of EZH2 and SUZ12 (Additional file 2: Table S9), which are subunits of polycomb repressive complex 2 (PRC2), consistent with the observed enrichment of methylation differences in PRC2 repressed regions (Additional file 1:   S3c) and previous observations that DNA methylation often interact with PRC2 binding [52,53]. PRC2 is a type of polycomb group (PcG) protein and plays important roles in multiple biological processes including proliferation and differentiation as well as maintenance of cellular identity through regulation of gene expression. Of particular relevance to AD, PRC2 silences genes involved in neurodegeneration and its deficiency leads to the de-repression of developmental regulators such as the Hox gene clusters, which manifest in progressive and fatal neurodegeneration in mice [54].

Gene ontology and pathway analysis
Because of the relatively smaller number of gene sets being tested, a 25% FDR significance threshold, instead of the conventional 5% FDR, was suggested for GSEA [55]. At 25% FDR, the significant DNA methylation differences in females were enriched in TYROBP causal network (FDR = 0.014), TCR signaling in naïve CD4 + T cells (FDR = 0.130) and ROBO receptors bind AKAP5 (FDR = 0.160) gene sets, and significant methylation differences in males were enriched in the initial triggering of complement gene set (FDR = 0.245). The TYROBP causal network was previously inferred from a large-scale network analysis of human late-onset AD brains [56]; it was FDR significant (P-value < 0.001, FDR = 0.014) in females (Additional file 1: Fig. S4a), compared to a nominal association in males (P-value < 0.001, FDR = 0.620). Interestingly, the core enrichment subset of genes identified by GSEA in the female and male networks regulated by TYROBP involved DNA methylation differences at different genes (Additional file 1: Fig. S4b), highlighting different regulatory mechanisms for this gene network in males and females.
The comparison with gene ontology (GO) terms showed at 25% FDR, significant methylation differences in females were enriched in 25 GO biological processes ( Table 4, Additional file 2: Table S10), many of which are involved in inflammatory responses associated with AD pathology including CD8 positive alpha beta T cell activation and interferon alpha production, as well as other biological processes critical for AD pathogenesis such as response to platelet derived growth factor and positive regulation of axon extension. For males, we did not identify any significant GO terms at 25% FDR; the strongest enrichment with nominal P-value less than 0.001 involved immune responses to the accumulation of amyloid-β (Aβ) in the brain, such as regulation of T cell activation via T cell receptor contact with antigen bound to MHC molecule on antigen presenting cell, and other biological processes recently implicated in AD such as response to angiotensin [57,58] and cell redox homeostasis [59,60].

Correlation of sex-specific DNA methylation differences in AD neuropathology with expression levels of nearby genes
At 5% FDR, among FDR significant CpGs in females, all 381 CpGs were linked to a nearby gene by GREAT software, in which 14 were significantly associated with target gene expression levels (Additional file 2: Table S11), and half of them (n = 7) had effects in the negative direction. Among FDR significant CpGs in males, out of the 46 CpGs that were linked to a nearby gene, 2 were significantly associated with target gene expression and both were in the negative direction. Notably, in females, several of the most significant CpG methylation-gene expression associations were observed for the HLA-DPA1 gene, which encodes microglia receptors involved in antigen presentation and is regulated by PU.1 [61]. In males, the most significant CpG-gene expression was for HLA-DRB1, another PU.1 target gene [61]. For the 14 CpGs identified by our sex-by-Braak stage interaction analysis, only one CpG (cg24917065) was significantly associated with target gene (SLC25A37) expressions.

Correlation and overlap of sex-specific DNA methylation differences in AD neuropathology with genetic susceptibility loci
To evaluate if the significant methylation differences are located in the vicinity of sex-specific genetic variants implicated in AD, we compared our sex-specific CpGs and DMRs with the recently identified sex-specific SNPs associated with AD biomarkers [62] or AD pathology [63]. We found only 5 CpGs, mapped to the SERP2, KCNE1, TNKS1BP1, FAM165B, PLCB4 genes were located within 500 kb of the sex-specific SNPs (Additional file 2: Table S12).
To search for mQTLs, which are genetic variants associated with DNA methylations, we next tested associations between the sex-specific CpGs and DMRs with SNPs that are located within 500 kb from them using 434 female samples, 254 male samples from the ROSMAP study, which had both genotype and DNA methylation data. While no mQTL-DMR pairs reached 5% FDR significance, we did identify 572 and 284 FDR-significant mQTL-CpG pairs associated with the sex-specific CpGs in females and males, respectively (Additional file 2: Tables S13, S14). Among the 381 and 76 sex-specific CpGs identified in female and male samples, respectively, 41 (11%) and 15 (20%) had at least one corresponding mQTL in brain samples. Among the 14 CpGs identified in our sex-by-Braak stage interaction analysis, 2 and 7 CpGs with at least one brain mQTL, corresponding to 21 and 236 significant mQTL-CpG pairs, were identified at 5% FDR in females and males, respectively (Additional file 2: Table S15). These mQTL-CpG pairs point to important potential molecular mechanisms of diseaseassociated genetic variants that might due, at least in part, to their influences on DNA methylations, which can be further validated in mechanistic studies.

Drug target analysis of sex-specific DNA methylation differences
To investigate the clinical impact of the sex-specific DNA methylation differences, we next compared them with targets of drugs in the ChEMBL database [48] that are annotated to Alzheimer's disease, many of which are antipsychotic medications commonly prescribed to AD patients for treating psychiatric symptoms that accompany AD. We found that 13 CpGs and 2 DMRs, mapped to 20 genes, had overlap with targets of 16 different drugs (Additional file 2: Table S16). Among them, CACNA1C encodes a voltage-dependent calcium channel, which is a target of cholinesterase inhibitor donepezil. Previously, drug responses for donepezil were shown to be modulated by the sex hormone estrogen receptor alpha (ESR1) genotype [64]. Several CpGs and one DMR are mapped to targets of valproic acid, a mood stabilizer often prescribed for AD patients and was shown to have different pharmacokinetic profiles between male and female subjects [65]. Interestingly, two CpGs and 1 DMR also mapped to targets of caffeine, which was included in cocktail therapy in AD clinical trials [66,67]. Although caffeine reduces the risk for AD [68,69] in both men and women, the protective effect seems to be greater in women [70]. Also, CHRM3 encodes muscarinic acetylcholine receptor, which is targeted by two commonly prescribed antipsychotic drugs for AD patients, trazodone and haloperidol. In both human and animal models, it has been observed treatment with haloperidol induces sex-specific DNA methylation differences [71,72]. While this study could not establish the mechanisms at which DNA methylation interacts with drugs that AD patients take, we hypothesize that these might include the influence of DNA methylation on drug responses or the differences in DNA methylation resulted from drug actions.

Discussion
To identify sex-specific differences in AD neuropathology, we performed a sex-stratified analysis and a sex-by-Braak stage interaction analysis for each cohort and then used a meta-analysis strategy to combine the cohort-specific association signals. In the sex-stratified analysis, as discussed above, a substantial number of the significant loci showed the same direction but attenuation of effect size for methylation-Braak stage association in a different sex (Tables 1, 2). Therefore, it is not surprising that many of these significant CpGs were identified previously in sex-combined meta-analysis [22]. Among FDR significant methylation differences in females, 325 CpGs (85%) and 40 DMRs (56%), mapped to genes such as HOXA3, AZU1, and MBP were also previously identified in our sex-combined meta-analysis [22] (Additional file 2:  Tables S3, S5). Similarly, in the analysis of male samples, among FDR significant differences, 58 (76%) CpGs and 15 DMRs (56%), mapped to genes such as MAMSTR, RHBDF2, and AGAP2, overlapped with significant hits from sex-combined meta-analysis [22] (Additional file 2:  Tables S4 and S6). However, our sex-specific analysis provided the new insight that the effects of these known AD genes appear to be predominately driven by effects in only one sex (Tables 1, 2).
On the other hand, our sex-specific analysis also uncovered novel methylation differences at 84 CpGs and 42 DMRs that were not identified previously by sex-combined analyses [22], which may had reduced power due to heterogeneity between the sexes. For example, among the top 10 CpGs in the sex-stratified analysis (Table 1), a new locus at cg22632947, which mapped to the gene body of the PRKCA gene, was highly significant in female samples (estimate = − 0.139, P-value = 1.50 × 10 -7 , FDR = 3.00 × 10 -3 ), but not significant in male samples (estimate = − 0.005, P-value = 0.857, FDR = 0.995) (Additional file 1: Fig. S5). The PRKCA gene encodes protein kinase Cα (PKCα), which participates in synaptic loss resulting from the accumulation of amyloid-β (Aβ) in AD neuropathology [73,74]. Another novel locus is at cg18942110 in the promoter of the CRTC3 gene, where methylation-Braak stage association was highly significant in male samples (estimate = − 0.164, P-value = 2.23 × 10 -6 , FDR = 3.19 × 10 -2 ), but not significant in female samples (estimate = − 0.031, P-value = 0.306, FDR = 0.952) (Additional file 1: Fig. S5). CRTC3 is a member of the CRTC family, which are coactivators of the transcription factor CREB (cAMPresponse element binding protein). In addition to its crucial role in maintaining synaptic plasticity and facilitation of short-term memory to long-term memory, the CREB signaling pathway also mediates synapse loss induced by Aβ in AD [75]. Notably, synapse loss significantly correlates with cognitive impairment [76,77] and has been observed to be an early feature of AD pathogenesis [78,79].
The sex-by-Braak stage interaction analysis also uncovered several additional novel methylation loci that affected AD neuropathology in a sex-specific manner. Notably, none of the 14 CpGs detected in our interaction analysis was identified in previous large-scale DNA methylation studies [18][19][20][21][22], suggesting that sex-specific differences such as these can be missed by conventional studies that do not consider the impact of sex. This is likely due to the cancelation of effects in sex-combined analysis because the majority of these 14 CpGs had different directions of methylation-Braak stage effects in male and female samples (Table 3). Among genes mapped to these 14 CpGs, TMEM39A is a member of the transmembrane (TMEM) protein family. In recent GWAS, a genetic variant on TMEM39A was discovered and replicated as an important risk locus for multiple sclerosis, an autoimmune condition of the central nervous system [80,81]. While relatively little is known about the role of TMEM39A in AD, given its important contributions to inflammation, dysregulated type I interferon responses, and other immune processes [82] which are also implicated in AD, methylation differences affecting this gene are particularly relevant. Another noteworthy gene is TNXB and its pseudogene TNXA, which are located in the major histocompatibility complex (MHC) class III region on chromosome 6. TNXB encodes tenascin proteins, which are extracellular matrix glycoproteins demonstrated to modulate synaptic plasticity in the brain [83]. In particular, genetic variants at the HLA-DQB1 locus discovered in the recent AD genetic meta-analysis [84] included eQTLs for TNXB/TNXA in brain tissues [84,85].
Consistent with previous studies [18,19,86,87], we observed the majority of these sex-specific differences were hyper-methylated in samples with AD neuropathology, for which methylation levels increased as the AD Braak stage increased (Additional file 2: Tables S17, S18). More specifically, 59% of the significant CpGs and 69% of the significant DMRs in females, along with 66% of the significant CpGs and 89% of the significant DMRs in males were hyper-methylated in samples with AD neuropathology (Additional file 2: Tables S3-S6).
To better understand the relevance of these Braakassociated sex-specific differences, we also compared our results with several previous studies. The comparison with Xia et al. [16] and Xu et al. [17], which examined differential methylation between males and females in the prefrontal cortex, but without considering AD neuropathology [16,17], showed our results were largely distinct. Among 451 unique CpGs identified in our sex-stratified analysis or sex-by-Braak stage interaction analysis, only 16 were also identified in Xia et al. [16] and none were identified in Xu et al. [17] (Additional file 2: Tables S3-S6). This is probably due to different hypotheses tested in our study and the sexual dimorphism studies -while our study examined the impact of sex on methylation-Braak stage association, the previous studies examined differential methylation between the sexes, regardless of AD neuropathology. The comparison of our results with sex-specific DNA methylation differences in fetal brain development [88,89] also showed very little overlap (Additional file 2: Table S19); one hypothesis could be that the Braak-associated sex-specific DNA methylation differences identified in this study might be influenced by environmental risk factors for AD, such as diet and exercise.
The results of our gene set analysis highlighted a number of critical sex-specific biological processes in AD neuropathology. Notably, the TYROBP causal network reached the FDR significance threshold in females (FDR = 0.014) but was only nominally significant in males. Interestingly, Braak-associated CpG methylation differences that drove pathway associations (core enrichment genes) occurred at different genes in females and males (Additional file 1: Fig. S4), indicating a potentially sex-specific regulatory mechanism for this network. TYROBP (TYRO protein tyrosine kinase-binding protein) is a key regulator of the complement pathway in the immune/microglia network, which is activated as Aβ accumulates in LOAD brains [56,90]. TYROBP is a transmembrane adaptor protein for TREM2, SIRPβ1, and CR3 receptors, which are known to be involved in AD pathogenesis [90][91][92]. Also, TYROBP is regulated by SPI1, a central hub for the network of genes involved in myeloid immune response in neurodegeneration [93]. In patients with LOAD, TYROBP was observed to be up-regulated in the brains in multiple cohorts [56]. Recent studies suggested TYROBP-mediated signaling is involved in multiple important functions as aggregating Aβ activates microglia, including enhanced phagocytosis of damaged neurons [56,90] and suppression of inflammatory responses [94], as well as neuronal pruning activity [56]. Interestingly, in gene ontology (GO) analysis, among the most significant GO Biological Process terms (P-value < 0.001) in females and males, none of them overlapped (Additional file 2: Table S10), even though the relevancy of all the top biological processes were supported by recent AD literature (Table 4). These results suggest different biological processes are associated with AD pathology in males and females.
Importantly, a number of these sex-specific biological processes pointed to important potential biomarkers and therapeutic targets for the treatment of AD. For example, one of the top biological processes enriched with significant methylation differences in female samples is response to platelet derived growth factor. Recently, multiple studies have shown that reduced levels of plateletderived growth factors (PDGFs) in plasma significantly correlate with mild cognitive impairment and have proposed PDGFs as a potential biomarker for AD [95,96]. For the significant methylation differences in male samples, one of the top biological processes highlighted by our enrichment analysis is dysregulation in the complement system. Recently, a number of novel agents targeting the complement system are being developed and tested in clinical trials for potential effective therapy for AD [97]. Therefore, clinical trials testing potential treatment for AD patients might have more power for detecting treatment effects by considering sex and targeting the subgroup with the higher predicted benefit based on patient molecular profiles such as DNA methylation.
There are several limitations for this study. The methylation levels in the studies analyzed here were measured on the bulk prefrontal cortex, which contains a complex mixture of cell types. To reduce confounding effects due to cellular heterogeneities, we included the estimated neuron proportion of each brain sample as a covariate variable in all our analyses. Currently, a challenge with cell-type-specific studies is that they are often limited to smaller sample sizes due to labor-intensive sample preparation procedures and therefore have limited statistical power. Also, we did not identify any CpGs or DMRs from chromosome X, this might suggest that sex differences in AD neuropathology are not primarily due to chromosome X. Alternatively, the lack of association might also be due to the limited coverage by the 450k array. Future studies utilizing high throughput sequencing that provides better coverage of the epigenome will help clarify the role of the X chromosome in AD neuropathology. Finally, the associations we identified in this study do not necessarily reflect causal relationships. Future studies that employ longitudinal designs are needed to identify causal changes in DNA methylation as AD initiates and progresses.
In summary, our study highlighted the importance of stratifying on sex and analyzing sex-by-disease interaction in the analysis of DNA methylation data to discover the epigenetic architectures underlying AD neuropathology. Our meta-analysis discovered many novel sex-specific DNA methylation differences consistently associated with the AD Braak stage in multiple studies. Because of cancelation of effects in different directions, or dilution from samples with no effect, these sex-specific effects would be missed by sex-combined analysis. Moreover, for many genes previously linked to AD neuropathology, our work provided evidence that the DNA methylation differences at these genes were predominately driven by effects in only one sex. Our enrichment analysis highlighted divergent biological processes in males and females, which underscored sex-specific regulatory mechanisms involved in AD neuropathology. Finally, our results also have important implications for precision medicine-many of the sex-specific DNA methylation differences also pointed to important potential AD biomarkers and therapeutic targets, suggesting a pressing need for developing and applying sex-specific treatment strategies for AD.
Additional file 1: Supplementary Figures. Figure S1. Quantile-quantile (QQ) plots of observed and expected distributions of P-values in Gasparoni, London, Mount Sinai, and ROSMAP cohorts. Figure S2. Comparison of methylation-Braak stage associations in female samples and male samples. Figure S3. Enrichment of FDR-significant CpGs and CpGs located within FDR-significant DMRs with positive and negative effect estimates in various genomic features and chromatin states. Figure S4. Gene Set Enrichment of the TYROBP causal network with sex-specific Braak-associated DNA methylation differences. Figure S5. Forest plots for several top CpGs identified in sex-stratified analysis and sex-by-Braak stage interaction analysis.
Additional file 2: Supplementary Tables. Table S1. Sample information of the brain samples datasets included in the meta-analysis. Table S2.
Quality control (QC) information on DNA methylation samples and probes for each cohort contributing to this meta-analysis. Table S3. In females, a total of 381 differentially methylated CpGs were significantly associated with the AD Braak stage at 5% FDR in the meta-analysis of four brain samples cohorts (Gasparoni, London, Mt. Sinai, ROSMAP). Table S4. In males, a total of 76 differentially methylated CpGs were significantly associated with the AD Braak stage at 5% FDR in the meta-analysis of four brain samples cohorts (Gasparoni, London, Mt. Sinai, ROSMAP). Table S5. In females, at 5% FDR, a total of 72 co-methylated DMRs were significantly associated with AD Braak stage in the meta-analysis of the four brain samples cohorts. Table S6. In males, at 5% FDR, a total of 27 co-methylated DMRs were significantly associated with AD Braak stage in the meta-analysis of four brain samples cohorts. Table S7. Enrichment of sex-specific DNA methylation changes in different genomic features. Table S8. Enrichment of sex-specific DNA methylation changes in different chromatin states. Table S9. Enrichment of sex-specific DNA methylation changes in binding sites of transcription factors and chromatin proteins assayed by ENCODE and CODEX projects. Table S10. Results of Gene Set Enrichement analysis of sex-specific CpGs and DMRs in female and males. Table S11. Significant association between sex-specific differences and expression levels of nearby genes. Table S12. The sex-specific CpGs and DMRs located within 500 kb of sex-specific SNPs associated with AD biomarkers or AD neuropathology. Table S13. mQTLs assoicated with FDR-significant CpGs in females. Table S14. mQTLs assoicated with FDR-significant CpGs in males. Table S15. mQTLs associated with significant CpGs identified by sex-by-Braak stage interaction analysis. Table S16. Sex-specific DNA methylation differences overlapped with AD drug targets in ChEMBL database.