Time Course of Changes in Peripheral Blood Gene Expression During Medication Treatment for Major Depressive Disorder

Changes in gene expression (GE) during antidepressant treatment may increase understanding of the action of antidepressant medications and serve as biomarkers of efficacy. GE changes in peripheral blood are desirable because they can be assessed easily on multiple occasions during treatment. We report here on GE changes in 68 individuals who were treated for 8 weeks with either escitalopram alone, or escitalopram followed by bupropion. GE changes were assessed after 1, 2, and 8 weeks of treatment, with significant changes observed in 156, 121, and 585 peripheral blood gene transcripts, respectively. Thirty-one transcript changes were shared between the 1- and 8-week time points (seven upregulated, 24 downregulated). Differences were detected between the escitalopram- and bupropion-treated subjects, although there was no significant association between GE changes and clinical outcome. A subset of 18 genes overlapped with those previously identified as differentially expressed in subjects with MDD compared with healthy control subjects. There was statistically significant overlap between genes differentially expressed in the current and previous studies, with 10 genes overlapping in at least two previous studies. There was no enrichment for genes overexpressed in nervous system cell types, but there was a trend toward enrichment for genes in the WNT/β-catenin pathway in the anterior thalamus; three genes in this pathway showed differential expression in the present and in three previous studies. Our dataset and other similar studies will provide an important source of information about potential biomarkers of recovery and for potential dysregulation of GE in MDD.


INTRODUCTION
The mechanism of action (MOA) of antidepressant medications remains incompletely understood. Antidepressant medications are remarkably pleiotropic in their effects. In addition to binding to serotonin and/or norepinephrine transporters (SERT and NET, respectively), their affinity for other neurotransmitter receptors and ligand-gated ion channels may be relevant to their effects (Bianchi and Botzolakis, 2010). In addition, the acute effect on aminergic signaling has the potential to influence a wide range of downstream genetic pathways implicated in neurogenesis, synaptic plasticity, neuronal excitability, and metabolism (Baudry et al., 2011). To date, the identity of genes that are consistently up-or down-regulated by antidepressants remains obscure. The application of genomic techniques to patients undergoing treatment with antidepressants has the potential to identify clusters of genes and/or pathways that undergo transcriptional regulation in response to these drugs. Identification of reproducible gene expression (GE) changes during antidepressant treatment will be the first step toward the long-term goal of elucidating which changes are required for the relief of depressive symptoms (Bartova et al., 2010;Baudry et al., 2011;Redei et al., 2014).
Because of the limited accessibility of neuronal tissue, studies of the transcriptome in patients undergoing treatment for MDD have focused upon GE changes in peripheral blood (Tylee et al., 2013). Studies often have been of limited size but have examined changes during treatment with a variety of medications (Belzeaux et al., 2010;Mullins et al., 2014;Watanabe et al., 2015;Miyata et al., 2016) or cognitive behavioral therapy (Redei et al., 2014). To our knowledge, no previous study has systematically examined the changes in peripheral blood GE at multiple time points in an 8-week course of antidepressant treatment. We report here a study of the transcriptional changes that occur in response to treatment with escitalopram or bupropion in a cohort of 68 subjects at 1, 2, and 8 weeks of medication treatment, and the association of these changes with clinical outcome.

Subjects
Samples were collected from a subset of the adult subjects (age 21-75) who participated in a double-blind controlled trial conducted to evaluate a neurophysiologic biomarker for use in assigning antidepressant medication treatment (NCT00917059 "Personalized Response Indicators of SSRI Effectiveness in Major Depression" [PRISE-MD]). The objective of the overall trial was to assess outcomes when subjects were prospectively assigned to antidepressant treatment based upon biomarker status, using a quantitative electroencephalographic (qEEG) biomarker that previously had been shown to be predictive of response and remission during treatment with escitalopram (ESC) and bupropion (BUP) (Leuchter et al., 2009a;Leuchter et al., 2009b;Cook et al., 2013). All subjects started with 1 week of singleblind ESC (10 mg each morning) for measuring the biomarker and then were randomized to continue on ESC or switch to BUP (150 mg each morning) under double-blind conditions for the remainder of the protocol. Randomization was performed using a biomarker-based stratification to ensure that biomarkerpositive and biomarker-negative subjects received balanced allocation to the treatments. The study tested the hypothesis that treatment concordant with the biomarker would lead to better outcomes than biomarker-discordant assignment. ESC outcomes were assessed after 7 weeks of continuous treatment (study week 7) while BUP outcomes were assessed after 7 weeks of continuous BUP (study week 8). The study was statistically powered to detect a difference for subjects assigned to ESC but not for those receiving BUP. In accordance with principles of the Helsinki Declaration, all protocols were approved by the UCLA Institutional Review Board, and all subjects provided written informed consent.
All subjects were diagnosed with MDD based upon structured interview data (Mini-International Neuropsychiatric Interview, or MINI) (Sheehan et al., 1998) and had a symptom severity score of 12 or greater on the Quick Inventory of Depressive Symptomatology-Self Rated version (QIDS-SR16) (Rush et al., 2003). Subjects were excluded for serious medical or psychiatric comorbidities, recent exposure to electroconvulsive therapy, and other factors that would either pose safety concerns or render data uninterpretable (detailed at clinicaltrials.gov). Subjects were evaluated in an outpatient setting every week for evidence of symptomatic changes and/or side effects. The primary outcome measure was the score on the 17-item Hamilton Depression Rating Scale (Ham-D 17 ) (Hamilton, 1960). Degree of improvement was defined as the change in the Ham-D 17 score from pretreatment baseline to endpoint, with response defined as a decrease in Ham-D 17 of ≥50% and remission as a final Ham-D 17 of ≤7.
A total of 274 adults were screened, 180 were enrolled and randomized, and 133 completed at least 4 weeks of the trial and had evaluable clinical data, with usable samples for genomic assays collected from a subset of 68 individuals (40 ESC, 28 BUP) at up to four time points during the trial: pre-treatment baseline and 1, 2, and 8 weeks of treatment in the protocol (see Table 1 for a breakdown of samples by treatment and time point).

RNA Extraction and Measurement of Gene Expression
Whole blood was collected using PAXgene Blood RNA tubes (PreAnalytiX GmbH, Hombrechtikon, Switzerland) and frozen until time of extraction. Samples were extracted in one batch using a PAXgene Blood RNA Kit (Qiagen, Valencia CA); quantity and quality of RNA samples were checked using an Illumina expression assays were run on blood samples meeting standard QC criteria. After assay, raw data were inspected to remove samples with missing data for the majority of probes. Data were inspected after log 2 transformation to screen for outliers, using visual inspection of correlations among samples and probes and calculating z-scores to find samples with low correlations with the rest of the dataset. A total of 14 samples were excluded based on QC of raw data, and the resulting 246 samples were carried forward for analysis.
Raw GE data were preprocessed using Bioconductor packages for R following protocols developed by the UCLA Informatics Center for Neurogenetics and Neurogenomics (ICNN). Expression data were log 2 transformed and quantile normalized; no batch correction was necessary as all samples were processed at the same time. As multiple probes may tag the same gene on the array, the collapseRows function was used to select a single probe per gene, using the absolute maximum mean expression values of the multiple probes. Quality-control analysis was performed by examining the inter-array Pearson correlations, clustering based on variance, and the mean absolute deviation (MAD) using the top 1,000 most variant probes. Twelve outliers were identified based on GE data and removed from further analyses. After quality control steps were applied, 28,105 gene transcripts were present in 246 samples from 68 unique subjects.

Data Analysis
Analyses of demographic and clinical outcome data were performed using SPSS (SPSS, Inc.; Chicago, IL). Outcome and demographic differences were compared between outcome groups using Student's T-test. Analysis of differential GE was performed using a linear model fitting (LIMMA package), accounting for within-subject variance. Age, sex, and RNA integrity number (RIN) were included as covariates in all analyses. After linear model fitting, a Bayesian estimate of differential expression was calculated applying a nominal p-value threshold of 0.005. The resulting gene lists were interrogated using the data-driven tools described below.

Gene Ontology and Gene Pathways
We used the anRichment function within the WGCNA package in R (Langfelder and Horvath, 2008;Miller et al., 2011) to test for enrichment of: (1) gene ontology (GO) biological process terms and (2) pathway-based categories (e.g., NCBI BioSystems Reactome pathways). Enrichment analyses are conducted using a hypergeometric test, which can be used to identify whether a sub-population (e.g., differentially expressed genes) is overrepresented in a given sample (e.g., genes annotated with a specific function), when drawing from the larger population of probes. This test is based on a discrete probability distribution, representing the likelihood of obtaining a number of successful draws, out of N possible draws, without replacement from a finite population. To adjust for multiple comparisons, we applied a local false discovery rate (FDR) correction and report q-value estimates using the R q-value package for all enrichment tests.

Cell Type-Specific Enrichment Analyses
We queried data available in the RNA-seq database (Zhang et al., 2014) to identify gene sets preferentially enriched in cell types highly represented in the nervous system, including astrocytes, endothelial cells, microglia, neurons, and oligodendrocytes (further subdivided by precursors, newly formed, and myelinated oligodendrocytes). Following the authors' recommendations, the top 500 enriched genes (using FPKM > 20) were selected for each cell type. Because of evidence that the anterior thalamus may constitute a site for antidepressant medication action, we tested for enrichment of genes preferentially expressed in the thalamic structures (Hawrylycz et al., 2012), as well as genes involved in the Wnt/β-catenin pathway in the anterior thalamus (Wisniewska et al., 2012).
After a comprehensive review of the literature, published studies relevant to the current data set were compiled. This includes studies focusing on GE differences in MDD patients vs. controls, or within MDD patients as a function of treatment (Mamdani et al., 2011;Belzeaux et al., 2012;Menke et al., 2012;Liu et al., 2014;Mamdani et al., 2014;Guilloux et al., 2015;Hennings et al., 2015;Hodgson et al., 2016;Jansen et al., 2016). This includes studies of peripheral blood GE changes in humans undergoing treatment with citalopram and psychotherapy (n = 34) (Guilloux et al., 2015), citalopram alone (n = 77) (Mamdani et al., 2014), a variety of different psychotropic agents as well as ECT (n = 16) (Belzeaux et al., 2012), nortriptyline or ESC (n = 136) (Hodgson et al., 2016), or citalopram (n = 34) (Belzeaux et al., 2016). We examined genes sets differentially expressed in these previous studies to determine whether there was significant overlap with the present study, as well as to determine whether individual genes were reported in more than one previous study of antidepressant treatment, as well as the Wnt/β-catenin pathway.

Drug-Induced Gene Changes in Mice
Using the Genomic Signature Identification tool in the genes2mind database, we extracted lists of genes that represent drug-specific genomic signatures in the mouse brain following antidepressant administration. Specifically, we extracted the top 100 genes associated with BUP or fluoxetine at 1, 2, 4, and 8 hours and all time points.

Exploratory Pathway Analysis
Ingenuity Pathway Analysis (IPA) was used to explore biological networks enriched by genes identified through the differential expression analysis. Follow-up exploratory analyses were conducted to examine differences as a function of treatment type and response status.

Clinical Outcome
Of the 68 subjects, 67.7% responded and 47% remitted with treatment with a mean decrease in HamD 17 score of 60% (Table 1); by treatment group, 60.0% (37.5%) of ESC subjects and 78.6% (60.7%) of BUP subjects responded (remitted), respectively. There was a trend toward a greater proportion of male patients and a greater degree of improvement in the group of patients treated with BUP, although these differences were not statistically significant between the two treatment conditions. Thirty-six of the ESC-and 28 of the BUP-treated subjects completed 8 weeks of treatment. The number of subjects who completed treatment who had blood drawn for GE studies at each time point is shown in Table 1.

Gene Expression Analyses
The comparison of baseline to week 1, 2, and 8 samples revealed 156, 121, and 585 differentially expressed probes, respectively, at an uncorrected threshold of p < 0.005 (Figure 1). We did not explore differences between groups treated with ESC versus BUP given the small sample sizes when stratifying by treatment. GE data therefore were collapsed across treatments to examine the effects of treatment on GE irrespective of the specific drug used.
Only a subset of gene transcripts showed a change in expression at more than one time point, with 7 shared between the first and second weeks, 31 between the first and eighth week, 14 between the second and eighth week, and only 4 across first-, second-, and eight-week time points. To limit the number of comparisons in further data analyses, we restricted our detailed investigations to the most robust effect of those transcripts showing differential expression: 8 weeks compared to baseline (585 genes).
An examination of GO.BP categories revealed increased expression for several process terms, although none of these terms were statistically significant at the 0.05 level, after correction for multiple comparisons with FDR. The top 10 GO.BP process terms are shown in Table 2A. To supplement enrichment analyses of GO terms, we also tested for enrichment of gene pathways using a compiled set of lists from CHDI. While there was increased expression of several of these pathways, none were significantly enriched after correction for multiple comparisons (FDR p > 0.05). The top 10 pathway terms are presented in Table 2B.
To further interrogate these results, we submitted the list of 585 differentially expressed genes (at uncorrected p < 0.05) for annotation enrichment to DAVID (at an uncorrected threshold of p < 0.05). Five biological process terms were significantly enriched: regulation of cell proliferation, GE, cellular biosynthetic process, positive regulation of cell proliferation, and positive regulation of biological process. Two molecular function terms were significantly enriched, specifically RNA binding and cistrans isomerase activity.
Examination of gene sets preferentially enriched for specific cell types within the nervous system did not reveal a significant overlap between our global treatment effect gene list and genes preferentially expressed in any of these cell types.
Examination of the results of nine previous peripheral GE studies of patients undergoing antidepressant treatment and with published gene lists available for comparison revealed a significant overlap of 18 genes (FDR q = 0.046) between our global treatment effect gene list and one other published gene list (Guilloux et al., 2015) comparing non-remitters to remitters at pretreatment baseline. The list of overlapping genes is provided in Table 3. There also was not a significant overlap with genes expressed in the any region of the thalamus, but there was a trend (p < 0.07) toward overlap with genes involved in the Wnt/β-catenin pathway in the anterior thalamus (Table 4). There were a total of 10 individual genes that were identified in the present study and at least two previous studies. This list is shown in Table 5.

DISCUSSION
The results reported here are consistent with those of previous studies and indicate a significant effect of antidepressant medication treatment on GE in peripheral blood. We identified significant transcriptional changes in genes that were up-or down-regulated over the course of 8 weeks of antidepressant treatment with either ESC or BUP, two medications in widespread clinical use (American Psychiatric Association, 2010). The number of genes showing differential expression increased over the course of treatment, with 585 genes showing transcriptional changes by week 8. This emphasizes the importance of considering GE not just at end of treatment but at a series of time points during antidepressant treatment, ranging from the first week through the end of a full course of 8 weeks. There was limited overlap in the genes that were up-or down-regulated at each time point, suggesting that there may be changes in the specific pattern of gene regulation at each time point. In contrast to some previous studies, there was no relationship between changes in GE and clinical outcomes of antidepressant treatment, which may reflect heterogeneity among subjects, so that a change in expression in a particular subset of genes is sufficient to lead to clinical improvement in some individuals but not others.
The purpose of this study was to generate hypotheses and to identify genes to examine using bioinformatics tools. By comparing GE at weeks 1, 2, and 8 to baseline, we obtained 156, 121, and 585 significant differences in expression using an uncorrected threshold of p < 0.005. We selected this threshold to be in line with previous recommendations for preventing bias in this type of analysis (Carvalho et al., 2016). However, there is no convention for selecting an appropriate statistical threshold in genome-wide GE analyses of such samples, because the preprocessing and analysis steps tend to vary widely across groups. While it is common to use a threshold of p < 0.05 with a fold-change greater than 1.5, for example, this is arbitrary. Moreover, the use of 0.05 as the threshold in our study would likely have identified a very large number of genes and included potentially spurious findings. Our use of a threshold of p < 0.005 was motivated partly by our previous experience with similar datasets, and also reasoning that this threshold, while more conservative than p < 0.05, would not restrict the lists of significant genes to a number too small to allow further interrogation with bioinformatics tools.
This dataset is among of the largest of a growing number of studies that have examined GE in peripheral blood associated with antidepressant treatment (Mamdani et al., 2011;Belzeaux et al., 2012;Mamdani et al., 2014;Redei et al., 2014;Guilloux et al., 2015;Hodgson et al., 2016;Belzeaux et al., 2016;Jansen et al., 2016;Pettai et al., 2016). In contrast to several of these studies, we did not detect a significant association between changes in GE and clinical outcome (Mamdani et al., 2011;Mamdani et al., 2014;Redei et al., 2014;Belzeaux et al., 2016;Hodgson et al., 2016;Pettai et al., 2016) or between exposure to different antidepressant medications (Hodgson et al., 2016). It is possible that a larger sample size would be needed to detect these types of differences in GE. It also is possible that the therapeutic effects of antidepressants are not reflected by global GE, or that such changes are not reliably or consistently reflected in peripheral blood cells (as opposed to nervous system tissue). The use of peripheral samples is currently standard in the field and ethical and technical factors prevent ready access to nervous system cell types, except for cadaveric tissue or nasal epithelium (as recently reported for subjects with schizophrenia and autism) (Lepagnol-Bestel et al., 2008;Dong et al., 2012;Abdolmaleky et al., 2014;Choi et al., 2014;Gandal et al., 2018). This approach is less likely to be useful for studies of MDD or the response to antidepressants; while MDD is a chronic and recurrent condition, the phenotype of illness is not always expressed consistently, and response to antidepressant frequently is time-limited.
A subset of the 585 genes showing differential expression at 8 weeks showed similar changes at week 1 (31) and week 2 (14) comparisons with baseline (Figure 1). The differences in Ribosomal protein S10 SP2 Sp2 transcription factor SUGP2 SURP and G-patch domain-containing 2 TMEM208 Transmembrane protein 208   GE over time suggests that there may be different phases of gene regulation over the course of antidepressant treatment, with some genes showing only early transcription changes and fewer showing consistent change over time. GE studies of peripheral blood would be particularly well-suited to detect such transient state-dependent changes. It also is possible, however, that a number of the changes detected at only one or two time points during treatment represent false positive changes in GE. Future studies should examine changes in GE over time to verify the present finding. To characterize the nature of genes identified in our resulting global-treatment-effect gene list, we used a number of datadriven tools, which allowed us to test whether there was an overlap between our gene list and reference gene list(s) greater than expected by chance. Our analyses using GO terms or pathways did not show a significant overlap; by contrast, annotation enrichment to DAVID revealed an enrichment of terms related to several basic biological processes (e.g., GE, cellular biosynthetic process) as well as cell proliferation. This finding is consistent with the observation that antidepressant drugs can affect a variety of different cellular processes (Blier et al., 1998;Miller et al., 2009;Sharp, 2012;Duhr et al., 2014;Stokes et al., 2015;Sun et al., 2015;Björkholm and Monteggia, 2016).
We did not detect an enrichment for any genes that were preferentially expressed in nervous system cell types. The lack of a detectable signal for genes enriched in neurons may at first pass be surprising, because the putative targets of ESC and BUP are neuronally expressed monoamine transporters. The clinical response to these drugs, however, is unlikely to be confined to the particular receptor to which the drugs initially bind. Rather, the long-term, therapeutic response is likely to involve neurons in pathways that are downstream from the initial drug target (Duman et al., 2016). It is also possible that the lack of a signal associated with neuronal or glial cell types may reflect the limitations of using a non-neuronal cell type for our genomic analysis. Peripheral blood cells provide a convenient source for human mRNA and are a common tissue source from living patients. It is likely that some transcriptional pathways are shared between white blood cells and neuronal tissue, but it is even more likely that pathways would be shared with microglia, as both function as elements of the immune system.
We did detect a significant overlap between the genes differentially expressed in the present study and one previous study (Guilloux et al., 2015). These genes involve a variety of cellular processes, most frequently involving gene regulation and transcription (Table 3). This finding is consistent with a putative role for antidepressant medications in increasing neuroplasticity through changes in GE.
We also detected a trend toward an overlap with genes that are differentially expressed in the Wnt/β-catenin pathway in the anterior thalamus (Table 4). Wnt signaling is increasingly recognized as playing an important role in the mature central nervous system, and abnormalities in this signaling pathway have been implicated in MDD (Voleti and Duman, 2012;Zhou et al., 2016;Tayyab et al., 2018). Nuclear β-catenin accumulates in mature neurons throughout the nervous system, particularly anterior thalamic cells where it plays a role in regulating expression of voltage-and ligand-gated ion channels (Wisniewska et al., 2012). Although the precise anatomic site(s) of antidepressant activity are not known, the anterior thalamus has the highest concentrations of serotonin transporter receptor sites in the brain and has been postulated to be a primary site of action of SSRI antidepressant medications such as ESC. Multiple reports have shown changes in cerebral oscillatory activity in prefrontal region during antidepressant treatment on a time frame similar to the GE changes reported here (Cook et al., 2002;Cook et al., 2005;Bruder et al., 2008;Bares et al., 2010;Leuchter et al., 2015;Leuchter et al., 2017).
There were 10 genes showing differential expression in the current study that showed differential expression in at least two previous treatment studies or in the Wnt/βcatenin pathway ( Table 5). Two of these genes are involved in cell signaling pathways, four in gene regulation, three in the Wnt/β-catenin pathway, and one of unknown function. The fact that Wnt/β-catenin pathway genes show expression changes in multiple studies is intriguing because changes in prefrontal rhythmic oscillatory activity have been reported as a reproducible specific biomarker of ESC efficacy. It is tempting to speculate that genes regulating excitability such as voltage-and ligand-gated ion channels in some brain regions might be differentially expressed in patients treated with antidepressants. This possibility remains speculative but may be relevant to GWAS studies linking loci near several channel genes to several psychiatric disorders.
We did not examine GE changes in patients treated with ESC versus BUP because our sample size was likely too small to reliably detect differences in expression associated with particular antidepressants, such as have been reported in larger studies of ESC and nortriptyline (Hodgson et al., 2016). Our findings of a robust change in GE when collapsing across treatments suggest it is possible that the major initial or downstream effects of current antidepressants at a molecular level are relatively similar, despite variability in the response of particular patients to each drug (Rush et al., 2011). Despite differences in the immediate effect of antidepressants on extracellular amine concentrations, it is possible that the downstream therapeutic effects of antidepressants occur via common pathways. If so, the molecular signatures of divergent antidepressants could be quite similar regardless of which molecules they bind to initiate these events.
These results must be interpreted within the context of several limitations of this study. First, while this is among the larger studies of GE during antidepressant treatment, the sample size of 64 subjects limits the statistical power to detect differential expression as well as association with treatment efficacy. Second, the unusual design of the study in which treatment was changed from ESC to BUP in some subjects restricts our ability to attribute changes in expression to any particular antidepressant treatment. Nevertheless, the consistency in GE changes over time and the overlap with one relevant previous study suggests that there may be a shared signal underlying treatment response that is detectable in peripheral blood across studies. Further studies of drug-induced GE will help us to triangulate on loci that are relevant to depression. The power afforded by integrating results across additional datasets may help overcome the limited sample sizes in most extant studies. It is possible that reliable detection of significant changes will require peripheral blood samples from thousands of patients during treatment. If so, very large-scale efforts or alternative methods may be needed to address these questions.

ETHICS STATEMENT
This study was carried out in accordance with the recommendations of the UCLA Institutional Review Board (Medical IRB 3), with written informed consent obtained from all subjects prior to any experimental procedures. All subjects gave written informed consent in accordance with the Declaration of Helsinki; our protocol was reviewed and approved by the UCLA Institutional Review Board and protected this vulnerable population of individuals with depression by means including providing them with the opportunity to review their decision about participation with family, friends, and advisors in advance.

FUNDING
This work was supported by ARRA grant R01MH085925-02S1 (IAC: PI) and intramural resources at the University of California (AFL, SPH).