Sample overview and genome-wide DNA methylation variation
Table 1 provides a demographic and clinical description of the cohort. As intended, the ADHD group had significantly higher clinical symptoms, and as expected, had higher rates of psychiatric comorbidity, and more lifetime exposure to psychiatric medications. In the full sample, the ADHD group had slightly lower family income than controls (p = 0.03), but the difference was not significant in the European-ancestry subsample used for the PRS analyses. As is common in studies of ADHD, the ADHD group had a slightly lower IQ (WISC-IV FSIQ mean difference = 7.1; p = 5.6e–10), but both groups had mean IQ in the average range. The groups were of similar age, although age was included as a covariate as a precaution. Boys were overrepresented in the ADHD group and we adjusted all results for sex.
We calculated a predicted age for each child, based on age-associated DNA methylation sites10. As expected, reported ages and predicted ages were highly correlated (r = 0.58, p = 1.3e–54), providing reassurance as to the validity of the DNA methylation data obtained from saliva. There was no evidence for differential epigenetic age acceleration in the ADHD group (p = 0.98) (Fig. S4). We also tested for a “global” difference in DNA methylation levels between ADHD cases and controls, determined by averaging cell-type-corrected beta values across all 568,281 probes (see the “Methods” section, Eq. (1)), and found no significant difference (p = 0.74). This result indicates that, as expected, ADHD is not associated with any systemic methylomic differences across probes included on the Illumina EPIC array.
DNA methylation associated with ADHD
Differential DNA methylation between ADHD cases and controls at individual sites was investigated by regressing methylation levels at each probe on ADHD status (Eq. (2)) with our standard covariates. EWAS results were well-controlled with an inflation factor λ = 1.02 (calculated as the ratio of the median observed log10(p-value) to the median expected log10(p-value))53. A Q–Q plot and Manhattan plot are presented in Fig. S5.
No DMPs passed EWAS significance (p = 8.8e–8), which corrects for all 568,281 probes tested. Table 2 shows the 7 DMPs associated with ADHD at our a priori suggestive significance level of p < 1e–5. Methylation differences between ADHD cases and controls ranged between 0.3% and 1.4%. The top 100 probes associated with ADHD are provided in Table S5 for descriptive purposes. The top-ranked probe, cg17478313 (p = 1.54e–6), shows higher DNA methylation in ADHD cases compared with controls (Δβ = 0.93%) and is located in the promotor region of SLC7A8 (Figs. 1a and 2a). Probe cg21609804 (p = 2.82e–6), also with higher methylation (Δβ = 1.38%) in ADHD cases, is located in the 3′-UTR of MARK2 (Fig. 1b). Additional suggestive DMPs were annotated to PDLIM5, VPS28, ZNF706, and FAM59A.
DNA methylation associated with ADHD polygenic risk
We have previously shown that the ADHD PRS is associated with ADHD status (explaining ~4.5% of the variance) in this sample4. Here we use the ADHD PRS to investigate the relationship between overall ADHD genetic burden and DNA methylation (Eq. (3)). Again, EWAS results for the PRS are well-controlled (λ = 1.08). A Q–Q plot and Manhattan plot are presented in Fig. S6. One probe met our genome-wide significance threshold, cg15472673 (p = 6.71E–8), characterized by reduced DNA methylation with higher PRS (Fig. 2b). The association remains (p = 9.76e–8), when including ADHD status in the regression model, indicating that the effect is not driven by elevated polygenic burden in ADHD cases. This probe is located in a CpG island of a bivariate promoter between the GART and SON genes. None of the SNPs included in the PRS are direct mQTLs (see below) for cg15472673 (all mQTL p-values > 1e–5), indicating that the association with the PRS is not driven by a simple genetic effect on DNA methylation. DNA methylation levels at 12 other probes were associated with the PRS at p < 1.0e–5, 10 of which showed increased methylation with higher PRS (Table 2b). A summary of all findings for ADHD diagnosis and the ADHD PRS is presented in Fig. S7.
Sex-specific variation in DNA methylation
Because we previously reported an association between ADHD and DNA methylation at sites annotated to VIPR2 and MYT1L specifically in boys31, we first examined sex-by-diagnosis interaction effects among all probes annotated to these two genes (239 probes) at an adjusted significance threshold of p < 0.0002 (0.05/239). The two strongest sex-by-diagnosis interactions annotated to VIPR2 were for cg26975193 and cg20998127 (sex-by-diagnosis interaction p = 7.51e–6 and p = 0.000459). Supporting the finding from our previous study, males with ADHD had lower methylation compared with male controls at both sites (Δβ = −0.22%, p = 0.0185; Δβ = −0.51%, p = 3.99e–5). However, among females, ADHD cases had higher methylation levels than controls at cg26975193 (Δβ = 0.35%, p = 0.00597) and were not significantly different from controls at cg20998127 (Δβ = 0.14%, p = 0.384). To ensure replication of the hypomethylation effect in boys with ADHD reported by Wilmot et al.31, we removed the samples included in that previous report (n = 73 that survived QC for the current probe set). Effects among males were consistent at both sites (for cg26975193: Δβ = −0.19%, p = 0.0819; for cg20998127: Δβ = −0.51%, p = 4.14e–4).
For MYT1L, our previously reported finding did not survive multiple-testing correction (minimum sex-by-diagnosis interaction p-value = 0.0039 for cg02870147; males Δβ = −0.49%), although the direction of effect at this site (lower methylation among males with ADHD) was consistent with our previous data, even with prior samples (n = 73) removed (males Δβ = −0.42%). Likewise, no probes annotated to MYT1L showed main effects that passed the Bonferroni threshold (minimum p-value of 0.00395 for cg22140907; Δβ = −0.31%). All results from the analysis of these genes are reported in Table S6.
We next performed an EWAS of sex-specific DMPs associated with both ADHD and the ADHD PRS. No sex interactions were significant at the EWAS-wide threshold, although nine probes show sex-by-diagnosis interactions (Table 2c) and five probes show sex-by-PRS interactions (Table 2d) at our suggestive significance threshold (p < 1e–5).
Differentially methylated regions
The results from all DMP analyses were used to investigate differentially methylated regions (DMRs) using the Comb-p software tool41. A single significant DMR on chromosome 6 was identified (Šidák-corrected p = 3.4e–5), which contained eight probes associated with the ADHD PRS in a sex-specific manner (Fig. S8). Specifically, among females a higher PRS was associated with higher methylation levels, and an opposite (though much weaker) relationship was seen among males. This DMR (chr6: 31148383–31148553) lies within the major histocompatibility complex, ~3 kb upstream of PSORS1C3.
DNA methylation quantitative trait loci
We identified methylation quantitative trait loci (mQTL) associated with all DMPs suggestively associated with (a) ADHD or (b) the ADHD PRS (Table 2). DNA methylation at the two top-ranked ADHD-associated DMPs (cg17478313 annotated to SLC7A8 and cg21609804 annotated to MARK2, Table 2a) was significantly associated with genotypes at nearby SNPs (Fig. 3). The SLC7A8 probe (cg17478313) was associated with SNP rs7141505 (p = 1.2e–46), located in the gene’s promoter region. The MARK2 probe (cg21609804) was associated with an intronic SNP rs928948 (38 kb upstream of the CpG) at p = 2.07e–59 (Fig. 1b). Furthermore, for both mQTLs, the relationship between genotype and DNA methylation level is similar in both ADHD cases and controls (Fig. 3b, d). Also, for both cg17478313 and cg21609804, both ADHD status (p = 2.2e–4 and p = 1.85e–4) and genotype (p = 4.4e–45 and p = 8.9e–58) are associated with methylation levels when included together in the regression models. These results indicate that the mQTL effects are not simply due to allele frequency differences between ADHD cases and controls. An additional 279 SNPs are associated with methylation levels at these two DMPs (31 for cg17478313 and 248 for cg21609804), and pass our experiment-wide mQTL significance threshold (p < 1.8e–10).
For the ADHD PRS, of all the associated DMPs (at p < 1e–5), significant mQTLs were found only for cg04453792 annotated to USP31 (p < 1.8e–10). All mQTL analysis results are provided in Table S7.
To determine whether mQTL variants associated with the top-ranked DMPs (for both ADHD and the ADHD PRS) were also implicated in GWAS for ADHD (i.e., pleiotropic for methylation and ADHD), we examined the results from the PGC ADHD GWAS meta-analysis3 in the regions of SLC7A8, MARK2, GART, SON, USP31, and LOC100130015 (Table 2 and Table S7). We assessed all GWAS SNPs included in our mQTL analysis that were in linkage disequilibrium with these six genes or their putative regulatory regions (20 kb upstream or downstream of each gene, see Supplementary Materials). No genome-wide significant (or suggestive) ADHD-associated SNPs assigned to MARK2 (minimum p = 0.02), SLC7A8 (minimum p = 0.3), GART/SON (minimum p = 0.029), USP31 (minimum p = 0.069), or LOC100130015 (minimum p = 0.011) were seen in the GWAS meta-analysis, indicating there is no evidence of pleiotropy for these mQTLs.
To expand the investigation of pleiotropy, all mQTLs observed in our cohort (not just those associated with DMPs) were tested for colocalization with variants associated with ADHD in the recent ADHD GWAS meta-analysis conducted by the PGC3. Of the 12 ADHD-associated regions identified in the GWAS meta-analysis, 11 contained genome-wide significant cis-mQTLs in our cohort (Table S3). Evidence for colocalization/pleiotropy was found for variants in 5 of the 12 ADHD-associated regions (Table S4). For two of the regions, 12q21.33 and 15q21.1 (Figs. S9 and S10), both the SMR and coloc methods identify the same causal SNPs (rs2279574, SMR p = 2.3e–8, coloc posterior probability = 0.98, and rs1656622, SMR p = 3.6e–5, coloc posterior probability = 0.82). SNP rs2279574 is a missense variant within DUSP6, and rs1656622 lies within an intron of SEMA6D. All significant mQTLs in the 12 ADHD-associated regions are included in Table S10.
Gene set enrichment with DM probes
For exploratory and hypothesis-generating purposes, we report gene set findings. DMPs nominally associated with the ADHD PRS (p < 0.001, the recommended default in the methylGSA package)49 are annotated to 91 genes with elevated expression (five-fold higher-than-average expression in all other tissues, see Methods) in the brain (unadjusted enrichment p-value = 0.0097). Several (Fig. S11) relate to ion channels (e.g., KCNIP1, KCNK10, CACNA1E, and CACNB4) or are involved in cell adhesion (e.g., NCAM2, NRXN1, CNTNAP2, and CDH22), both previously implicated in ADHD or other psychiatric traits54,55,56,57,58. No GO categories were significantly enriched with either ADHD-associated or ADHD PRS-associated DMPs after multiple testing correction. All enrichment analysis results are shown in Table S8.