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.

Methods

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 single-blind 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 biomarker-positive 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-D17) (Hamilton, 1960). Degree of improvement was defined as the change in the Ham-D17 score from pretreatment baseline to endpoint, with response defined as a decrease in Ham-D17 of ≥50% and remission as a final Ham-D17 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 for a breakdown of samples by treatment and time point).

Table 1

Number of samples by time point and medication.

	Escitalopram (ESC)	Bupropion (BUP)	Total
Baseline	40	24	64
Week 1	34	23	57
Week 2	35	26	61
Week 8	36	28	64
Total	145	101	246

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 Agilent Bioanalyzer 2100 (Agilent, Santa Clara CA) in the UCSF Genome Core Facility. Total mRNA was amplified and labeled using the Illumina TotalPrep RNA Amplification Kit before being hybridized to the Illumina Expression RefSeq HT-12 BeadChip, querying the expression of ∼47,000 RefSeq transcripts, as per manufacturer’s protocol. The BeadChips were then scanned with the Illumina iScan system and signal was extracted using the Illumina BeadStudio Software (Illumina, San Diego, CA).

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 log2 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 log2 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 over-represented 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).

Reference Gene Sets From Relevant Published Studies

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.

Results

Clinical Outcome

Of the 68 subjects, 67.7% responded and 47% remitted with treatment with a mean decrease in HamD17 score of 60% (); 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 .

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 (). 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.

Figure 1

Venn diagrams of genes differentially expressed over time, and the intersection of genes between each comparison. Differential expression analyses were conducted to compare expression levels within patients and between the indicated timepoints: BL = baseline, TP1 = week 1, TP2 = week 2, TP8 = week 8 after pharmacotherapy; e.g., TP1vs BL indicates those genes differentially expressed between TP1 as compared to BL. The number of genes differentially expressed at a threshold of p < 0.005 between each pair of time points are indicated. Those genes showing a change in expression across multiple comparisons are indicated in the intersections of each diagram. In (A), the number of genes showing increased expression is indicated in red, while the number of genes showing decreased expression is indicated in green. The total number of differentially expressed genes is shown in (B).

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 . 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 2A

Enrichment of top 10 gene ontology biological process terms in the global treatment effect gene set.

Gene Ontology biological process terms (N genes overlapping)	FDR corrected q-value
Translational termination (14)	0.46
Viral life cycle (21)	0.64
Cellular metabolic process (244)	0.64
Cellular component disassembly (31)	0.64
Cellular protein complex disassembly (16)	0.64
Nuclear-transcribed mRNA catabolic process (13)	0.64
mRNA metabolic process (27)	0.64
Translational elongation (13)	0.64
mRNA catabolic process (13)	0.64
Submandibular salivary gland formation (2)	0.64

Table 2B

Enrichment of top 10 gene pathways in the global treatment effect gene set.

Gene pathways (N genes overlapping)	FDR corrected q-value
Reactome pathway: influenza viral RNA transcription and replication (10)	0.08
Reactome pathway: eukaryotic translation and elongation (9)	0.08
Reactome pathway: influenza life cycle (11)	0.08
Reactome pathway: influenza infection (11)	0.08
Reactome pathway: translation (11)	0.08
Reactome pathway: GTP hydrolysis and joining of the 60S ribosomal subunit (9)	0.08
Reactome pathway: peptide chain elongation (8)	0.08
Reactome pathway: eukaryotic translation and termination (8)	0.08
Reactome pathway: viral mRNA translation (8)	0.08
Reactome pathway: SRP-dependent co-translational protein targeting to membrane (9)	0.08

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 cis–trans 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 . 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 (). 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 3

Eighteen genes overlapping between treatment-related gene list and Guilloux et al., 2015.

Gene symbol
ABCC3	ATP-binding cassette, sub-family C (CFTR/MRP), member 3
ADAR	Adenosine deaminase, RNA-specific
ADIPOR1	Adiponectin receptor 1
EPOR	Erythropoietin receptor
FAM136A	Family with sequence similarity 136, member A
GSR	Glutathione reductase
LSM7	LSM7 homolog, U6 small nuclear RNA and MRNA degradation associated
MARCKS	Myristoylated alanine-rich protein kinase C substrate
NDUFB2	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 2, 8kDa
PBX1	Pre-B-cell leukemia homeobox 1
R3HDM4	R3H domain containing 4
RASSF7	Ras association (RalGDS/AF-6) domain family (N-terminal) member 7
ROMO1	Reactive oxygen species modulator 1
RPL27	Ribosomal protein L27
RPS10	Ribosomal protein S10
SP2	Sp2 transcription factor
SUGP2	SURP and G-patch domain-containing 2
TMEM208	Transmembrane protein 208

Table 4

Thirteen genes overlapping between treatment-related gene list and genes selected by Wisniewska et al. (2012).

Gene symbol
DPYSL5	Dihydropyrimidinase-like 5
FGF10	Fibroblast growth factor 10
GABRG2	Gamma-aminobutyric acid (GABA) A receptor, gamma 2
HOXC6	Homeobox C6
ITGA7	Integrin, alpha 7
KDM2B	Lysine (K)-specific demethylase 2B
MARCKS	Myristoylated alanine-rich protein kinase C substrate
QRFP	Pyroglutamylated RFamide peptide
RAP2B	RAP2B, member Of RAS oncogene family
SDPR	Serum deprivation response
SMYD2	SET and MYND domain-containing 2
SYT5	Synaptotagmin V
TPM1	Tropomyosin 1 (alpha)

Table 5

Genes showing differential expression in current study and at least two previous studies.

Gene name	Number of previous studies showing differential expression	Study 1	Study 2	Study 3	Study 4	Gene function
map2k3	2	Belzeaux2012suppletable.txt	Pettai et al., 2016 mmc1.cvs			Signaling pathway (dual specificity kinase)
sp2	3	Belzeaux2012suppletable.txt	Belzeaux et al., 2016 mmc3.csv	Guilloux et al., 2015		Transcriptional factor (activates mRNA synthesis)
atxn10	2	Belzeaux2012suppletable.txt	Mamdani et al., 2011 tp201112x1.csvtp201112x2.csv			Signaling pathway (activates the Ras-MAP kinase pathway)
chchd1	2	Belzeaux et al., 2016 mmc3.csv	Mamdani et al., 2011 tp201112x1.csvtp201112x3.csv			Mitochondrial translational factor
serac1	2	Belzeaux et al., 2016 mmc3.csv	Mamdani et al., 2011 tp201112x1.csvtp201112x3.csv			Mitochondrial phospholipid remodeling
c6orf106	2	Belzeaux et al., 2016 mmc3.csv	Mamdani et al., 2011 tp201112x1.csvtp201112x2.csv			Brain function unknown(expression associated with poor human cancer prognosis)
kdm2b	2	Belzeaux et al., 2016 mmc3.csv	Wisniewska et al., 2012 1471-2164-13-635-S1.xls			Wnt/β-catenin pathway in thalamusRegulation of histone methylation
smyd2	2	Belzeaux et al., 2016 mmc3.csv	Wisniewska et al., 2012 1471-2164-13-635-S1.xls			Wnt/β-catenin pathway in thalamusRegulation of histone methylation
MARCH8	2	Hodgson et al., 2016 Table S1	Mamdani el al., 2011 2011-12-16203-Supplementary_table4.xls			Ubiquitin protein ligase
MARCKS	4	Mamdani et al., 2011 tp201112x2.csv	Wisniewska et al., 2012 1471-2164-13-635-S1.xls	Redei et al., 2014.pdf	Guilloux et al., 2015	Wnt/β-catenin pathway in thalamusSignaling (most prominent substrate for protein kinase C)

There was no significant overlap between our gene list and previous preclinical studies (above) of drug-specific genomic signatures in the mouse brain examining the effects of either BUP or fluoxetine.

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 (). The differences in 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 data-driven 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 (). 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 (). 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 (). 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.

Author Contributions

Performed treatment study: IAC, AFL. Performed genetic analyses: SPH. Performed data analyses: IAC, EC, GC, AMH. Wrote manuscript: AFL, DEK, IAC.

Funding

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

Conflict of Interest Statement

IC discloses that he has received research support from Covidien (formerly Aspect Medical Systems), National Institutes of Health, and NeoSync, Inc.; he has been an advisor/consultant/reviewer for Arctica Health, Cerêve, HeartCloud, NeuroDetect, NeuroSigma, NIH (ITVA), U.S. Departments of Defense and Justice, and the VA (DSMB); he is editor of the Patient Management section of the American Psychiatric Association’s FOCUS journal; his biomedical intellectual property is assigned to the Regents of the University of California, and he has stock options in NeuroSigma, where he has served as Chief Medical Officer. GC discloses that he has received research funding from Takeda Pharmaceutical Company Ltd. for an unrelated project. DK discloses that he has received research support from the National Institutes of Health. He currently serves as an investigator on a study sponsored by Neosync, Inc. AL discloses that he has received research support from the National Institutes of Health, Neuronetics, Department of Defense, CHDI Foundation, and NeuroSigma, Inc. He has served as a consultant to NeoSync, Inc., Ionis Pharmaceuticals, Inc., and ElMindA. He is Chief Scientific Officer of Brain Biomarker Analytics LLC (BBA). AL owns stock options in NeoSync, Inc. and has equity interest in BBA. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.