Methylome-wide change associated with response to electroconvulsive therapy in depressed patients

Electroconvulsive therapy (ECT) is a quick-acting and powerful antidepressant treatment considered to be effective in treating severe and pharmacotherapy-resistant forms of depression. Recent studies have suggested that epigenetic mechanisms can mediate treatment response and investigations about the relationship between the effects of ECT and DNA methylation have so far largely taken candidate approaches. In the present study, we examined the effects of ECT on the methylome associated with response in depressed patients (n = 34), testing for differentially methylated CpG sites before the first and after the last ECT treatment. We identified one differentially methylated CpG site associated with the effect of ECT response (defined as >50% decrease in Hamilton Depression Rating Scale score, HDRS), TNKS (q < 0.05; p = 7.15 × 10−8). When defining response continuously (ΔHDRS), the top suggestive differentially methylated CpG site was in FKBP5 (p = 3.94 × 10−7). Regional analyses identified two differentially methylated regions on chromosomes 8 (Šídák’s p = 0.0031) and 20 (Šídák’s p = 4.2 × 10−5) associated with ΔHDRS. Functional pathway analysis did not identify any significant pathways. A confirmatory look at candidates previously proposed to be involved in ECT mechanisms found CpG sites associated with response only at the nominally significant level (p < 0.05). Despite the limited sample size, the present study was able to identify epigenetic change associated with ECT response suggesting that this approach, especially when involving larger samples, has the potential to inform the study of mechanisms involved in ECT and severe and treatment-resistant depression.


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 response 1,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 depression 3,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 function 5,6 . There is evidence that pre-treatment methylation profiles may predict likelihood of achieving remission 7 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, see [7][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 forms 11 .
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 ECT 12 . Another candidate gene is brain-derived neurotrophic factor (BDNF) which is involved in neuroplastic changes associated with stress and depression 13,14 . Methylation differences of the promoter region of BDNF have been observed in many psychiatric disorders 15 , and also after use of antidepressant medication 16 . 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 rats 17 , and in humans, differences in expression of BDNF have been described 18,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 determined 20,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-userelated 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 patients 11 . Descriptive and clinical characteristics of the participants are described in Table 1.

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 fluctuations 23,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 in 27 , 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 to 28 . 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 crosshybridizing 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 in 29 .

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 package 30 . 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 field 31,32 and results from simulation experiments 33 .

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 .
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 in 11 , 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.

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  were nominally significant but did not remain significant after Šídák correction (see Table 3 and Fig. 1).
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.

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 trafficking 34 . Telomere length, a marker associated with aging is also known to be associated with psychiatric disorders 35 including major depressive disorder 36 , as well as depressive symptoms 37 . In several genome-wide association studies (GWASs) of depressionrelated 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 patients 43 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. RAP-GEF2, a protein-coding gene suggested to be involved in signal transmission, in BDNF receptor pathway signaling 45 , in schizophrenia 46 , and is found to be a target for regulated miRNAs in MDD 47 . PCM1, located in a chromosomal region on 8p, which has been implicated in various neuropsychiatric disorders including schizophrenia and depression 48 , is a protein-coding gene critical for cell division, and is involved in the proliferation and neurogenesis of neuroprecursors 49 .
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 depression 50 . FKBP5 demethylation resulting from childhood trauma has been linked to longterm stress hormone system deregulation and effects on immune function and brain areas associated with stress regulation 51 . Depressive phenotypes are shown to be associated with the age-related decrease in FKBP5 methylation 52 . 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 depression 53 . 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 remission 54 . 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 syndromerelated proteins are involved in the development of mental disorders 55 .
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 cancer 56,57 . In a recent large genome-wide geneenvironment analysis, its paralog, LRATD1, was observed to be associated with unipolar depression and response to trauma exposure 58 . 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 apoptosis 59 , and NNAT is involved in brain development and neuronal differentiation 60 . 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 action [61][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 stress 51,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 findings 65 . 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.