Pilot study of genome-wide DNA methylation and gene expression for treatment response to escitalopram in panic disorder

Abstract

BACKGROUND

Antidepressants, particularly selective serotonin reuptake inhibitors, are currently considered the first-line treatment for panic disorder (PD). However, little is known about the relationship between the biomarkers that may predict better treatment.

AIM

To compare genome-wide methylation and gene expression patterns between responsive and non-responsive patients with PD after 4 wk of escitalopram treatment.

METHODS

Thirty patients with PD were enrolled in this study (responders = 13; non-responders = 17). All patients were assessed using the PD Severity Scale-Chinese version before and after treatment. The Illumina Infinium MethylationEPIC (850k) BeadChip for genome-wide methylation screening and mRNA sequencing was used in all patients with PD.

RESULTS

A total of 701 differentially methylated positions (DMPs) were found between responders and non-responders (|Δβ| ≥ 0.06, q < 0.05), and the hyper- and hypomethylated CpG sites were 511 (72.9%) and 190 (27.1%), respectively. Relative to non-responders, there were 59 differential transcripts, of which 20 were downregulated and 39 were upregulated (q < 0.05). However, no differentially expressed genes were identified by mRNA sequencing after correcting for multiple testing (|log2(FC)| > 1, q > 0.05).

CONCLUSION

This preliminary study showed that DMPs might be associated with the treatment response to escitalopram in PD; however, these DMPs need to be verified in large samples.

Key Words: Panic disorder, Methylation, Expression profiling, Escitalopram

Core Tip: No genome-wide methylation studies or mRNA sequencing have been conducted to identify early response biomarkers in patients with panic disorder (PD). This study aimed to compare genome-wide methylation and gene expression patterns between responsive and non-responsive patients with PD after 4 wk of escitalopram treatment. A total of 701 differentially methylated positions (DMPs) were found between responders and non-responders, and the hyper- and hypomethylated CpG sites were 511 (72.9%) and 190 (27.1%), respectively. This preliminary study showed that DMPs might be associated with the treatment response to escitalopram in PD.

INTRODUCTION

Panic disorder (PD) is a common anxiety disorder characterized by recurrent and unexpected panic attacks. The estimated 12-month and lifetime prevalences of PD are 2.4% and 3.8%, respectively[1,2]. A cross-national epidemiological study reported that the lifetime prevalence of panic attacks was 13.2%[3]. Patients with PD experience symptoms such as tachycardia, chest pain, breathlessness, and dizziness. Consequently, they most frequently seek care in medical settings, such as emergency departments, and their condition is often undiagnosed[4]. PD causes substantial suffering and increases economic costs for both patients and society[5]. Hence, the availability of effective treatment may not only benefit individual patients but also provide economic returns to society.

Antidepressants, particularly selective serotonin reuptake inhibitors (SSRIs), are currently considered first-line treatments for PD. A meta-analysis suggested that SSRIs result in high remission rates with a minimal risk of adverse events in the treatment of PD[6,7]. However, not all patients benefit from the antidepressant therapy. For example, in a 12-week clinical trial of sertraline, citalopram, escitalopram, and paroxetine, 11.1%–18% of patients with PD did not respond to treatment[8,9]. Thus, predictive biomarkers of antidepressant response could greatly benefit clinical practice by decreasing the duration of drug efficacy evaluations.

Previous studies have found that common genetic variants could explain 42% of the individual differences in antidepressant response[10]. In addition, gene polymorphism studies focusing on 5-HTT, 5-HTR1A, 5-HTR2A, COMT, BDNF, and P450 (CYP) have been conducted in several large-scale studies on antidepressant drug responses in PD[11-14]. However, no consistent findings have been obtained from these studies. For example, with regard to treatment response to SSRIs in patients with PD, Zou et al[13] found that 5-HTTLPR polymorphism, rather than 5-HTR1A, may be an early predictor of response to sertraline in 2020, whereas Yevtushenko et al[11] indicated the importance of a 5-HT1A receptor gene polymorphism in 2010. Therefore, genetic variation in a single gene cannot fully explain individual differences in the response to treatment.

Emerging evidence from human and animal studies suggests a key role of epigenetic markers, including DNA methylation and histone modifications, in the prediction of antidepressant response[15]. Methylation of some candidate genes, such as SLC6A4, BDNF, and IL11, has shown promising results as a biomarker for predicting antidepressant responses in major depressive disorder (MDD). However, the research methods and results have been heterogeneous[16]. Genome-wide methylation analysis can accurately determine the location of DNA methylation in the genome and screen for effect-related differentially methylated genes. Currently, 850k methylation BeadChip is a new generation of DNA methylation chips developed based on the original 450k methylation chip[17]. The newly developed Human Methylation 850 BeadChip covers over 850k CpG methylation sites and is the most useful tool for analyzing the DNA methylation profile of the human genome. In addition, DNA methylation is a key epigenetic mechanism involved in the developmental regulation of gene expression[18]. Nevertheless, no genome-wide methylation studies or mRNA sequencing have been conducted to identify early response biomarkers in patients with PD.

Hence, a pilot study of genome-wide DNA methylation and gene expression analysis was conducted in patients with PD using the Illumina 850k BeadChip and mRNA sequencing to identify possible predictors of treatment response to escitalopram.

MATERIALS AND METHODS

Study sample

All patients with PD were recruited from the inpatient and outpatient departments of Sichuan Provincial People’s Hospital between March 2019 and December 2020. All patients satisfied the following inclusion criteria: (1) Primary diagnosis of PD in line with the Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition (DSM-IV) performed by a professional psychiatrist through a standardized structured clinical interview for DSM-IV Axis I disorders[19], with none of the patients having any other psychiatric disorders; (2) All PD patients were medication-naïve and had no history of any antidepressant or other psychotropic medication intake; (3) All subjects were without suicide ideations and attempts; (4) Subjects had no family history of psychiatric disorders in a first-degree relative; (5) All patients were Han Chinese aged 18 to 60 years; (6) All subjects were free of acute or chronic somatic disorders, head trauma, or neurological illnesses, and all subjects were free of alcohol consumption within 2 wk before their examination. Women were non-pregnant and non-nursing; and (7) All patients received escitalopram (10–20 mg qd) for 4 wk. The study was approved by the Sichuan Provincial People’s Hospital ethics committee [reference number: (2021) Ethics Review (313)]. All individuals provided written informed consent prior to the initiation of study procedures.

Clinical assessment

Patients were assessed at baseline and after 4 wk of treatment using the PD Severity Scale (PDSS). The PDSS comprises seven items, and participants are instructed to rate each item from 0 (none) to 4 (extremely severe) based on the severity of each symptom[20]. The PDSS-Chinese version has good internal consistency (Cronbach’s alpha) with an overall score of 0.83[21]. Treatment response for PD was defined as a reduction in the pretreatment PDSS score of at least 40% at 4 wk[22]. The PD samples were divided into responder and non-responder groups.

DNA and RNA extraction

A 4-mL ethylenediaminetetraacetic acid-anticoagulated peripheral blood sample was collected from all subjects in a fasting state via venipuncture between 7:00 a.m. and 8:30 a.m. DNA and RNA were extracted from whole blood samples from PD patients at baseline using DNeasy Blood and miRNeasy Mini kits (Qiagen, Hilden, Germany), respectively, according to the manufacturer’s protocol.

Genome-wide DNA methylation analysis

DNA was subjected to sodium bisulfite conversion using the EZ DNA Methylation Kit (Zymo Research, Irvine, CA, United States). The Illumina Infinium Human Methylation 850K BeadChip (Illumina, San Diego, CA, United States) was used to assess genome-wide DNA methylation, with a genome-wide coverage of over 850k CpG methylation sites per sample according to the manufacturer’s instructions. DNA quality control, bisulfite conversion, genome-wide methylation analysis, and initial methylation signal detection quality control were performed by Sinotech Genomics Co., Ltd. (Shanghai, China).

The raw intensity data were imported into R (v4.0.0; R Foundation for Statistical Computing, Vienna, Austria) and analyzed using the ChIP (ChAMP) package (v2.18.2) for data preprocessing, normalization, and comparison. β-values (ranging from 0 to 1) were used to determine the DNA methylation levels at each CpG site. Probes were filtered, including probes with a detection p-value of < 0.01, probes with less than three beads in at least 5% of the samples, non-CpG probes, and multi-hit probes. Beta-Mixture Quantile was used to normalize β-value matrices to adjust for type I and II probe biases[23]. In addition, we used singular value decomposition analysis to analyze the batch effect caused by the BeadChip Slide and Array and applied Combat to correct this batch effect[24]. Finally, all CpG sites were annotated using EPICanno software. ilm10b5. hg19.

Gene expression analysis

An Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) was used to assess RNA quality. A Qubit^® 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States) and NanoDrop One spectrophotometer (Thermo Fisher Scientific) were used to quantify total RNA, after which the library was constructed. The library for polymerase chain reaction product purification by paired-end libraries was synthesized using the mRNA-seq Lib Prep Kit for Illumina (ABclonal, Wuhan, China) following the sample preparation guide. Libraries were quantified using a Qubit^® 3.0 Fluorometer (Thermo Fisher Scientific) and validated using an Agilent 2100 Bioanalyzer (Agilent Technologies) to calculate the molar concentration and insert size. Clusters were generated using the cBot user guide and sequenced using an Illumina NovaSeq 6000 system (Illumina). Library construction and sequencing were performed by Sinotech Genomics Co. Ltd.

Hisat2 (2.1.0) was used to map the paired-end sequence files (fastq) to the reference genome (GRCh38.91). The output sequencing alignment/map files were converted into binary alignment/map files and sorted using SAM tools (v1.7). Gene abundance was expressed as fragments per kilobase of exons per million reads mapped (FPKM). StringTie software was used to count the fragments within each gene, and the TMM algorithm was used for normalization (https://www.ncbi.nlm.nih.gov/COG/). Finally, the FPKM value of each gene was calculated. EdgeR software was used to conduct differential expression analysis of the mRNA[25].

Statistical analysis

Data were analyzed using SPSS software (version 18.0; SPSS Inc., Chicago, IL, United States). Intergroup comparisons of continuous variables were conducted using Student’s t-test, and categorical variables were analyzed using Pearson’s chi-square test. Differentially methylated CpG positions were identified using ChAMP. Δβ was calculated as the difference of mean β-values between responders and non-responders. We set the criteria for differential methylation positions (DMPs) as the calling significance of an absolute change in Δβ-value between groups (|Δβ| ≥ 0.06 and q < 0.05). Differentially expressed mRNAs were defined as |log2(FC)| > 1 and q < 0.05. To identify genes with the same function and pathway in DNA methylation and gene expression profiling, we performed gene ontology (GO) (http://geneontology.org/) analysis for biological processes, cellular components, and molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://www.kegg.jp/) for differentially methylated and expressed genes. For all statistical analyses, a two-tailed p-value < 0.05 indicated statistical significance. The q-value cut-off of 0.05 was corrected using the Benjamini–Hochberg method for multiple hypotheses[26].

RESULTS

Demographic and clinical characteristics

PD samples were divided into 13 responders (6 men and 7 women; 33.54 ± 11.64 years old) and 17 non-responders (8 men and 9 women; 33.88 ±9.65 years old) at 4 wk. No statistically significant differences were found between the responders and non-responders in terms of sex (χ² = 0.002, p > 0.05) or age (t = -0.088, p > 0.05). There were 4 smokers (30.8%) and 9 nonsmokers (69.2%) in the responder group, and 5 smokers and 12 nonsmokers in the non-responder group. The number of drinkers and nondrinkers among the responders was 3 (23.1%) and 10 (76.9%) in the responder group, and 5 drinkers and 12 nondrinkers in the non-responder group. No statistically significant group difference was observed in terms of smoking and drinking (χ² = 0.006, 0.151; p > 0.05). In addition, the average PDSS total score was 13.08 ± 4.25 in responders and 14.94 ± 3.25 in non-responders at the start of treatment, and no statistically significant group difference was observed (t = -1.363, p > 0.05). The average PDSS total score for responders was 4.54 ± 2.11, which was significantly lower than that for non-responders 10.35 ± 2.00 after 4 wk of treatment (p < 0.001).

Differential methylation analysis

Compared with non-responders, there were 701 DMPs (|Δβ| ≥ 0.06, q < 0.05), and the hyper- and hypomethylated CpG sites were 511 (72.9%) and 190 (27.1%), respectively (Supplementary Table 1). The distribution of DMPs in the CpG islands was as follows: N-shore (34.6%), S-shore (24.5%), N-shelf (16.5%), island (14.4%), and S-shelf (10.1%). These DMPs in different regions of the gene were in the following order: Body region (64.4%), transcriptional start site (TSS) 1500 (15.3%), 5′-untranslated region (UTR) region (12.5%), 1^st Exon (3.4%), TSS 200 (3.1%), and the 3′-UTR region (1.3%). According to the University of California Santa Cruz annotation, 437 unique genes were identified, including 311 hypermethylated and 122 hypomethylated genes; the coexistence of hyper- and hypomethylation sites was found in four genes (Supplementary Table 2).

Enrichment analysis of differentially methylated genes

GO analysis revealed distinct functional categories for the associated genes, and approximately 226 biological processes, 14 cellular components, and 22 molecular functions were identified in differentially methylated genes (q < 0.05). Figure 1 shows the top 30 GO annotations for the differentially methylated genes. KEGG analysis showed that these differentially methylated genes were involved in 43 significant pathways (p < 0.05). However, no statistically significant pathways were found after adjusting using the Benjamini–Hochberg method (q > 0.05) (Supplementary Table 3).

Figure 1 Top 30 of gene ontology enrichment in differentially methylated genes. GO: Gene ontology.

Differential mRNA expression analysis

Compared with non-responders, there were 59 differential transcripts after adjusting using the Benjamini–Hochberg method, of which 20 were downregulated and 39 were upregulated (q < 0.05) (Figure 2). After the data were analyzed using the screening procedure, the results of mRNA sequencing showed that 132 differentially expressed genes were identified (p < 0.05). However, no statistically significant group differences between responders and non-responders were observed after adjusting using the Benjamini–Hochberg method (q > 0.05) (Supplementary Table 4).

Figure 2 Differential transcripts between responders and non-responders.

DISCUSSION

DNA methylation plays an important role not only in the diagnosis of diseases but also in the prediction of efficacy[16,27]. In this study, a pilot investigation of genome-wide DNA methylation in the early response to antidepressants in patients with PD was conducted using the Illumina 850k BeadChip. To our knowledge, this is the first genome-wide DNA methylation study of the treatment response to escitalopram in PD. According to the screening criteria for differential methylation sites, 701 DMPs were screened among responders and non-responders in patients with PD, and these DMPs were located within 437 unique genes. Similarly, one study assessed genome-wide DNA methylation using the Infinium MethylationEPIC BeadChip in patients with MDD for escitalopram treatment response. They identified 2571 significant DMPs, and 303 DMPs with an absolute change in Δβ-value between groups larger than 0.2[28]. These findings indicate that DMPs are potential peripheral predictors of antidepressant treatment response, and present an important opportunity to improve symptoms through prediction of medication response. Not only that, there is also some emerging evidence to suggest that PD patients have aberrant DNA methylations[29]. Psychiatry as a medical discipline, a diagnosis identifying a disorder should lead to an effective therapy[30]. Hence, epigenetic factors contributing to antidepressant response will be a unique and promising opportunity to implement personalized medicine in PD treatment. In addition, in the present study, most DMPs with an absolute change in value between groups were lower than 0.1. Furthermore, Ju et al[28] found that many DMPs showed very small differences in methylation (Δβ < 0.5%), making it difficult to fully explain efficacy prediction from a single CpG site. The treatment efficacy may be largely determined by the combined effects of multiple loci.

Notably, the current study identified 10 DMPs with an absolute change in Δβ-value larger than 0.2 (|Δβ| ≥ 0.2) between groups, and these genes included HLA-DPB1, HLA-DPA1, PDE1A, COL23A1, RUFY4, FRMD5, and SOHLH1. For example, the HLA-DPB1 site (cg12865025) is located in the TSS 1500 region. The major histocompatibility complex class II molecule consists of non-covalently associated alpha and beta chains, and this class II molecule is a heterodimer consisting of an alpha (DPA) and a beta chain (DPB)[31]. HLA-DPA1 and HLA-DPB1 are expressed on the surfaces of antigen-presenting cells; they play a central role in the adaptive immune system[32] and are associated with the pathophysiology of psychiatric disorders[17]. In recent years, increasing evidence has suggested that aberrations in immune-related pathways contribute to the pathophysiology and prediction of treatment responses in patients with PD[33,34]. Moreover, escitalopram led to a decrease in immune system activation[35], augmented tumor necrosis factor- peripheral secretion, induced faster kinetics of interleukin-1β secretion[36], and decreased Interleukin-17 levels[37]. Therefore, immune-related pathways may play an important role in the response to antidepressants in patients with PD. However, whether these DMPs can be used as molecular markers needs to be verified in larger samples.

In this study, functional enrichment analysis of differentially expressed genes was performed. GO analysis identified 262 GO terms in differentially methylated genes. Notably, the top 30 GO annotations, such as neuron recognition, neuron maturation, and cerebral cortex cell migration, may be associated with the treatment response. These enriched genes, including DAB1, NRP1, DIXDC1, ROBO1, PEX5, APP, ANKS1A, EDNRB, BCL2, and NTM, play central roles in brain development, neuronal remodeling, neurite outgrowth, and adhesion. For example, antidepressants such as citalopram increase APP secretion in primary rat neuronal cultures[38]. In addition, previous findings suggest that citalopram reduces mutant APP, Aβ, and mitochondrial toxicities and may have a protective role against mutant APP and Aβ-induced injuries in patients with depression, anxiety, and Alzheimer’s disease[39]. In addition, these genes have a potential role in the development of many neuropsychiatric diseases. For example, a previous study found that Dab1 knockout mice exhibited behavioral abnormalities, including hyperactivity, decreased anxiety-like behavior, and impaired working memory[40]. Hence, these differentially methylated genes in the GO terms may be involved in the pathophysiology and therapeutic response in PD, and the functions of these genes need further exploration. Further studies are needed to explore whether differential methylation in peripheral blood is similar to that in the brain.

Few studies have explored the association between gene expression and the prediction of therapeutic response in patients with PD. In this study, mRNA sequencing was compared between responders and non-responders among patients with PD. Unfortunately, no statistically significant group differences were found in differentially expressed genes. Furthermore, no candidate genes were identified based on integrative analysis of differential DNA methylation and expression. This may be due to the following reasons. First, the sample size was small. The different clinical features of patients with PD are closely related to the heterogeneity of heredity. Second, the use of whole blood, which is a mixed cell type sample, may limit the identification of mRNA expression changes. Finally, gene expression is affected by many factors, such as environmental and other epigenetic factors, especially non-coding RNA, which represents a promising source of peripheral biomarkers of antidepressant responses[41]. Ultimately, it is likely that combinations and not individual biomarkers will have the greatest utility in predicting antidepressant responses. In the future, researchers should integrate multiple types of information, including genetic, epigenetic, and gene expression data, to identify the most meaningful panel of biomarkers.

CONCLUSION

In conclusion, this preliminary study showed that DMPs may be associated with an early response to antidepressants in patients with PD. However, the results of our study should be considered in light of the following limitations. Since this is the first study to investigate genome-wide DNA methylation and gene expression in the prediction of early response to antidepressants in PD, it would be valuable to replicate our findings in a larger cohort. Furthermore, these DMPs must be verified using larger sample sizes. Lastly, it is necessary to test the association between treatment outcomes and changes in DNA methylation and gene expression from pre- to post-treatment with a long-term follow-up.

ARTICLE HIGHLIGHTS

Research background

Selective serotonin reuptake inhibitors are currently considered the first-line treatment for panic disorder (PD). However, not all patients benefit from the antidepressant therapy.

Research motivation

No genome-wide methylation studies or mRNA sequencing have been conducted to identify early response biomarkers in patients with PD.

Research objectives

To compare genome-wide methylation and gene expression patterns between responsive and non-responsive patients with PD after 4 wk of escitalopram treatment.

Research methods

Thirty patients with PD were enrolled in this study (responders = 13; non-responders = 17). All patients were assessed using the PD Severity Scale-Chinese version before and after treatment. The Illumina Infinium MethylationEPIC (850k) BeadChip for genome-wide methylation screening and mRNA sequencing was used in all patients with PD.

Research results

A total of 701 differentially methylated positions (DMPs) were found between responders and non-responders (|Δβ| ≥ 0.06, q < 0.05), and the hyper- and hypomethylated CpG sites were 511 (72.9%) and 190 (27.1%), respectively. Relative to non-responders, there were 59 differential transcripts, of which 20 were downregulated and 39 were upregulated (q < 0.05).

Research conclusions

This preliminary study showed that DMPs might be associated with the treatment response to escitalopram in PD, however, these DMPs need to be verified in large samples.

Research perspectives

DNA methylation contributing to antidepressant response will be a unique and promising opportunity to implement personalized medicine in PD treatment.

ACKNOWLEDGEMENTS

The authors would like to thank all participants in this study.