Convergence of evidence from a methylome-wide CpG-SNP association study and GWAS of major depressive disorder

Translational Psychiatry, Aug 2018

DNA methylation is an epigenetic modification that provides stability and diversity to the cellular phenotype. It is influenced by both genetic sequence variation and environmental factors, and can therefore potentially account for variation of heritable phenotypes and disorders. Therefore, methylome-wide association studies (MWAS) are promising complements to genome-wide association studies (GWAS) of genetic variants. Of particular interest are methylation sites (CpGs) that are created or destroyed by the alleles of single-nucleotide polymorphisms (SNPs), as these so-called CpG-SNPs may show variation in methylation levels on top of what can be explained by the sequence variation. Using sequencing-based data from 1132 major depressive disorder (MDD) cases and controls, we performed a MWAS of 970,414 common CpG-SNPs. The analysis identified 27 suggestively significant (P < 1.00 × 10−5) CpG-SNPs associations. Furthermore, the MWAS results were over-represented (odds ratios ranging 1.36–5.00; P ranging 4.9 × 10−3–8.1 × 10−2) among findings from three recent GWAS for MDD-related phenotypes. Overlapping loci included, e.g., ROBO2, ASIC2, and DCC. As the CpG-SNP analysis accounts for the number of alleles that creates CpGs, the methylation differences could not be explained by differences in allele frequencies. Thus, the results show that the MWAS and GWASs provide independent lines of evidence for the involvement of these loci in MDD. In conclusion, our methylation study of MDD contributes novel information about loci of relevance that complements previous findings and generates new hypothesis about MDD etiology, such as that the functional effects of genetic association may be partly mediated and/or enhanced by the methylation status in these loci.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://www.nature.com/articles/s41398-018-0205-8.pdf

Convergence of evidence from a methylome-wide CpG-SNP association study and GWAS of major depressive disorder

Abstract DNA methylation is an epigenetic modification that provides stability and diversity to the cellular phenotype. It is influenced by both genetic sequence variation and environmental factors, and can therefore potentially account for variation of heritable phenotypes and disorders. Therefore, methylome-wide association studies (MWAS) are promising complements to genome-wide association studies (GWAS) of genetic variants. Of particular interest are methylation sites (CpGs) that are created or destroyed by the alleles of single-nucleotide polymorphisms (SNPs), as these so-called CpG-SNPs may show variation in methylation levels on top of what can be explained by the sequence variation. Using sequencing-based data from 1132 major depressive disorder (MDD) cases and controls, we performed a MWAS of 970,414 common CpG-SNPs. The analysis identified 27 suggestively significant (P < 1.00 × 10−5) CpG-SNPs associations. Furthermore, the MWAS results were over-represented (odds ratios ranging 1.36–5.00; P ranging 4.9 × 10−3–8.1 × 10−2) among findings from three recent GWAS for MDD-related phenotypes. Overlapping loci included, e.g., ROBO2, ASIC2, and DCC. As the CpG-SNP analysis accounts for the number of alleles that creates CpGs, the methylation differences could not be explained by differences in allele frequencies. Thus, the results show that the MWAS and GWASs provide independent lines of evidence for the involvement of these loci in MDD. In conclusion, our methylation study of MDD contributes novel information about loci of relevance that complements previous findings and generates new hypothesis about MDD etiology, such as that the functional effects of genetic association may be partly mediated and/or enhanced by the methylation status in these loci. Introduction Major depressive disorder (MDD) is a complex disorder that is characterized by persistent dysphoria and is often accompanied by considerable morbidity1,2,3 and mortality2,4. Because MDD has a lifetime prevalence of almost 15%5, tends to start early in life6, and is often chronic7,8, it is the leading contributor to disability worldwide9,10. In comparison with other (psychiatric) disorders, discerning the biological basis of MDD has been difficult. Only very recently, a number of genetic variants were identified and replicated11,12. However, these variants had small effect sizes and explained only a small proportion of the disease risk. DNA methylation is an epigenetic modification that provides stability and diversity to the cellular phenotype. Because methylation is dynamic in nature and can be altered by environmental factors, it can potentially account for key clinical features of MDD such as its episodic nature or mediate the effects of environmental risk factors such as stress13,14,15. Therefore, methylome-wide association studies (MWAS), which test a genome-wide set of methylation sites for association with an outcome of interest, are promising complements to genome-wide association studies (GWAS) of genetic variants. Of particular interest are methylation sites (CpGs) that are created or destroyed by single-nucleotide polymorphisms (SNPs). These sites, commonly referred to as CpG-SNPs, may show variation in both methylation and sequence, and may therefore convey information beyond either of the two data types alone. However, methylation-dependent association signals in CpG-SNPs are not captured by GWAS and are very poorly captured by a regular MWAS. Therefore, a specific CpG-SNP analysis is needed to detect these signals. Regular GWAS studies detect differences in allele frequencies between cases and controls. In contrast, a CpG-SNP analysis tests whether groups of cases and controls with the same genotype show differences in methylation at these sites. Thus, these two analyses capture different signals. Similarly, while a regular MWAS detects differences in methylation it does not account for differences in genotype and will therefore often lack the statistical power to detect association signals for CpG-SNPs. A CpG-SNP MWAS remedies this by including information on the actual genotypes of each subject. Furthermore, the link between sequence variation and methylation levels at these sites may allow CpG-SNPs to function as important cis-regulatory polymorphisms that connect genetic variation to variation in methylation. For example, the alleles present and the methylation levels observed at a specific CpG-SNP have been associated with a variety of regulatory functions16,17,18,19. In addition, in a high-density analysis of methylation quantitative trait locus (meQTL), CpG-SNPs were involved in the majority of all identified meQTLs20. To study whether CpG-SNPs contribute to MDD disease risk, we used a sequencing-based approach that provides nearly complete coverage of all CpGs21,22, including close to 1 million CpG-SNPs. To further explore the MWAS findings and their potential relevance for MDD, we also tested for their overlap with results from recent GWASs. Materials and methods Description of the NESDA sample DNA from blood was obtained from 1200 individuals from the Netherlands Study of Depression and Anxiety (NESDA). MDD was diagnosed using the DSM-IV-based Composite International Diagnostic Interview (CIDI version 2.1) that was administered by specially trained research staff23. In addition, to a current MDD diagnosis, cases had a score >14 on the IDS-SR3024, a 30-item self-report measure of depression symptoms. Controls had no lifetime psychiatric disorders and an IDS-SR30 score <14. The sample selection was further based on good quality GWAS genotype information available from a previous investigation25 (for a summary description, see the Supplementary Note). For further details about NESDA, and demographic and clinical characteristics of participants used for the present study, see Table S1. The study was approved by the ethical committees of all participating locations, and participants provided written informed consent. Assaying the methylome with MBD-Seq We assayed the methylome using an optimized protocol for methyl-CG binding domain sequencing (MBD-Seq) that provides almost complete coverage of all CpGs in the genome21. In short, we used ultrasonication to shear genomic DNA into, on average, 150 bp fragments followed by enrichment with MethylMiner™ (Invitrogen) to capture the methylated fraction of the genome. The captured fragments were eluted and used to create a barcoded sequencing library for each methylation capture. Labeled sequencing fragment libraries were pooled in equal molarities and sequenced on a NextSeq500 instrument (Illumina). To ensure consistency in the sample preparation, MethylMiner captures and library constructions were both performed using Biomek NxP robotics (Beckman Coulter). Samples were performed in a randomized order and all labtecnical procedures were performed blind to any phenotype information. The sequence reads were aligned to the human reference genome (hg19/GRCh37) using Bowtie226. Data processing and quality control Quality control and data processing (Supplementary Note) were performed using our RaMWAS Bioconductor package, which is specifically designed for large-scale methylation studies. After rigorous quality control of samples, reads, and CpGs, 1132 subjects (320 controls and 812 cases) with an average of 48.7 million reads per sample (=81.9% of all reads) remained. For each of these individuals, our dataset included commonly methylated high-quality methylation information for 21,869,561 CpGs27. Among these, 970,414 were common CpG-SNPs (CpGs created/destroyed by SNPs with minor allele frequency > 10%) that were used for MWAS. To identify the CpG-SNPs we used directly genotyped and imputed genotype information (Supplement) from the NESDA participants. The imputed SNPs were filtered by imputation R2 ≥ 0.9 and minor allele frequency ≥ 0.1 in cases and controls. Finally, an in silico experiment described elsewhere28 was used to remove CpG-SNPs in loci showing alignment problems. MWAS of CpG-SNPs To test for association between the methylation level at each CpG-SNP and MDD, we performed a regression analysis with four sets of covariates. First, we regressed out 19 assay-related variables (i.e., potential technical artifacts) including the quantity of methylation-enriched DNA captured, sample batches, and peak location22. Second, we regressed out the demographic variables age and sex. Third, to avoid confounding due to cell-type heterogeneity, we regressed out blood cell type proportions as estimated by the methylation data29 using MBD-Seq “reference methylomes” we generated after isolating all common cell types in blood30. Fourth, principle component analysis was used to capture any remaining unmeasured source of variation. Specifically, using a scree test we selected the first principle component. The MWAS was performed by fitting the following regression equation:31 $$Y = b_0 + b_1{\mathrm{CpG}} - {\mathrm{SNP}} + b_2\left( {{\mathrm{CpG}} - {\mathrm{SNP}} \times {\mathrm{MDD}}} \right) + b_3{\mathrm{MDD}} + b_4X_1 + \ldots + b_kX_k + E,$$ where Y are the CpG scores, b0 is the intercept of the regression line, b4…bk the effects of covariates, and E are the residual effects. The CpG-SNP is coded as 0, 1, and 2, which corresponds to having 0, 1, or 2 copies of the SNP allele that creates/destroys the CpG relative to the reference genome. MDD is coded 0 for controls and 1 for cases. Figure 1 describes nine scenarios for how the regression lines change with alterations of the b2 and b3 parameters, when b1 is equal to 1. A non-zero value of parameter b1 indicates that the site is methylated with the amount of methylation being proportional to the number of CpGs (i.e., has a methylation quantitative locus or “meQTL effect”). Parameter b2 estimates the case-control difference at the CpG-SNP site that is proportional to the number of CpGs (i.e., the “CpG-SNP dose effect”). Parameter b3 captures case-control differences from nearby sites and thus do not depend on the number of CpG creating alleles of the SNP (i.e., a “local effect”). MBD-Seq assays the methylation of regions that are about the size of the sequenced fragments (~150 bp). Therefore, part of the differences observed at the CpG-SNP may reflect the effects of nearby CpGs resulting in non-zero values b3. In the overall association (i.e., “CpG-SNP MWAS”) we tested the null-hypothesis, H0: b2 = b3 = 0. Fig. 1: Overview of possible scenarios for the regression lines. Keeping the b0 and b1 parameters constant while altering the b2 and b3 parameters result in nine possible scenarios for the regression lines. The b0 parameter is kept at 0.2 and b1 (the meQTL effect) is equal to 1 for all plots. The value for the b2 parameter (the dose effect) was altered between a value equal to zero (0; left column), a positive value (0.5; middle column), and a negative value (−0.5; right column). Similarly, the value for the b3 parameter (the local effect) was altered between a value equal to zero (0; top row), a positive value (0.25; middle row), and a negative value (−0.25; bottom row). A non-zero value for b3 means that the locus is affected by a case-control differences from nearby CpGs (i.e., a “local effect”). This effects is independent of the effect from b2 and can either enhance or (partly) diminish the dose effect Full size image Permutation of CpG-SNP MWAS to study the null distribution To test if the lambda observed for the MDD CpG-SNP MWAS was caused by associations to the outcome variable, or if it was caused by that the test statistic distribution did not follow the theoretical null we used permutations. Using exactly the same dataset, we performed MWAS for 100 permutations of the MDD outcome variable and recorded the lambdas. Next, the observed association P-values from the MDD CpG-SNP MWAS were corrected for the average permutation-obtained lambda. Replication of cumulative MWAS signals by resampling To study the significance of the cumulative MWAS signals, we used the “ramwas7riskScoreCV” function in RaMWAS. Specifically, the function uses elastic nets32,33,34 as implemented in the R Glmnet package. Elastic nets are akin to multiple regression analysis but suitable for our scenario where the number of predictors is much larger than the number of observations. Elastic nets were fitted by setting the alpha parameter to zero (i.e., ridge regression that retains all predictive sites in the model). To avoid overfitting, k-fold cross-validation is used35. That is, the sample was randomly partitioned into k = 10 equal sized subsamples. Of the k subsamples, k−1 are used as a “training” set to fit the elastic net and obtain weights for each predictive methylation site. The estimated weights are then used in the remaining “test” set to predict the outcome from the methylation data. By alternating the subjects used in the training and test sets, predictions are obtained for all subjects in the study. RaMWAS repeats the entire cycle of CpG-SNP selection through MWAS followed by estimation of prediction weights using elastic nets for each of the k-folds. Because both the selection of CpG-SNPs and estimation of their weights are not affected by the participants in the test set, we obtain unbiased predictions of the outcome for each subject. Furthermore, the score of CpG-SNPs is for an important part determined by the number of CpGs. To capture only effects associated with MDD, we removed the effect of the number of CpGs from the methylation score prior to conducting the “in sample” replication. By testing whether these methylation predictions are significantly correlated with actual MDD status, we performed an “in sample” replication of the cumulative MWAS signal. Permutation-based enrichment test of overlap To perform enrichment tests of the overlap between datasets we used the “shiftR” R-package. shiftR first maps the two datasets to each other based on chromosomal location. In our analyses, no flanking regions were used. Thus, for SNPs we considered a single base position and for CpGs we considered two bases. Next, the P-values are used to cross-classify each mapped marker in the two datasets as being in the top or bottom. Based on the resulting 2 by 2 tables as input, shiftR tests the null hypothesis that the enrichment odds ratio equals 1. To perform these test, shiftR uses circular permutations36. Specifically, through fast bitwise operations, it shifts the mapping of the two datasets by a single random integer in each permutation. This approach to generate the empirical test statistic distribution under the null hypothesis preserves the correlational structure of the data. We used 1 million permutations for each test. Multiple thresholds can be specified to define “top findings” (i.e., for our analyses we used the top 1 and 5%). To account for this “multiple testing”, the same thresholds are used in the permutations where the test statistic distribution under the null hypothesis is generated from the most significant (combination of) thresholds. Three GWAS Three independent (meta-analysis of) GWASs were recently reported for MDD or related phenotypes. Similar to the phenotyping in the NESDA sample, the 23andMe study12 and the study by the Converge Consortium37 determined phenotype status using information about current or prior MDD diagnosis. In contrast, the GWAS meta-analysis performed by the Social Science Genetics Association Consortium (SSGAC)11 studied depressive symptoms, which for the majority of the individuals (>105,000 individuals out of 161,460) were assessed based on self-reported frequency an individual had experienced feelings of unenthusiasm/disinterest and depression/hopelessness during the past 2 weeks. Thus, this assessment was not a clinical diagnosis of depression nor a validated method for assessing depression symptoms. In contrast, when SSGAC studied neuroticism, an MDD-related phenotype, the status for the majority of individuals was assessed using a validated questionnaire that applied different harmonized neuroticism assessment batteries (n = 63,661) and a 12-item version of the Eysenck Personality Inventory Neuroticism38 (n = 107,245). Therefore, for the purpose of comparison with our MWAS for MDD, we used the SSGAC GWAS meta-analysis results of neuroticism11. For calculating the enrichment test statistic, shiftR classifies markers as being among the top vs. bottom results. However, from the 23andMe study, we could only get access to the P-values from the top 10,000 SNPs. To address this restriction we used SNPs retained in the multiple Psychiatric Genetic Consortia (http://www.med.unc.edu/pgc) studies after quality control. After removing the 10,000 top 23andMe SNPs, we assumed that these common and QC’ed SNPs were likely tested or were in LD with tested SNPs in the 23andMe study but yielded P-values lower than those of the top 10,000. The top 10,000 SNPs all had P-values < 10−5. To define a second threshold for the 23andMe study, we also selected the 745 SNPs with P-values < 10−8. To account for this “multiple testing”, the same two thresholds were used in the permutations. To maximize the compatibility of the analysis, all GWAS datasets were subjected to the same procedure as used for the 23andMe study. Results Methylome-wide CpG-SNP analysis We utilized the methylation data in combination with genotype information from the same individuals to perform a MWAS involving 970,414 common CpG-SNPs. Permutations of the MWAS generated an average lambda of 1.02 with a 95% confidence interval from 1.0087 to 1.0321. Thus, as shown in the Q-Q plot (Fig. 2a), the slightly inflated lambda (lambda = 1.062) observed for the MDD CpG-SNP MWAS is likely caused by a combination of true associations and by that the test statistic distribution did not follow the theoretical null distribution. As it would be practically non-feasible (too time-consuming) to obtain permutation P-values for each site we instead control for the deviation in the theoretical null distribution. Thus, the P-values were corrected for the average permutation-obtained lambda (Fig. 2b). Fig. 2: Q-Q plots and Manhattan plot of CpG-SNP MWAS. a Quantile-Quantile plot of the CpG-SNP MWAS before correction. The observed P-values, on a –log10 scale, are plotted against their expected values (gray main diagonal line) under the null hypothesis assuming none of the sites have an effect. Yellow lines indicate the 95% confidence intervals (CI). b Quantile-Quantile plot of the CpG-SNP MWAS after correction for permutation-obtained lambda. The deviation of P-values from the main diagonal indicates that, even after correction, there are potentially many markers associated with MDD. c Manhattan plot of the CpG-SNP MWAS. The plot shows the MWAS P-values on a –log10 scale (y-axis) by their chromosomal location (x-axis). The dashed line marks the threshold for suggestively significant findings (P = 1 × 10−5) Full size image The Manhattan plot (Fig. 2c) shows 27 suggestively significant loci (P < 1.00 × 10−5 after lambda correction) across the genome (Table 1). In Fig. S1, we show the regression plots for all the 27 sites. Twenty-five (92.6%) of the sites showed that the methylation levels were dependent on the number of CpG alleles (i.e., there was a significant meQTL effect) and 23 sites (85.2%) showed that this effect was different between cases and controls (i.e., there was a significant CpG-SNP dose effect). Thus, the associations observed for the two sites lacking CpG-SNP dose effects, as well as for the two sites that did not show significant meQTL effects, are likely caused by local effects from nearby CpGs. Table 1 CpG-SNP MWAS findings with P < 1.00e-5 Full size table Focusing on the 23 CpG-SNPs with both meQTL and CpG-SNP dose effects, we identified five sites (21.7%) with a positive dose effect. These sites showed a consistent pattern where the case-control difference gets bigger with more CpG-creating alleles but where the cases show higher methylation levels than the controls. The reaming 18 sites (78.3%) showed a negative dose effect. Thus, the negative dose effect occurred significantly more often (P = 0.0040) than expected by chance. A negative dose effect translates to (Fig. 1, right column) a case-control difference that gets bigger with more CpG-creating alleles and where the cases show lower methylation levels than the controls. As was shown in Fig. 1, the CpG-SNP MWAS associations (detected with MBD-Seq data) is in addition to a meQTL effect and a dose effect, also influenced by the local effect from nearby CpGs. This local effect can both enhance or diminish the CpG-SNP dose effect. The deviation of the observed P-values from the main diagonal, observed in the Q-Q plot (Fig. 2b) after correction for artificial inflation, suggests multiple sites are potentially associated with MDD. To study the significance of the cumulative CpG-SNP MWAS signal for large portions of the top markers, we used a resampling approach that fits elastic nets and employs k-fold cross-validation to avoid overfitting and obtain an unbiased estimate of the cumulative effect across markers. Results showed that the cumulative association was significant (P = 4.01 × 10−8), with the signal coming from the top 15,000 markers. Overlap between CpG-SNP MWAS and GWAS When comparing our MWAS results with three recent GWAS, a significant, or trend toward, enrichment was observed for all three GWASs when using the top 1% of results (the 1% threshold) for the CpG-SNP MWAS and the most stringent threshold for each of the three GWAS (Table 2). The highest enrichment was observed with the 23andMe study (P = 4.9 × 10−3, OR = 5.00) followed by SSGAC (P = 3.8 × 10−2, OR = 1.42) and Converge (P = 8.1 × 10−2, OR = 1.36). The overlap included 55 CpG-SNP sites (Table S2). The most significant site (P = 4.40 × 10−3) in the CpG-SNP MWAS that overlapped with the GWAS data was located in the Roundabout, axon guidance receptor, homolog 2 gene (ROBO2). The overlapping CpG-SNPs included 26 genes present in GO. These genes were overrepresented (P < 0.01) in 12 level-5 GO terms (Table 3). The most significant term (P = 4.57 × 10−4) was “Regulation of synapse organization” which, among other genes, included ROBO2. Table 2 Results from permutation-based enrichment tests of CpG-SNP MWAS and recent GWAS Full size table Table 3 Over-represented gene ontology terms (P < 0.01) among genes detected from overlapping CpG-SNP MWAS and GWAS findings Full size table Discussion Here we present the first MWAS of common CpG-SNPs (CpGs created/destroyed by SNPs with minor allele frequency > 10%) in MDD cases and controls. The methylation data were generated using a sequencing-based approach and involved 970,414 CpG-SNPs and 1132 individuals. Furthermore, we investigated the overlap of this study with recent GWAS for MDD, or related phenotypes. The MWAS suggested that multiple sites are potentially associated with MDD and resampling showed that the cumulative signal replicated. Furthermore, permutation-based enrichment tests suggested significant overlap with top findings from the MWAS and recent GWAS. Methylome-wide CpG-SNP analysis The majority of the associated CpG-SNPs that expressed a significant meQTL effect and a significant dose effect in the MWAS showed a distinct pattern where methylation increased with the number of CpG alleles present, but where this increase was attenuated in MDD cases compared to controls. Thus, cases often showed less methylation than controls at the differently methylated loci. Many possible explanation may exist for this pattern. However, consistent with a general function of DNA methylation that protects the integrity of the genome by inactivating DNA elements39,40, this pattern would be in agreement with that a portion of potentially damaging mutations might not be properly silenced in MDD cases. Interestingly, the same pattern with less methylation observed in cases than in controls has previously been observed also in CpG-SNP studies for psychosis using both blood and brain tissue31. Overlap between CpG-SNP MWAS and GWAS Many of the genes implicated by both the MWAS and the GWASs are of critical importance for neuronal function. Some of the overlapping gens have previously been associated with psychiatric disorders. For example, ROBO2 (roundabout, axon guidance receptor, homolog 2) is critical for the maintenance of inhibitory synapses in the adult ventral tegmental area, a brain region important for the production of dopamine41, and has been implicated in schizophrenia42,43,44 and bipolar depression45. ASIC2 (acid-sensing, proton-gated, ion channel 2) plays a role in neurotransmission46. DCC (deleted in colorectal carcinoma—netrin 1 receptor) upregulation in prefrontal cortex pyramidal neurons causes vulnerability to stress-induced social avoidance and anhedonia in mouse, and mutations in DCC have been associated with brain malformation47. Furthermore, DCC has been suggested to confer susceptibility to depression-like behaviors in mice and humans48 and was recently associated with mood instability, which has a strong genetic correlation to MDD49. In addition, the netrin 1 pathway, which involves DCC, has been identified as a candidate pathway for MDD50. Critically, both ROBO2 and DCC interact in opposing fashion and have strong roles in directing axon pathfinding in developing neurons51,52. In summary, several of the genes detected in the MWAS-GWAS overlap serve critical biological functions of likely relevance to MDD etiology. The overlap between the CpG-SNP MWAS and GWAS cannot be explained by the allele frequency differences between cases and controls that produce GWAS signals. It is true that methylation levels will be higher in the group with the higher frequency of the SNP allele that creates the CpG-SNP. However, these methylation differences are fully accounted for by the effect of the SNP as a “covariate” in the model we used for the CpG-SNP MWAS. Indeed, performing a GWAS with only the SNPs that were included in the CpG-SNP MWAS showed a lambda of 0.995 without any strong association signals (smallest P-value = 5.28 × 10−6). Thus, the CpG-SNP MWAS and GWAS provide additional and independent lines of evidence for the involvement of these loci in MDD. Conclusion In the first CpG-SNP MWAS for MDD, we identified 27 suggestively significant sites. A significant number of these sites showed a negative CpG-SNP dose effect with less methylation in cases than controls. Furthermore, the MWAS results were over-represented among findings from three recent GWASs, which for example added additional support for the involvement of DCC in MDD. As the analysis approach prevents the methylation results to be driven by allele frequency differences between cases and controls, these results show that MWAS and GWAS provide additional and independent lines of evidence for the involvement of these loci in MDD. In conclusion, CpG-SNP methylation studies of MDD can contribute novel and biologically relevant information that complements previous findings detected by regular MWAS or GWAS alone. Availability of data and materials Following local IRB approval individual level methylation data will be made available via dbGap (submission in preparation). Additional information Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References 1. Judd, L. L. The clinical course of unipolar major depressive disorders. Arch. General. Psychiatry 54, 989–991 (1997). Article Google Scholar2. Lopez, A. D., Mathers, C. D., Ezzati, M., Jamison, D. T. & Murray, C. J. Global and regional burden of disease and risk factors, 2001: systematic analysis of population health data. Lancet 367, 1747–1757 (2006). ArticlePubMed Google Scholar3. Wittchen, H. U., et al. The size and burden of mental disorders and other disorders of the brain in Europe 2010. Eur. Neuropsychopharmacol. 21, 655–679 (2011). ArticlePubMed Google Scholar4. Angst, F., Stassen, H. H., Clayton, P. J. & Angst, J. Mortality of patients with mood disorders: follow-up over 34-38 years. J. Affect. Disord. 68, 167–181 (2002). ArticlePubMed Google Scholar5. Kessler, R. C., et al. The epidemiology of major depressive disorder: results from the National Comorbidity Survey Replication (NCS-R). J. Am. Med. Assoc. 289, 3095–3105 (2003). Article Google Scholar6. Avenevoli, S., Swendsen, J., He, J. P., Burstein, M. & Merikangas, K. R. Major depression in the national comorbidity survey-adolescent supplement: prevalence, correlates, and treatment. J. Am. Acad. Child Adolesc. Psychiatry 54, 37–44 e32 (2015). ArticlePubMed Google Scholar7. Hardeveld, F., Spijker, J., De Graaf, R., Nolen, W. A., & Beekman, A. T. Prevalence and predictors of recurrence of major depressive disorder in the adult population. Acta Psychiatr. Scand. 122(3), 184–91 (2010). ArticlePubMed Google Scholar8. Mueller, T. I., et al. Recurrence after recovery from major depressive disorder during 15 years of observational follow-up. Am. J. Psychiatry 156, 1000–1006 (1999). PubMed Google Scholar9. Whiteford, H. A., Ferrari, A. J., Degenhardt, L., Feigin, V. & Vos, T. The global burden of mental, neurological and substance use disorders: an analysis from the Global Burden of Disease Study 2010. PLoS ONE 10, e0116820 (2015). ArticlePubMedPubMed Central Google Scholar10. World Health Organization. The Global Burden of Disease: 2004 Update (World Health Organization, Geneva, 2008). 11. Okbay, A., et al. Genetic variants associated with subjective well-being, depressive symptoms, and neuroticism identified through genome-wide analyses. Nat. Genet. 48, 624–633 (2016). ArticlePubMedPubMed Central Google Scholar12. Hyde, C. L., et al. Identification of 15 genetic loci associated with risk of major depression in individuals of European descent. Nat. Genet. 48, 1031–1036 (2016). ArticlePubMedPubMed Central Google Scholar13. Kaffman, A. & Meaney, M. J. Neurodevelopmental sequelae of postnatal maternal care in rodents: clinical and research implications of molecular insights. J. Child Psychol. Psychiatry 48, 224–244 (2007). ArticlePubMed Google Scholar14. Szyf, M., Weaver, I. C., Champagne, F. A., Diorio, J. & Meaney, M. J. Maternal programming of steroid receptor expression and phenotype through DNA methylation in the rat. Front. Neuroendocrinol. 26, 139–162 (2005). ArticlePubMed Google Scholar15. Abdolmaleky, H. M. et al. Methylomics in psychiatry: modulation of gene-environment interactions may be through DNA methylation. Am. J. Med. Genet. B Neuropsychiatr. Genet. 127B, 51–59 (2004). ArticlePubMed Google Scholar16. Kumar, D. et al. A functional SNP associated with atopic dermatitis controls cell type-specific methylation of the VSTM1 gene locus. Genome Med. 9, 18 (2017). ArticlePubMedPubMed Central Google Scholar17. Izzi, B. et al. Allele-specific DNA methylation reinforces PEAR1 enhancer activity. Blood 128, 1003–1012 (2016). ArticlePubMed Google Scholar18. Shilpi, A., Bi, Y., Jung, S., Patra, S. K. & Davuluri, R. V. Identification of genetic and epigenetic variants associated with breast cancer prognosis by integrative bioinformatics analysis. Cancer Inform. 16, 1–13 (2017). ArticlePubMedPubMed Central Google Scholar19. Onuma, H. et al. Dual effects of a RETN single nucleotide polymorphism (SNP) at -420 on plasma resistin: genotype and DNA methylation. J. Clin. Endocrinol. Metab. 102, 884–892 (2017). PubMed Google Scholar20. McClay, J. L. et al. High density methylation QTL analysis in human blood via next-generation sequencing of the methylated genomic DNA fraction. Genome Biol. 16, 291 (2015). ArticlePubMedPubMed Central Google Scholar21. Chan, R. F. et al. Enrichment methods provide a feasible approach to comprehensive and adequately powered investigations of the brain methylome. Nucleic Acids Res. 45, e97 (2017). ArticlePubMedPubMed Central Google Scholar22. Aberg, K. A. et al. A MBD-seq protocol for large-scale methylome-wide studies with (very) low amounts of DNA. Epigenetics 12, 743–750 (2017). ArticlePubMed Google Scholar23. Wittchen, H. U. Reliability and validity studies of the WHO-composite International Diagnostic Interview (CIDI): a critical review. J. Psychiatr. Res. 28, 57–84 (1994). ArticlePubMed Google Scholar24. Rush, A. J. et al. The Inventory for Depressive Symptomatology (IDS): preliminary findings. Psychiatry Res. 18, 65–87 (1986). ArticlePubMed Google Scholar25. Boomsma, D. I. et al. Genome-wide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects. Eur. J. Hum. Genet. 16, 335–342 (2008). ArticlePubMed Google Scholar26. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012). ArticlePubMedPubMed Central Google Scholar27. Shabalin, AA., Hattab, MW., Clark, SL., Chan, RF., Kumar, G., Aberg, KA., & van den Oord, EJCG. RaMWAS: fast methylome-wide association study pipeline for enrichment platforms. Bioinformatics. 34(13), 2283–2285 (2018). ArticlePubMed Google Scholar28. Aberg, K. A. et al. MBD-seq as a cost-effective approach for methylome-wide association studies: demonstration in 1500 case-control samples. Epigenomics 4, 605–621 (2012). ArticlePubMedPubMed Central Google Scholar29. Houseman, E. A. et al. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinformatics 13, 86 (2012). ArticlePubMedPubMed Central Google Scholar30. Hattab, M. W. et al. Correcting for cell-type effects in DNA methylation studies: reference-based method outperforms latent variable approaches in empirical studies. Genome Biol. 18, 24 (2017). ArticlePubMedPubMed Central Google Scholar31. van den Oord, E. J. et al. A whole methylome CpG-SNP association study of psychosis in blood and brain tissue. Schizophr. Bull. 42, 1018–1026 (2016). ArticlePubMed Google Scholar32. Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22 (2010). ArticlePubMedPubMed Central Google Scholar33. Simon, N., Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for Cox’s proportional hazards model via coordinate descent. J. Stat. Softw. 39, 1–13 (2011). ArticlePubMedPubMed Central Google Scholar34. Tibshirani, R. et al. Strong rules for discarding predictors in lasso-type problems. J. R. Stat. Soc. Ser. B Stat. Methodol. 74, 245–266 (2012). Article Google Scholar35. Hastie, T., Tibshirani, R., Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction (Springer Verlag, New York, 2001). 36. Cabrera, C. P. et al. Uncovering networks from genome-wide association studies via circular genomic permutation. G3 2, 1067–1075 (2012). ArticlePubMed Google Scholar37. Converge Consortium. Sparse whole-genome sequencing identifies two loci for major depressive disorder. Nature 523, 588–591 (2015). ArticlePubMed Central Google Scholar38. Eysenck, H. J., & Eysenck, S. B. G. (1975). Manual of the Eysenck Personality Questionnaire (Junior and Adult). Kent, UK: Hodder & Stoughton. 39. Chan, S. W., Henderson, I. R. & Jacobsen, S. E. Gardening the genome: DNA methylation in Arabidopsis thaliana. Nat. Rev. Genet. 6, 351–360 (2005). ArticlePubMed Google Scholar40. Kapoor, A., Agius, F. & Zhu, J. K. Preventing transcriptional gene silencing by active DNA demethylation. FEBS Lett. 579, 5889–5898 (2005). ArticlePubMed Google Scholar41. Gore, B. B. et al. Roundabout receptor 2 maintains inhibitory control of the adult midbrain.eLife. Apr 10(6), e23858 (2017). Google Scholar42. Potkin, S. G. et al. A genome-wide association study of schizophrenia using brain activation as a quantitative phenotype. Schizophr. Bull. 35, 96–108 (2009). ArticlePubMed Google Scholar43. Potkin, S. G. et al. Identifying gene regulatory networks in schizophrenia. Neuroimage 53, 839–847 (2010). ArticlePubMedPubMed Central Google Scholar44. Aberg, K. A., et al. A comprehensive family-based replication study of schizophrenia genes. JAMA Psychiatry 70, 573–581 (2013). ArticlePubMedPubMed Central Google Scholar45. Meda, S. A. et al. Multivariate analysis reveals genetic associations of the resting default mode network in psychotic bipolar disorder and schizophrenia. Proc. Natl Acad. Sci. USA 111, E2066–E2075 (2014). ArticlePubMed Google Scholar46. Kreple, C. J. et al. Acid-sensing ion channels contribute to synaptic transmission and inhibit cocaine-evoked plasticity. Nat. Neurosci. 17, 1083–1091 (2014). ArticlePubMedPubMed Central Google Scholar47. Marsh, A. P. et al. Mutations in DCC cause isolated agenesis of the corpus callosum with incomplete penetrance. Nat. Genet. 49, 511–514 (2017). ArticlePubMedPubMed Central Google Scholar48. Torres-Berrio, A. et al. DCC confers susceptibility to depression-like behaviors in humans and mice and is regulated by miR-218. Biol. Psychiatry 81, 306–315 (2017). ArticlePubMed Google Scholar49. Ward, J. et al. Genome-wide analysis in UK Biobank identifies four loci associated with mood instability and genetic correlation with major depressive disorder, anxiety disorder and schizophrenia. Transl. Psychiatry 7, 1264 (2017). ArticlePubMedPubMed Central Google Scholar50. Zeng, Y. et al. A combined pathway and regional heritability analysis indicates NETRIN1 pathway is associated with major depressive disorder. Biol. Psychiatry 81, 336–346 (2017). ArticlePubMedPubMed Central Google Scholar51. Tessier-Lavigne, M. & Goodman, C. S. The molecular biology of axon guidance. Science 274, 1123–1133 (1996). ArticlePubMed Google Scholar52. Zhang, C., Gao, J., Zhang, H., Sun, L. & Peng, G. Robo2–Slit and Dcc–Netrin1 coordinate neuron axonal pathfinding within the embryonic axon tracts. J. Neurosci. 32, 12589–12602 (2012). ArticlePubMed Google Scholar Download references Acknowledgements The current methylation study was supported by grant R01MH099110 from the National Institute of Mental Health. Funding for NESDA was obtained from the Netherlands Organization for Scientific Research (Geestkracht program grant 10-000-1002); the Center for Medical Systems Biology (CSMB, NWO Genomics), Biobanking and Biomolecular Resources Research Infrastructure (BBMRI-NL), VU University’s Institutes for Health and Care Research (EMGO+) and Neuroscience Campus Amsterdam, University Medical Center Groningen, Leiden University Medical Center, National Institutes of Health (NIH, R01D0042157-01A, MH081802, 1RC2MH089951, and 1RC2MH089995). Part of the genotyping and analyses were funded by the Genetic Association Information Network (GAIN) of the Foundation for the National Institutes of Health, for which computing was supported by BiG Grid, the Dutch e-Science Grid, which is financially supported by NWO. Author information AffiliationsCenter for Biomarker Research and Precision Medicine, Virginia Commonwealth University, Richmond, VA, USAKarolina A. Aberg, Andrey A. Shabalin, Robin F. Chan, Min Zhao, Gaurav Kumar, Shaunna L. Clark, Lin Y. Xie & Edwin J. C. G. van den OordDepartment of Psychiatry, Amsterdam Neuroscience, VU University Medical Center/GGZ inGeest, Amsterdam, The NetherlandsGerard van Grootheest, Yuri Milaneschi & Brenda W. J. H. Penninx AuthorsSearch for Karolina A. Aberg in:Nature Research journals • PubMed • Google Scholar Search for Andrey A. Shabalin in:Nature Research journals • PubMed • Google Scholar Search for Robin F. Chan in:Nature Research journals • PubMed • Google Scholar Search for Min Zhao in:Nature Research journals • PubMed • Google Scholar Search for Gaurav Kumar in:Nature Research journals • PubMed • Google Scholar Search for Gerard van Grootheest in:Nature Research journals • PubMed • Google Scholar Search for Shaunna L. Clark in:Nature Research journals • PubMed • Google Scholar Search for Lin Y. Xie in:Nature Research journals • PubMed • Google Scholar Search for Yuri Milaneschi in:Nature Research journals • PubMed • Google Scholar Search for Brenda W. J. H. Penninx in:Nature Research journals • PubMed • Google Scholar Search for Edwin J. C. G. van den Oord in:Nature Research journals • PubMed • Google Scholar Conflict of interest B.W.J.H.P. has received research funding (non-related) from Jansen Research and Boehringer Ingelheim. The remaining authors declare that they have no conflict of interest. Corresponding author Correspondence to Karolina A. Aberg. Electronic supplementary material Supplementary material Rights and permissions Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. About this article Publication history Received 15 December 2017 Revised 04 June 2018 Accepted 10 June 2018 Published 22 August 2018 DOI https://doi.org/10.1038/s41398-018-0205-8


This is a preview of a remote PDF: https://www.nature.com/articles/s41398-018-0205-8.pdf

Karolina A. Aberg, Andrey A. Shabalin, Robin F. Chan, Min Zhao, Gaurav Kumar, Gerard van Grootheest, Shaunna L. Clark, Lin Y. Xie, Yuri Milaneschi, Brenda W. J. H. Penninx, Edwin J. C. G. van den Oord. Convergence of evidence from a methylome-wide CpG-SNP association study and GWAS of major depressive disorder, Translational Psychiatry, 2018, DOI: 10.1038/s41398-018-0205-8