Integrative DNA Methylation and Gene Expression Analysis of Cognitive Behavioral Therapy Response in Children and Adolescents with Obsessive-Compulsive Disorder; a Pilot Study

Introduction

Obsessive-compulsive disorder (OCD) is a serious and common mental disorder with an estimated prevalence of 1% to 3% in children and adolescents.¹ Pediatric OCD is associated with significant distress and marked interpersonal, academic, and occupational impairments,² which can have a detrimental impact on psychosocial development.³

Treatments for OCD include pharmacological approaches using selective serotonin reuptake inhibitors (SSRI) and psychological approaches involving cognitive behavioral therapy (CBT). The recommended first-line treatment for pediatric OCD is CBT. It has been proven to be more effective than no intervention and while comparable to treatment with SSRIs, CBT has a lower risk-to-benefit ratio compared to medication and a higher acceptability among patients and their families.⁴ However, there is significant variability in how children and adolescents with OCD respond to CBT, with 39% of patients showing adequate remission of their symptoms.⁵ Similar variability is observed when patients with early-onset OCD are treated with SSRI monotherapy (22% remission rates) or when CBT is combined with an SSRI (54% remission rates).⁶

Clinical guidelines recommend CBT as the first-line treatment for patients with mild to moderate symptoms. It can be combined with an SSRI as the initial treatment in more severe cases or when there is no adequate response to CBT alone.^7,8 Although the severity of OCD symptoms may help in guiding treatment selection, the observed variability in remission rates highlights the importance of identifying moderators and predictors of response to help clinicians optimize the initial treatment for a particular patient.

Several factors have been proposed as predictors of a poorer outcome to CBT in pediatric OCD, such as an older age, the severity of symptoms and impairment, co-morbidities and a family history of OCD.⁹ However, their importance and validity as predictors remain controversial. Genetic variants represent a potential source of predictors, with the study of such variants referred to as therapy genetics. The first evidence of the contribution of genetic variants to psychological therapy response came from candidate gene studies; however, these findings have proven to be difficult to replicate.^10,11 Recently, genome-wide association studies (GWAS) on outcomes following psychological therapy were published for both children and adults with anxiety disorders.^12,13 However, these studies were underpowered to detect the small effect size of single genetic variants with genome-wide significance.

Several investigations have explored the epigenetic process of DNA methylation and differential gene expression in treatment response. Early studies using candidate gene (BDNF, NGF, FKBP5, and MAO-A) approaches have demonstrated that changes in DNA methylation may underlie response to psychological therapies in a variety of disorders including OCD.^14–21 A small number of studies have examined the role of gene expression and the response to psychological therapy: two studies using FKBP5 as a candidate gene in post-traumatic stress disorder and two studies using genome-wide expression analysis in anxiety disorders.^12,22–24

Here, we performed a genome-wide methylation analysis using peripheral blood obtained after eight weeks of CBT from a cohort of children and adolescents with a diagnosis of OCD who were drug-naïve and never previously treated with psychological interventions. Furthermore, we integrated the methylation and gene expression profiles using samples from the same individuals. The main objective of the present study was to provide new insight into the biological mechanisms of CBT and to identify new candidate biomarkers of CBT response.

Materials and Methods

Subjects

Twelve children and adolescents aged between 8 and 16 years who were receiving CBT for the first time participated in the present study. None of the participants had been treated previously with antidepressants or other psychotropic drugs, or with psychological therapies. Patients were diagnosed using the Diagnostic and Statistical Manual of Mental Disorders-V (DSM-V).²⁵ The study was carried out at the Child and Adolescent Psychiatry and Psychology Service of the Institute of Neuroscience at the Hospital Clinic of Barcelona. The study was naturalistic and the treatment was established according to the Clinical Guidelines for the Treatment of Obsessive-Compulsive Disorder of the Child and Adolescent Psychiatry and Psychology Service. All procedures were approved by the Hospital Clínic ethics committee. Written informed consent was obtained from all the parents and verbal informed consent was given by all the participants following explanation of the procedures involved. All experiments were performed in accordance with relevant guidelines and regulations. This study was conducted in accordance with the Declaration of Helsinki.

Cognitive-behavioral therapy counseling consisted of one session that covered the psycho-educational aspects of OCD (nature of OCD, clinical characteristics and principles of behavior therapy) and twelve sessions (two sessions every week during the first month and a weekly session during the second month) based on exposure and response prevention.

Information on illness severity was obtained during the initial phase of the study using the Children’s Yale-Brown Obsessive Compulsive Scale (CYBOCS).²⁶ The same scale was administered after 8 weeks of CBT. Treatment response was evaluated using the percentage of improvement as follows: ((CYBOCS_8weeks- CYBOCS_basal)/ CYBOCS _basal)*100. Patients were classified as responders or non-responders according to the percentage of improvement after 8 weeks of CBT. Responders were patients with an improvement > 35%, while non-responders were those with an improvement < 25%. Patients with an improvement > 25% and < 35% were classified as partial responders.²⁷

Biological Samples

Two blood samples from each participant were collected before the start of CBT: one in EDTA (BD Vacutainer K2EDTA tubes; Becton Dickinson, Franklin Lakes, New Jersey, USA) and another in PAXgene Blood RNA tubes (Qiagen, Valencia, CA). Genomic DNA was extracted using the MagNA Pure LC DNA Isolation Kit and a MagNA Pure LC 2.0 instrument (Roche Diagnostics GmbH, Mannheim, Germany). DNA concentration and quality were measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Surrey, CA). Genome-wide DNA methylation was profiled at the CEGEN-PRB3-ISCIII using the Illumina Infinium MethylationEPIC BeadChip Kit. Total RNA was isolated in accordance with the manufacturer’s protocol (PAXgene Blood RNA kit). RNA quality and quantity were measured by an Agilent Bioanalyzer 2100 (Santa Clara, CA). 1 μg of purified RNA from each of the samples was submitted to the Kompetenzzentrum fur Fluoreszente Bioanalytik Microarray Technology (KFB, BioPark Regensburg GmbH, Regensburg, Germany) for labeling and hybridization using Human Genome U219 array plates (Affymetrix, Santa Clara, CA, USA), following the manufacturer’s protocols.

Genome-Wide DNA Methylation Profiling

Raw intensity data (.IDAT) files were received and bioinformatics processes were conducted in house using the Chip Analysis Methylation Pipeline (ChAMP) Bioconductor package.²⁸ Raw IDAT files were used to load the data into the R environment with the champ.load function, which also allows for the probe QC and removal steps to occur simultaneously. Probes with weak signals (p < 0.01) (n = 3103), cross-reactive probes (n = 11), non-CpG probes (n = 2952), probes with < 3 beads in at least 5% of the samples per probe (n = 10,683), probes that bound to SNP sites (n = 96,500), and sex chromosome probes (n = 61,734) were all considered problematic for the accurate detection of downstream methylation. After removing these probes, 736,109 probes remained for downstream analysis. β values were then normalized using the champ.norm function, specifically with the beta-mixture quartile method (BMIQ function). Cell counts were measured using the champ.refbase function. The following cells were counted: CD8+ T cells, CD4+ T cells, natural killer (NK) cells, B cells, monocytes, and granulocytes. Next, the singular value decomposition (SVD) method was performed by champ.SVD to assess the amount and significance of the technical batch components, along with any potential confounding variables (eg, sex, age, and cell count) in our dataset. Using the champ.runCombat function, Combat algorithms were applied to correct for slide and array (significant components detected by the SVD method). No effect of sex, age or cell count was detected.

After filtering, normalization, and the detection of batches and covariates, differentially methylated positions (DMPs) were identified using the champ.DMP function, which implements the limma package to calculate the p-value for differential methylation using a linear model (FDR-adjusted p-values < 0.05). An absolute value of the difference between the β-value medians (Δβ) of responders and non-responders higher than 0.2 was set as the cut-off value to decrease the number of significant CpGs and identify sites with the more biologically relevant methylation differences.

Genome-Wide Expression Analysis and the WGCNA Procedure

Microarray data preprocessing was performed using the Babelomics 5 suite (http://www.babelomics.org/).²⁹ The data were standardized using robust multichip analysis. Multiple probes mapping to the same gene were merged using the average as the summary of the hybridization values. Co-expression modules were identified using the R software package for weighted gene co-expression network analysis (WGCNA).³⁰ The co-expression analysis involved constructing a matrix of pairwise correlations between all pairs of genes across all selected samples. Next, the matrix was raised to a soft-thresholding power (β = 8 in this study) to obtain an adjacency matrix. To identify modules of co-expressed genes, we constructed the topological overlap-based dissimilarity, which was then used as input to average linkage hierarchical clustering. This step resulted in a clustering tree (dendrogram) whose branches were identified for cutting based on their shape, using the dynamic tree-cutting algorithm. The above steps were performed using the automatic network construction and module detection function (blockwiseModules in WGCNA), with the following parameters: minModuleSize of 30, reassignThreshold of 0, and mergeCutHeight of 0.25. The modules were then tested for their associations with the trait by correlating module eigengenes with trait measurements.

We then used ClueGO v2.1, a Cytoscape plug-in, to perform a gene set enrichment analysis, as described previously.³¹ Briefly, we selected the unstructured terms of biological processes from Gene Ontology (GO). Only terms with an adjusted p-value < 0.05 and experimental evidence were selected for analysis. Genes involved in significant modules were mapped to their enriched term based on the hypergeometric test (two-sided), with the p-value being corrected by the Benjamini-Hochberg method. ClueGO created a functional module network in which the different GO terms were clustered according to the strength of the association between the terms calculated using chance-corrected kappa statistics.

Statistical Analysis

Data were analyzed using the SPSS 22.0 software (IBM, Chicago, IL, USA). The normality of continuous variables was tested using the Kolmogorov–Smirnov and Shapiro–Wilk tests, while the equality of the variance between the groups was assessed using Levene’s test. Two-tailed p-values < 0.05 were considered to be of statistical significance. In genes enriched with DMPs significantly associated with CBT response (Δβ > 0.2, FDR-adjusted p-values < 0.05), the β values of the most significant DMPs in each gene were tested for correlation, using Spearman’s rank correlation coefficient, with the eigengene values of the modules significantly associated with CBT response in the WGCNA.

Results

An overview of the study design is shown in Figure 1. Table 1 shows the demographic and clinical characteristics of the study participants. As can be observed, there were non-significant differences between the responders and non-responders for age, sex, symptom severity at baseline and family history of OCD. Although non-significant, a higher percentage of non-responders (100%) than responders (50%) presented co-morbidities.

Table 1 Demographic and Clinical Characteristics of the Study Participants

Figure 1 Overview of the study.

Genome-Wide Methylation Analysis

We classified 55,149 probes as significant DMPs (FDR-adjusted p-values < 0.05). However, this included DMPs with very small differences in methylation between responders and non-responders. Therefore, a cut-off of Δβ > ± 0.2 was applied, which identified 89 DMPs with methylation changes that were more likely to be biologically relevant (Supplementary Table 1).

The 89 significant CpGs mapped to 70 genes. Two of these genes were enriched with significant DMPs (FDR-adjusted p-value < 0.05, Δβ > ± 0.2) and were selected for subsequent analysis: PIWIL1 and MIR886. PIWIL1 was enriched with five CpGs that were significantly hypermethylated in the non-responders. These DMPs were upstream of the transcription start site (from +1500 to +200 bp), in a region that, according to the UCSC browser, includes a CpG island and a DNase hypersensitive site (Figure 2A). The most significant CpG in PIWIL1 (cg13861644) is included in the Blood Brain DNA Methylation Comparison Tool (https://epigenetics.essex.ac.uk/bloodbrain/),³² showing a significant correlation (p < 0.05) between methylation levels in the blood and in the prefrontal cortex (r = 0.76), entorhinal cortex (r = 0.83), superior temporal gyrus (r = 0.77) and cerebellum (r = 0.73) (Supplementary Figure 1).

Figure 2 (A) Distribution of significant DMPs (FDR-adjusted p-value < 0.05, Δβ > ± 0.2) in the PIWIL1 gene and methylation β values in responders and non-responders. (B) Distribution of significant DMPs (FDR-adjusted p-value < 0.05, Δβ > ± 0.2) in the MIR886 gene and methylation β values in responders and non-responders. (C) Module eigengene values (y-axis) for the yellowgreen module in individual samples (x-axis). Black bars indicate non-responders, while gray bars indicate responders. (D) Scatter plots showing correlations between yellowgreen module eigengene values (x-axis) and methylation β values of the cg13861644 in PIWIL1 (y-axis). Black points correspond to non-responders, while gray points correspond to responders. (E) Scatter plots showing correlations between yellowgreen module eigengene values (x-axis) and methylation β values of the cg04481923 in MIR886 (y-axis). Black points correspond to non-responders, while gray points correspond to responders.

MIR886 was enriched with four DMPs that were significantly hypomethylated in the responders. These CpGs were upstream of the transcription start site (from +1500 to +200 bp), a region that, according to the UCSC browser, includes a promoter region enriched with H3K27AC marks in all the cell lines considered by ENCODE (Figure 2B). The four CpGs in MIR886 are included in the Blood Brain DNA Methylation Comparison Tool, showing a significant correlation (p < 0.05) between methylation levels in the blood and in the prefrontal cortex (r > 0.89), entorhinal cortex (r = 0.95), superior temporal gyrus (r > 0.92) and cerebellum (r > 0.52) (Supplementary Figure 2).

Genome-Wide Expression Analysis

We applied WGCNA to genome-wide expression data, which identified 70 gene co-expression modules (Supplementary Figure 3). One module, the yellowgreen (197 genes), showed a significant correlation with CBT response (r = −0.85, FDR-corrected p-value = 0.0003). The yellowgreen module contained genes that were upregulated in non-responders to CBT (Figure 2C).

To explore the biological mechanism associated with the genes of the yellowgreen module, we performed a gene set enrichment analysis using the unstructured terms of biological processes from Gene Ontology (GO). We identified five clusters involving ten significant terms (Bonferroni-corrected p-value < 0.05) (Table 2) that were related to DNA replication, chemotaxis, hormone metabolism and catecholamine transport.

Table 2 Gene Set Enrichment Analysis of Biological Processes from Gene Ontology (GO) Obtained for the Yellowgreen Module. The Table Shows the GO Terms Identified, Their Cluster Distribution According to ClueGO, Their Bonferroni-Corrected p-values and the Associated Genes Found in the Yellowgreen Module

Integrative Analysis of Methylation and Gene Expression

We next investigated the possible relationship between the differences in DNA methylation between the responders and non-responders and the gene co-expression modules that were associated with CBT response. We analyzed the correlation between the β values of the most significant DMP in the PIWIL1 and MIR886 genes and the module eigengene values. There were significant correlations between the yellowgreen module and the cg13861644 in PIWIL1 (r = 0.74, p = 0.005) and the cg04481923 in MIR886 (r = 0.72, p = 0.008). Patients showing higher methylation in these CpGs showed an upregulation of the genes in the yellowgreen module (Figure 2D and E).

We also analyzed the correlation between the β values of the most significant DMP in the PIWIL1 and its expression. The Human Genome U219 array plates only includes probes for the PIWIL1 gene but not for the MIR886. Non-significant correlation between methylation and expression was detected between cg13861644, the most significant DMP in the PIWIL1 gene, and its expression in the microarray. Although non-responders showed lower gene expression of PIWIL1 (4.5±1.9), in agreement with the observed hypermethylation, than responders (5.3±0.3), the difference was not significant (p>0.05).

Discussion

To our knowledge, the present study is the first to analyze and integrate differences in DNA methylation and gene expression in association with CBT response in the peripheral blood of children and adolescents with early-onset OCD. We identified two genes, PIWIL1 and MIR886, that were enriched in significant CpG sites that showed meaningful differences (Δβ > ± 0.2) in the methylation level between responders and non-responders and also a strong correlation in DNA methylation between the blood and brain. These CpGs showed higher methylation levels in non-responders and were associated with a module of 197 genes that were co-expressed and upregulated in the non-responders. Interestingly, PIWIL1 and MIR886 are involved in the tight control of gene expression with non-coding RNAs (ncRNAs). Small ncRNAs have roles in neuronal function, cognition, learning and memory.³³

PIWIL1 encodes a Piwi-like protein that forms an evolutionarily-conserved gene regulatory mechanism together with the Piwi-interacting RNAs (piRNAs), a class of small ncRNAs. Piwi proteins and piRNAs are found primarily within the male germline, where they are necessary for germ cell maintenance and spermatogenesis because they protect the genome by silencing transposon expression at both the epigenetic and post-transcriptional levels.^34–37 In addition to their role in germline genome defence, there is growing recognition that the Piwi pathway is involved in neuronal development throughout the lifespan and in neuronal gene regulation in the adult brain.^38–42 Moreover, functional disruption of the Piwi pathway has indicated that it is also involved in learning and memory and in the regulation of behavioral responses to the environment.⁴³ These findings are consistent with the strong association between coding mutations in the Piwi genes in humans and autism.⁴⁴

MIR886 is a Pol III non-coding RNA 886 gene (nc886), which was previously proposed to encode a pre-miR-886 or an RNA component of the vault complex referred to as vtRNA2-1.⁴⁵ However, a later study did not find any evidence that nc886 gives rise to microRNAs or that it associates with the vault complex.⁴⁶ This gene was previously shown to be elevated in Friedreich’s ataxia and differentially methylated in Parkinson’s disease.^47–49 nc886 has a CpG island in its upstream region that is maternally imprinted.⁵⁰ Genomic imprinting is the monoallelic expression of a subset of genes in a conserved, parent-of-origin fashion. The frequency of imprinting of the nc886 CpG island in children has been associated with the genetic background and has also been linked to the mother’s age and season of conception, indicating that genetic and environmental factors may affect the establishment of imprinting, which is closely associated with human physiology.^50,51 Changes in gene expression of imprinted sites within the placenta, including of MIR886, that are suggestive of an altered imprinting status have been linked to newborn neurobehavioral outcomes.⁵²

The genes in the yellowgreen module are associated with several biological processes such as DNA replication, chemotaxis, hormone metabolism and catecholamine transport. These results agree with those of one of the two studies using genome-wide expression analysis of CBT response in anxiety disorders, which identified similar GO terms of “DNA transcription and elongation” and “positive regulation of chemotaxis”.¹²

Although correlations between DNA methylation in promoter regions and gene expression have been reported,⁵³ in our study we did not observed this effect. This could be due to the small sample size of our study. However, it could also be related to the complex mechanisms implicated in the epigenetic regulation of gene expression. The hypermethylation observed in the promoter region of PIWIL1 could not affect the basal expression of this gene but could modify its regulation by transcription factors that participate in the modulatory effects exerted by CBT therapy.

The findings of this study should be interpreted by bearing in mind several important limitations. The sample size limited the statistical power of the study and made it difficult to detect small or modest effects on DNA methylation or gene expression. Given that the study was hypothesis-driven and due to the small sample size, our results should be seen as preliminary and should be considered as exploratory findings that require further confirmation. However, it should be noted that our sample comprised patients with early-onset OCD. Thus, the sample represented a homogeneous clinical population who had not been previously treated and who were at the initial stages of the illness. Moreover, several potential confounders were controlled for, such as age, smoking status, pharmacological treatment and the course of the disease. We also controlled for blood cell composition, as DNA methylation is cell type-specific and different cell compositions between samples could affect the methylation data obtained.

In conclusion and despite the study limitations, our results provide evidence that the epigenetic regulation of ncRNAs could be a predictor of CBT response and might be related to differences in the expression of genes involved in biological processes associated with CBT response. Our results have to be replicated in large samples before using the methylation level of these specific genes as predictive biomarkers with clinical application.