Introduction

Depression is a major contributor to global burden of disease and despite worldwide research efforts, the heterogeneous nature of the disorder makes it difficult to definitively unravel its underlying etiology and the factors influencing treatment response1,2. Electroconvulsive therapy (ECT) is an intervention with rapid and striking antidepressant effects and is the treatment of choice for patients with severe and treatment-resistant depression3,4. Studying biological changes associated with response to ECT in such a subgroup of severely depressed patients is a promising approach to gain insights into the underlying mechanisms of depression and treatment response.

DNA methylation is thought to be involved in disease pathology through its influence on gene expression and cellular function5,6. There is evidence that pre-treatment methylation profiles may predict likelihood of achieving remission7 and the literature on DNA methylation in depression has pointed to several possible genes of interest (e.g., BDNF, SLC6A4, NR3C1, FKBP5, etc. for reviews, see7,8,9,10) as being related to antidepressant response.

Investigating change of methylation levels during treatment may inform the biological processes underlying both depression and antidepressant response. Examining these changes in ECT patients offers an optimal research setting as: (1) treatment effects are substantial and occur soon after the intervention, and (2) ECT patients represent a subgroup of patients with the most severe form of depression. It is likely that this subgroup is not only clinically but also genetically more homogenous, especially as these patients tend to show a higher genetic burden for major depression than those with less severe forms11.

The few studies exploring differences in methylation related to ECT are only beginning to give insight into the factors involved. A translational study observed that the response to electroconvulsive stimulation (ECS) in rats was predicted by higher p11 promoter methylation and found this also to be the case in two human samples (n = 11; n = 65) in response to ECT12. Another candidate gene is brain-derived neurotrophic factor (BDNF) which is involved in neuroplastic changes associated with stress and depression13,14. Methylation differences of the promoter region of BDNF have been observed in many psychiatric disorders15, and also after use of antidepressant medication16. Studies in animals and humans also propose an involvement of BDNF in response to ECT: after ECS/ECT, alterations of BDNF levels have been found in rats17, and in humans, differences in expression of BDNF have been described18,19. Furthermore, associations between pre-treatment BDNF levels and ECT outcome have been described, but whether BDNF levels are indicative of response remains to be determined20,21.

Research on ECT-related changes in methylation levels has largely focused on candidate genes. However, as ECT is a non-specific treatment with widespread, yet uncertain effects on biological processes, the selection of predefined candidates, relying on prior knowledge may not give the full picture. As such, an approach investigating the whole methylome is expected to yield new and relevant insights. Few studies have made use of these technologies in ECT samples to date. Moschny et al.22 examined longitudinal DNA methylation before and after ECT in a small group of patients (8 responders and 4 non-responders). They did not find any significant differences in global DNA methylation, but identified eight genes potentially implicated in either ECT response or its mechanism through analysis of single probe variance, and two genes whose methylation changed during treatment course.

In the present study, we aimed to identify changes in methylation levels associated with the effects of ECT and to find potential biomarkers for antidepressant response. We obtained and compared epigenome-wide DNA methylation levels of ECT patients (n = 34) before and after ECT. Differentially methylated CpG sites and regions associated with response were examined. Pathway analyses were employed to search for functional pathways affected by ECT. Finally, we took a targeted look at methylation in genes which have been previously implicated in ECT response and depression-related studies.

Materials and methods

This study was approved by the Ethics Committee (II), Medical Faculty Mannheim, University of Heidelberg. All patients provided written consent after a detailed explanation of the content of the study. All experimental procedures were performed in accordance with the Declaration of Helsinki.

Participants

Patients (n = 34) diagnosed with a major unipolar depressive episode (based on International Classification of Diseases version 10, ICD-10), older than 18 years, and assigned to ECT treatment were recruited at the Department of Psychiatry and Psychotherapy at the Central Institute of Mental Health (CIMH) in Mannheim, Germany, between 2014 and 2016. Criteria for assignment to ECT were treatment-resistant depression (i.e., failure of two adequate dose-duration antidepressants from different classes, or psychotherapy in the current episode), positive experience of ECT from a previous episode, or severe depression with (a) psychotic symptoms, (b) severe suicidality, or (c) the refusal of food or fluid intake. Patients were excluded if they had any substance-use-related disorders (other than tobacco and alcohol use disorders) or a lifetime diagnosis of schizophrenia. All participants were of Caucasian descent. The patients kept the same medication regimens throughout ECT treatment. This sample is a subset of a sample reported in a previous genetic study in ECT patients11. Descriptive and clinical characteristics of the participants are described in Table 1.

Table 1 Patient demographics.
Full size table

Depression and DNA methylation assessment

The 21-item version of the Hamilton Depression Rating Scale (HDRS) was administered and blood was collected from patients prior to the first (T1) and after the last (T2) session of ECT (average sessions: 10.82, range: 5–25 sessions). T2 was between 1 and 7 days after the last ECT session. Time of collection was kept constant between 8:30 and 9:30 a.m. at both timepoints to keep it close to the clinical interview conducted around the same time, as well as preclude any potential confounding effects arising from ward routines and circadian fluctuations23,24. Standard procedures for extraction and processing were followed. DNA extraction was performed using the Chemagic Magnetic Separation Module 1 (Chemagen Biopolymer-Technologie AG; Baesweiler, Germany). All genomic DNA samples were stored at −20 °C prior to analysis. Epigenome-wide DNA methylation was measured using the Illumina Infinium Methylation EPIC array (>850,000 CpG sites). The arrays were processed at the Genome Analysis Center at the HelmholtzZentrum München and Max-Planck-Institute for Psychiatry in Munich, Germany.

ECT treatment

ECT was conducted with a Thymatron IV device (Somatics, LLC. Lake Bluff, IL, USA). Anesthetic drugs used were: s-ketamine (~1.0 mg/kg)25,26 and succinylcholine for muscle relaxation (~1.0 mg/kg). Seizure thresholds were titrated at the initial session; dosing in subsequent sessions was given at >2.5 above this threshold. If patients did not show a clinical improvement or seizures were insufficient, energy used was increased. Patients received 2–3 ECT sessions per week. All patients started with unilateral stimulation with the possibility to change to bilateral stimulation at the discretion of the ECT supervisor.

Statistical analysis

Data processing, quality control (QC) and other statistical analyses were performed using R (versions 3.4.4 and 3.6.3) analysis software (https://cran.r-project.org/).

Data preprocessing, QC, and filtering

Methylation values were extracted using an updated version of the pipeline indicated in27, adapted in-house for use with the Illumina Infinium EPIC array. Illumina background correction was applied to all intensity values. A detection p-value threshold of p < 10−16 was used and intensity values with detection p ≥ 10−16 were designated as missing data. The proportion of missing data points was determined, allowing the calculation of sample and CpG site-specific call rates. Samples with insufficient DNA quality as denoted by a call rate of <95% were excluded. Intensity values were quantile normalized for each of the six probe types present on the array separately. Intensity values were converted to methylation beta values according to the manufacturer’s recommendation. White blood cell fractions were estimated according to28. Five of the six resulting estimates were subsequently included as covariates in downstream analyses to control for influences of cell type distribution on DNA methylation. The estimate for granulocytes showed the highest variance inflation factor (VIF) and was omitted to avoid collinearity issues. CpG sites were filtered by removing cross-hybridizing probes, probes with high missing rate (>0.02), and probes linked to X- and Y-chromosomes. Correction for batch effects and other technical parameters was done by performing a principal component analysis on internal control probe intensity values and including the first 10 extracted principal components (PCs) as covariates in the downstream analyses. Prior to analysis, all methylation beta values were logit transformed (base 2) to M-values, which were used in downstream analysis as recommended in29.

Differentially methylated single CpG sites

Association testing of methylation M-values for each CpG site was done using a mixed linear model approach as implemented in the Limma R package. Participant ID was used as a blocking factor, and estimated cell fractions, 10 control probe PCs, age, sex, and smoking status were included as additional covariates to adjust for confounding factors. It was observed that sex was highly correlated with the 5th control probe PC, thus this PC was not included in the statistical models to avoid collinearity issues.

The main effects of interest specified included: (1) response (responders vs. non-responders), (2) timepoint (change after intervention), and (3) the interaction between timepoint and response (the difference in change between response groups). Models were calculated with response specified both as a binary (>50% decrease in HDRS score) and continuous variable (i.e., change in HDRS score, ΔHDRS).

In addition, we examined the relationship between baseline (T1) methylation and response (both binary and continuous) in additional models.

Significance was defined as false discovery rate (FDR) q < 0.05. Results at a suggestive threshold of p < 10−5 were also reported.

Differentially methylated regions

Differentially methylated region (DMR) analysis was performed on the results of the above analyses using the comb-p package30. Comb-p parameters were specified as: seed p-value = 0.001 and a maximum distance between probes of 500 base pairs. These parameters follow those used in previous studies in the field31,32 and results from simulation experiments33.

Pathway analysis

A Gene Ontology (GO) enrichment analysis was performed on the results of the different models using the missMethyl (v1.12.0) R package. We examined CpG sites at the suggestive threshold of 1 × 10−5.

Targeted examination of methylation change in candidate CpG sites

Change in methylation of candidate genes from the literature was examined in an exploratory search. First, we selected: (1) candidates implicated in reviews of DNA methylation and antidepressant medication, i.e., BDNF, MAOA, SLC6A2, SLC6A4, HTR1A, HTR1B, IL6, IL117; SLC6A4, NR3C1, FKBP5, and OXTR8; and (2) candidates specific to ECT, i.e., S100A10 (p11)12, RNF175, RNF213, TBC1D14, TMC5, WSCD1, AC018685.2, AC098617.1, CLCN3P1, AQP10, and TRERF122.

Autosomal CpG sites which were annotated to these candidate genes in the UCSC Genome Browser NCBI curated RefSeq (retrieved: August 10, 2018) were extracted from the results of the single site analyses above (for each variable of interest in both binary and continuous models). Furthermore, to examine their predictive value in our sample, association between baseline methylation and response was also examined. A secondary FDR correction was applied to the list of all candidate CpG sites to control for false positives.

Results

Descriptive statistics of the sample are shown in Table 1. In the sample analyzed, using binary criteria, (defined as decrease of HDRS score of more than 50%) 25 were responders and 9 were non-responders to ECT. Levene’s test found no significant differences between group variances. Briefly, as also reported for the total sample in11, binary response to ECT was positively correlated with sex (being male) while continuous response (ΔHDRS score) was also associated with male sex and positively correlated with increased age.

Single CpG site analysis

Binary response

In the binary response model, one significantly differentially methylated site (q < 0.05), cg10005358, mapped to TNKS, was observed as an effect of response. Eight sites reached a suggestive threshold of p < 1 × 10−5 (see Table 2). No significantly differentially methylated CpG sites were observed for the effect of timepoint (at p < 1 × 10−5; 22 CpG sites), or interaction effect (at p < 1 × 10−5; 12 CpG sites). Several CpG sites annotated to the same gene appeared among the top hits of these effects of interest (see Tables 2 and S1.1–S1.3, e.g., TNKS, PCM1, RAPGEF2, RAB21; all suggestive at p < 1 × 10−5).

Table 2 Top 10 differentially methylated CpG sites associated with binary response and ΔHDRS.
Full size table

Continuous response (ΔHDRS)

In the continuous response model, no effects yielded significantly differentially methylated CpG sites at q < 0.05. At a suggestive threshold of p < 1 × 10−5, 7, 9, and 5 differentially methylated CpG sites were observed for the effects of ΔHDRS, timepoint, and ΔHDRS × timepoint, respectively (see Table 2 and Tables S2.1–S2.3). FKBP5 (CHR 6: cg01294490) was the top hit for both effect of ΔHDRS (p = 4.46 × 10−7) and effect of interaction of ΔHDRS × timepoint (p < 3.94 × 10−7). FXR2 (CHR 17: cg02936535) was also observed among the top hits for all effects of interest (ΔHDRS p = 3.79 × 10−6; timepoint p = 6.78 × 10−6; ΔHDRS × timepoint p = 7.38 × 10−6).

Baseline methylation and response

Methylation at baseline was not significantly associated with either binary or continuous response. At a suggestive threshold of p < 1 × 10−5 baseline methylation was associated with binary response at 9 CpG sites and with ΔHDRS at 6 CpG sites (see Tables S3.1 and S3.2).

Differentially methylated region analysis

Two DMRs were identified as associated with effect of ΔHDRS in the continuous response model. One significant DMR on chromosome 8 was identified (3 probes, Šídák’s corrected p = 0.0031) and another on chromosome 20 (13 probes, Šídák’s corrected p = 4.2 × 10−5). The DMR on chromosome 8 (CHR 8: 127568854-127569023) is located in the LRATD2 (FAM84B) gene, while the chromosome 20 DMR (CHR 20: 36148620-36148861) is located in the BLCAP gene and in the promoter region of NNAT (791 base pairs upstream of the transcription start site, TSS). Two other regions on chromosomes 14 and 19 were nominally significant but did not remain significant after Šídák correction (see Table 3 and Fig. 1).

Table 3 Differentially methylated regions associated with ΔHDRS.
Full size table
Fig. 1: Manhattan plot of differentially methylated regions for ΔHDRS.
figure 1

Significantly differentially methylated regions associated with the effect of ΔHDRS are indicated in red (chromosomes 8 and 20).

Full size image

Analysis of the results from timepoint and timepoint × ΔHDRS models, as well as all binary response models, did not identify significant DMRs.

Pathway analysis

No significant pathways were observed in any of the models.

Candidate analysis

In the binary model, at nominal significance (p < 0.05 uncorrected), 43, 37, and 36 CpG sites were associated with response, timepoint, and response × timepoint interaction, respectively. The top 3 candidate CpG sites associated with binary response were: FKBP5 (cg01294490, p = 8.74 × 10−5), BDNF-AS (cg02386995, p = 1.51 × 10−4), and NR3C1 (cg23273257, p = 0.0012), but no CpG site tested survived secondary FDR correction for multiple testing.

In the continuous response model, 41, 43, and 42 CpG sites were nominally associated with ΔHDRS, timepoint and ΔHDRS × timepoint interaction, respectively. The top 3 candidate CpG sites associated with continuous response were, FKBP5 (cg01294490, p = 4.46 × 10−7), BDNF (cg15710245, p = 0.0033), and BDNF-AS (cg09878183, p = 0.0034). The site annotated to FKBP5 survived the secondary correction for multiple testing (q = 0.0004).

Discussion

The present study examined treatment-associated changes of DNA methylation levels in 34 patients in an epigenome-wide manner. By investigating the relationship between response status and change in methylation levels, this study identified several potential CpG sites involved in ECT response and outlines potential differences between response groups.

The top CpG site associated with binary response is located in TNKS, which is a protein-coding gene associated with blood pressure, alcohol consumption, implicated in cancer pathology, and involved in various processes such as the Wnt signaling pathway, telomere length, and vesicle trafficking34. Telomere length, a marker associated with aging is also known to be associated with psychiatric disorders35 including major depressive disorder36, as well as depressive symptoms37. In several genome-wide association studies (GWASs) of depression-related traits, TNKS was found to be associated with (p = 7.68 × 10−10)38, bipolar disorder (p = 3 × 10−6)39, and positive affect (p = 0.0003)40.

Among the top 5 CpG sites suggestively associated with binary response, genes associated with processes such as cell adhesion, cell growth, apoptosis in malignant tumors, protein metabolism (RAB21)41, and signaling in glucose metabolism (AQR)42 were found. RAB21 was previously reported to be associated with remission (p = 0.0103) in a GWAS of selective serotonin reuptake inhibitors in MDD patients43 and AQR as related to neuroticism (p = 9.58 × 10−8) and worry (p = 2.06 × 10−6)44, a well-known symptom in depression. RAPGEF2, a protein-coding gene suggested to be involved in signal transmission, in BDNF receptor pathway signaling45, in schizophrenia46, and is found to be a target for regulated miRNAs in MDD47. PCM1, located in a chromosomal region on 8p, which has been implicated in various neuropsychiatric disorders including schizophrenia and depression48, is a protein-coding gene critical for cell division, and is involved in the proliferation and neurogenesis of neuroprecursors49.

The top site in the continuous response analysis was located in FKBP5, a gene which is known to be an important endogenous regulator of the stress hormone system possibly linked to stress-related psychiatric disorders such as depression50. FKBP5 demethylation resulting from childhood trauma has been linked to long-term stress hormone system deregulation and effects on immune function and brain areas associated with stress regulation51. Depressive phenotypes are shown to be associated with the age-related decrease in FKBP5 methylation52. Altered epigenetic and genetic FKBP5 regulation may contribute to stress-related disease risk. Findings related to FKBP5 have pointed to it as important in the interaction with environment in stress-related disorders such as major depression53. The present study found an association between methylation in a CpG site in FKBP5 and the reduction in patients’ HDRS scores. Although the CpG site identified here has not been associated with antidepressant response in previous studies, the findings support FKBP5 as an important gene requiring further investigation in the present context.

Among the top 5 CpG sites suggestively associated with continuous response were CpG sites annotated to D2HGDH and FXR2. D2HGDH encodes for the enzyme D-2-hydroxyglutarate dehydrogenase, and is suggested to be downregulated in depressed patients during remission54. Proteins of the FXR family have commonly been reported in autism spectrum disorders, and evidence from GWASs in mood disorders and schizophrenia suggests that Fragile X mental retardation syndrome-related proteins are involved in the development of mental disorders55.

Two DMRs associated with continuous response were identified. The DMR on chromosome 8 lies in LRATD2 (FAM84B), which is known to be involved in gastric and prostate cancer56,57. In a recent large genome-wide gene-environment analysis, its paralog, LRATD1, was observed to be associated with unipolar depression and response to trauma exposure58. The DMR on chromosome 20 is located in BLCAP and in the promoter region of NNAT. BLCAP encodes a protein that regulates cell proliferation and reduces cell growth by stimulating apoptosis59, and NNAT is involved in brain development and neuronal differentiation60. Together with the single site results, these findings are in line with previous works in the field of ECT; alterations in mechanisms such as neurogenesis and neuroinflammatory immune response are proposed to be among the mechanisms of ECT action61,62,63.

Several candidate genes proposed in the literature were found harboring CpG sites with nominally significant changes between T1 and T2, and T1 methylation values for a number of them were also nominally associated with response (both binary and continuous) (see Tables S4.1–S6.2); the roles they play remain unclear. The present results appear to lend support to previous research which has suggested the importance of these candidates but these results are preliminary and further investigation is warranted. Also, the identification of FKBP5 in the present study suggests that future research should assess and control for factors such as childhood trauma and stress51,64.

This study had several limitations. Although the largest study to date, the present sample size was limited, and it is expected that future studies using a similar approach in larger samples will be able to further clarify our results. Sample size notwithstanding, we identified a single significantly differentially methylated CpG site, as well as some suggestive ones which need further investigation. While we assessed methylation levels in whole blood, ECT is applied to the brain; both central and peripheral mechanisms may be affected by the global nature of the treatment and care should be taken with the interpretation of these findings65. The possible effect of anesthesia and pharmacotherapy is a potential confounding factor in methylation studies. However, medication in each patient was kept constant during the ECT course, and there were no differences between patients regarding anesthesia administration or treatment dosage. Therefore, the observed changes are unlikely to have resulted from these medications.

The genes implicated in our findings have been previously involved in the etiology of depression and treatment response, but confirmation in larger samples is needed. Multi-center approaches and collaborative efforts could help in obtaining the sample sizes required to allow a more robust characterization of ECT response and give insights into the biological processes underlying the striking antidepressant effects of ECT.