RNA expression profiling in depressed patients suggests retinoid-related orphan receptor alpha as a biomarker for antidepressant response.

Response to antidepressant treatment is highly variable with some patients responding within a few weeks, whereas others have to wait for months until the onset of clinical effects. Gene expression profiling may be a tool to identify markers of antidepressant treatment response and new potential drug targets. In a first step, we selected 12 male, age- and severity-matched pairs of remitters and nonresponders, and analyzed expression profiles in peripheral blood at admission and after 2 and 5 weeks of treatment using Illumina expression arrays. We identified 127 transcripts significantly associated with treatment response with a minimal P-value of 9.41 × 10(-)(4) (false discovery rate-corrected). Analysis of selected transcripts in an independent replication sample of 142 depressed inpatients confirmed that lower expression of retinoid-related orphan receptor alpha (RORa, P=6.23 × 10(-4)), germinal center expressed transcript 2 (GCET2, P=2.08 × 10(-2)) and chitinase 3-like protein 2 (CHI3L2, P=4.45 × 10(-2)) on admission were associated with beneficial treatment response. In addition, leukocyte-specific protein 1 (LSP1) significantly decreased after 5 weeks of treatment in responders (P=2.91 × 10(-2)). Additional genetic, in vivo stress responsitivity data and murine gene expression findings corroborate our finding of RORa as a transcriptional marker of antidepressant response. In summary, using a genome-wide transcriptomics approach and subsequent validation studies, we identified several transcripts including the circadian gene transcript RORa that may serve as biomarkers indicating antidepressant treatment response.

Introduction

Response to antidepressant treatment is highly variable with 30–40% patients showing insufficient improvement even after administration of a sequence of different antidepressants.[1, 2] Currently available antidepressant drugs act almost exclusively as enhancers of monoaminergic neurotransmission. Thus, the identification of biomarkers that help preselecting patients that will respond to a specific treatment is from the clinical perspective as important as the discovery of new potential drug targets that ultimately may lead to more specific treatment options.

Genotyping of single-nucleotide polymorphisms (SNPs) did not deliver robust antidepressant response predictors that entered clinical routine.[3] The use of high-throughput technologies such as gene expression profiling to identify biomarkers for antidepressant response is a relatively new approach. Apart from a first study in six elderly depressed patients showing various nominal associations,[4] Mamdami et al.[5] found Interferon regulatory factor 7 to be upregulated in citalopram responders of an outpatient depression sample (N=62). Another microarray study of the same group[6] with 77 patients identified SMAD7 and SIGLECP3 as potential predictors of citalopram response. Using a selection of 14 putative candidates emerging from a case–control mRNA and miRNA scan, Belzeaux et al.[7] found a set of four possible predictors of antidepressant treatment response in 16 depressed patients (PPT1, TNF, IL1b and HIST1H1e). Although pointing to some known and some new possible candidates as biomarkers for antidepressant response, all previous studies lack sufficient replications. In addition, particularly mRNA analysis and transcriptomics are subject to technical variability that may explain some of the inconsistencies of previous findings.[8, 9]

Here, we present the results of a genome-wide pharmacotranscriptomics study using highly selected male depressed inpatients matched by age and depression severity at study inclusion. They were selected according to clinical extremes with respect to their response to antidepressants within a 5-week treatment period. Transcripts identified as significantly regulated in this analysis were subsequently analyzed in an independent replication sample of 142 depressed inpatients. In addition, we analyzed a promising candidate gene from the discovery round, retinoid-related orphan receptor alpha (RORa), and performed in vivo stress responsitivity assays and a mouse model of stress-evoked gene expression.

Materials and methods

Sample description and study design

Twenty-four male depressed inpatients participating in the Munich Antidepressant Response Signature project[2] were selected as 12 age- and depression severity-matched pairs. Briefly, in the Munich Antidepressant Response Signature study, at least moderately depressed inpatients were followed during their hospital stay and received a doctor's choice antidepressant medication mimicking the everyday clinical situation. We assessed psychopathology with the 21-item HAMD (Hamilton Depression Rating Scale)[10] and obtained blood samples within 5 days after hospital admission, and after 2 and 5 weeks of antidepressant treatment in 12 patients defined as remitters (HAMD <10 after 5 weeks; comparable to the HAMD-cutoff of 7 frequently used for 17-item HAMD Rating Scale;[11] mean: 4.8±2.2  s.d.) and 12 nonresponders (HAMD improvement <25% from baseline after 5 weeks compared with admission; mean: 18.9±4.3  s.d.). The patients were treated with different types of antidepressants. Nevertheless, there was no significant difference between the compared groups within the first 5 weeks of treatment. One hundred and fourty-two additional unipolar depressed inpatients (80 responders, 62 nonresponders after 5 weeks) of the Munich Antidepressant Response Signature study with available RNA sampling served as a replication and underwent the same assessment protocol as the discovery sample.

The study was approved by the Ethics Committees of the Medical Faculties at the Ludwig Maximilians University, Munich, Germany. Written informed consent was obtained from all subjects, and the study was carried out in accordance with the latest revision of the Declaration of Helsinki.

RNA preparation and gene expression analysis

Blood was drawn between 0800 and 0900 h; 10 h before blood collection, eating, drinking (except water), caffeine and nicotine intake were not allowed. Whole blood was collected in PAXgene RNA Blood tubes and isolated using PAXgene Blood RNA Kits (Qiagen, Hilden, Germany) according to the manufacturer's instructions. RNA was quantified by UV absorption using the Nanophotometer (Implen, Munich, Germany) and quality was determined on the Agilent Bioanalyzer (mean RIN=7.64±0.82 (s.d.); Agilent, Santa Clara, CA, USA). Two hundred nanograms of total RNA was reverse transcribed and converted to cDNA and subsequently biotin labeled using the Ambion TotalPrep Amplification Kit (Ambion, Austin, TX, USA). Seven hundred and fifty nanograms of cDNA was hybridized to Illumina HT-12 v3.0 arrays (Illumina, San Diego, California, USA) and incubated overnight at 55 °C according to the manufacturer's protocol. To reduce batch effects, complete probe sets of matched pairs were run in the same array slot. Arrays were washed, stained with Cy3 labeled streptavidin, dried and scanned on the Illumina BeadScan confocal laser scanner according to the manufacturer's protocol. The generated data were analyzed using the Illumina Beadstudio program and loaded into R for downstream statistical analysis.[12] Principal validity of the chip expression data was confirmed by technical replication of selected transcripts using quantitative real-time PCR (qPCR). QPCR assays were designed using the Roche universal probe library (http://qpcr.probefinder.com/organism.jsp).

For gene expression analysis in the replication sample, selected transcripts (Supplementary Table S1) were analyzed using qPCR. Total RNA was reverse transcribed to cDNA using random primers and the Superscript II reverse transcriptase (Invitrogen, Darmstadt, Germany) for qPCR run on the Roche 480 LightCycler (Roche Applied Science, Mannheim, Germany). All the samples were run in duplicates and duplicates discordant in crossing points by more than 0.4 cycles were excluded from the analysis. RORa, GCET2, CHI3L2, LSP1, SMARCC2, VPS13A, RECQL and RAP1GAP were chosen for RT-PCR validation from the microarray experiment. We applied the standard curve method for relative quantification using PLA, TBP, GUSB and HPRT as housekeeping genes reaching a high level of interclass correlation (interclass correlation coefficient range 0.789−0.996).

Combined dexamethasone/corticotropin-releasing hormone test

The combined dexamethasone suppression/corticotropin-releasing hormone challenge (dexamethasone/CRH) test was available in 80 patients at admission (6.9±2.4 (s.d.) days after admission) and in 39 patients after 6 weeks of treatment (46.9±3.3 (s.d.) days after admission). The test was performed as described in detail previously.[13] Briefly, the patient received 1.5 mg dexamethasone orally at 1800 h the day before the test, and the following day, adrenocorticotropic hormone (ACTH) and cortisol plasma levels were determined at 1500 h, 1530 h, 1545 h, 1600 h and 1615 h. An intravenous bolus of 100 μg of human CRH (Ferring, Kiel, Germany) was given at 1502 h. A radioimmunoassay kit was applied for the measurement of plasma cortisol concentrations (CT Cortisol RIA, DRG Diagnostics, Marburg, Germany). Plasma ACTH concentrations were assessed by an immunometric assay without extraction (cobas ECLIA, Roche Diagnostics, Rotkreuz, Switzerland). For the analysis of hormonal tests, we correlated RNA expression levels at admission and after 5 weeks with ACTH and cortisol at basal (that is, before CRH application) and after CRH application (that is, maximal concentration; Abas, Amax, Cbas, Cmax, respectively) using age and sex as covariates.

Animal experiments

Blood samples and brains were collected from adult male mice of the ‘stress reactivity' mouse model (details in S2).[14] To exclude phenotypic extremes, 16 mice of the ‘intermediate reactivity' line were used. The animals were killed under unstressed conditions and circadian trough levels of plasma corticosterone were determined by radioimmunoassay as described previously.[14] Hippocampi were quickly dissected and analyzed for RORa expression using qPCR as described before.[15, 16] At the time of sampling, all mice were about 4 months of age and single housed at least 2 weeks before the experiment. Animal husbandry was performed under standard laboratory conditions.

Genotyping for association analysis

Genome-wide genotypes were available from a different project (Illumina 100 K, Illumina 300 K, Illumina 610 K and Illumina OmniExpress chips).[17] All genotypes underwent a stringent quality control (removal of SNPs with callrates <98%, MAF <5%, deviation from HWE (P-value <10−5), removal of IDs with callrates <98%). Genotypes for rs12912233 were not available for all individuals and were imputed using impute v2 (http://mathgen.stats.ox.ac.uk/impute/impute_v2.html). The imputation info score was 0.97 indicating a good imputation quality.

Statistical analysis

Microarray expression data were transformed and normalized using variance stabilizing normalization.[18] For differential expression analysis, significantly regulated genes were ranked using an empirical Bayes method implemented in R.[19] A total of 48 770 annotated gene probes were detected. A total of 4826 genes were not expressed in any sample (no signal above background) and were excluded from the subsequent analysis. We performed an analysis of variance of contrasting gene expression of remitters and nonresponders at admission, after 2 and 5 weeks of treatment. We corrected for multiple testing using the false discovery rate (FDR).[20]

Group differences in demographic and baseline clinical data were compared using Pearson χ test in case of qualitative data and with two-sided t-tests (if not otherwise specified) for independent samples in case of quantitative data. GLM was applied for the analysis of gene expression in the replication experiment with response and remission at week 5 as dependent variables and age and sex as covariates. These statistical analyses were conducted with SPSS for Windows (version 18.0, SPSS, Chicago, IL, USA). PLINK (http://pngu.mgh.harvard.edu/~purcell/plink/) was used for genetic association studies.

Results

Clinical characteristics

Twelve remitting and 12 nonresponding male patients matched by age and baseline depression severity did not differ with respect to number of previous episodes, age of disease onset, treatment resistance (at admission), psychotic symptoms and body mass index as well as psychotropic and other medication during the study period (Table 1).

Table 1

Sample characteristics

	N	All (N=12)		Non-responder (N=12)		Remitter (N=12)		Value	df	Pa
		N/mean	(%/s.d.)	N/mean	(%/s.d.)	N/mean	(%/s.d.)
Baseline and disease characteristics
Age, mean (s.d.)	24	44.9	13.5	45.5	15.2	44.3	12.3	0.21	22	0.8382
Duration of hospital stay (weeks), mean (s.d.)	24	12.1	6.1	16.6	5.4	7.6	2.0	5.32	22	0.0000
Number of previous episodes, mean (s.d.)	21	4.4	7.4	7.0	9.6	1.6	1.6	1.75	19	0.0964
Age at disease onset, mean (s.d.)	24	29.8	15.2	28.3	17.0	31.3	13.7	−0.48	22	0.6389
HAMD at admission, mean (s.d.)	24	28.8	6.1	28.3	5.6	29.3	6.7	−0.36	22	0.7206
HAMD at discharge, mean (s.d.)	23	6.5	5.4	9.3	5.8	3.5	2.8	2.96	21	0.0074
HAMD at week 5, mean (s.d.)	22	12.5	8.0	18.9	4.3	4.8	2.2	9.43	20	0.0000
BMI at admission, mean (s.d.)	23	25.7	4.6	26.7	5.6	24.8	3.4	0.97	21	0.3441
HAMA at admission, mean (s.d.)	24	26.7	9.5	27.5	9.1	25.8	10.2	0.42	22	0.6775
Psychotic depression, No. (%)	24	3	12.5	1	8.3	2	16.7	0.38	1	>0.999
Treatment resistance, No. (%)	24	2	8.3	2	16.7	0	0.0	2.18	1	0.4780
Antidepressants (within 5 weeks), no. of patients (%)
Tricyclics	24	8	33.3	5	41.8	3	25.0	0.75	1	0.6670
SSRI	24	7	29.2	3	25.0	4	33.3	0.20	1	>0.999
SNRI	24	10	41.7	7	58.3	3	25.0	2.74	1	0.2140
NaSSA (=Mirtazapine)	24	9	37.5	4	33.3	5	41.7	0.18	1	>0.999
NARI	24	0		0		0
MAOI	24	1	4.2	0	0.0	1	8.3	1.04	1	>0.999
Others
Bupropion	24	4	16.7	4	33.3	0	0.0	4.80	1	0.0930
Trazodon	24	3	12.5	3	25.0	0	0.0	3.43	1	0.2170
Buspiron	24	1	4.2	1	8.3	0	0.0	1.04	1	>0.999
Psychotropics, others (within 5 weeks), no. of patients (%)
Antipsychotics	24	15	62.5	9	75.0	6	50.0	1.60	1	0.4000
Mood stabilizers	24	14	58.3	6	50.0	8	66.7	0.69	1	0.4080
Lithium	24	3	12.5	2	16.7	1	8.3	0.38	1	>0.999
Benzodiazepines	24	13	54.2	7	58.3	6	50.0	0.17	1	0.6820
Internal medication (within 5 weeks), no. of patients (%)
Thyroid medication	24	4	16.7	3	25.0	1	8.3	1.20	1	0.5900
Analgetics	24	6	25.0	2	16.7	4	33.3	0.89	1	0.6400
Cardiovascular medication	24	8	33.3	6	50.0	2	16.7	3.00	1	0.1930
Antibiotics	24	1	4.2	0	0.0	1	8.3	1.04	1	>0.999

Abbreviations: BMI, body mass index; HAMA, Hamilton Anxiety Rating Scale; HAMD, Hamilton Depression Rating Scale.

For some variables, there are missing data and N does not equal number of total patients. Percentages are based on available data. Bold values signify P<0.05.

Two-sided t-test and Chi-square, respectively, comparing nonresponders and remitters; Fisher's exact test used when cell count was <5.

Whole-genome gene expression

Using a combined contrast analysis of variance (comparing remitters and nonresponders at admission, after 2 and after 5 weeks), we found 127 genes having adjusted (FDR) P-values <0.05 (Supplementary Table S3). The lowest P-values were reached for transcripts of the Rap1 GTPase-activating protein (RAP1GAP, P=2.14 × 10−8, PFDR=9.41 × 10−4), germinal center expressed transcript 2 (GCET2, P=6.73 × 10−8, PFDR=1.13 × 10−3) and chitinase 3-like protein 2 (CHI3L2, P=7.69 × 10−8, PFDR=1.13 × 10−3) with lower expression levels in remitters compared with nonresponders in each case.

Comparing gene expression after 2 or 5 weeks of treatment with admission within remitters and nonresponders, respectively, lowest P-values were observed for SWI/SNF-related, matrix-associated, actin-dependent regulator of chromatin, subfamily c, member 2 (SMARCC2, P=5.91 × 10−8, PFDR=2.42 × 10−3) and leukocyte-specific protein 1 (LSP1, P=1.49 × 10−7, PFDR=2.42 × 10−3; Supplementary Table S4). RORa recently implicated in depression and posttraumatic stress disorder[21, 22] ranked 12th in the combined contrast (P=2.48 × 10−6, PFDR=9.06 × 10−3) and showed up in additional single contrasts (Supplementary Table S4).

Replication of selected transcripts in 142 patients

One hundred and forty-two additional unipolar depressed inpatients served as a replication sample. Eighty responders and 62 nonresponders after 5 weeks did not differ in baseline characteristics including age, sex, baseline depression severity, age of disease onset, body mass index and psychotic symptoms (Supplementary Table S5). For qPCR, we included the most significant associated transcripts of the whole-genome expression experiment of the combined (RAP1GAP, GCET2, CHIi3L2) and single contrasts (SMARCC2, LSP1). In the replication experiment, we combined this hypothesis-free replication approach with a candidate-gene approach based on biological plausibility: we included RORa, because it was the highest associated gene in the combined contrast rank list with a previously described association with depression and stress-related disorders.[21, 22] To proof the robustness of the chip experiment, we randomly selected two further transcripts with sufficient average expression levels and fold changes from the top 120 hits, vacuolar protein sorting 13 homolog A (VPS13A) and RecQ protein-like (DNA helicase Q1-like; RECQL), to be included in the replication experiment.

Using a GLM with age and sex as covariates, we confirmed association of RORa (P=6.23 × 10−4), GCET2 (P=2.08 × 10−3) and CHI3L2 (P=4.45 × 10−2) with lower expression levels in responders compared with nonresponders at admission in each case, which is the same direction of effect as in the chip experiment (Table 2). In addition, LSP1 significantly decreased after 5 weeks of treatment in responders and remitters (P=2.91 × 10−2, P=4.71 × 10−3, respectively). Comparing responders and nonresponders, expression of RORa was lower in responders also after 5 weeks of treatment (P=3.29 × 10−3), but did not reach significance after 2 weeks (P=1.18 × 10−2; Figure 1).

Table 2

Replication sample (N=142)

		Nonrespondersa		Respondersa		Pb	Non-remittersa		Remittersa		Pb
		Mean	s.d.	Mean	s.d.		Mean	s.d.	Mean	s.d.
RORa	Admission	1.10	0.35	0.93	0.25	0.0006	1.05	0.32	0.92	0.26	0.0105
	After 2 weeks (% change)	19.79	33.92	25.42	45.79	0.5163	20.24	36.84	28.93	48.94	0.3127
	After 5 weeks (% change)	3.25	38.95	4.66	41.11	0.9689	1.03	35.13	9.87	48.03	0.2749
GCET2	Admission	0.95	0.32	0.78	0.28	0.0021	0.90	0.32	0.76	0.26	0.0162
	After 2 weeks (% change)	12.35	30.63	23.14	54.58	0.2825	13.89	34.47	28.45	63.62	0.1409
	After 5 weeks (% change)	1.24	34.10	6.68	41.04	0.5768	1.11	32.79	10.44	46.50	0.2536
CHI3L2	Admission	1.33	0.74	1.11	0.62	0.0446	1.35	0.76	1.18	0.67	0.0968
	After 2 weeks (% change)	13.73	33.00	11.71	32.54	0.7466	11.24	34.77	15.41	27.69	0.5049
	After 5 weeks (% change)	10.58	30.38	14.90	37.99	0.6186	10.18	31.86	18.45	39.67	0.2471
LSP1	Admission	0.76	0.39	0.82	0.27	0.2631	0.75	0.33	0.86	0.30	0.0427
	After 2 weeks (% change)	10.31	33.27	15.97	42.13	0.2446	15.67	38.53	9.04	38.72	0.5172
	After 5 weeks (% change)	12.98	41.09	−1.95	34.48	0.0291	11.27	40.47	−8.28	29.42	0.0047
SMARCC2	Admission	0.87	0.32	0.90	0.20	0.5490	0.86	0.27	0.92	0.23	0.2106
	After 2 weeks (% change)	7.59	25.73	9.56	27.22	0.5133	9.67	24.90	6.68	29.86	0.6903
	After 5 weeks (% change)	2.66	29.74	−2.25	21.38	0.3179	2.19	27.06	−4.52	21.48	0.1619
VPS13A	Admission	1.09	0.22	1.02	0.21	0.0919	1.07	0.21	1.00	0.22	0.0636
	After 2 weeks (% change)	5.94	20.96	8.11	32.24	0.7332	4.30	22.63	13.32	36.31	0.1094
	After 5 weeks (% change)	0.42	26.24	2.08	28.96	0.8927	-0.80	24.51	5.51	32.93	0.2800
RECQL	Admission	1.25	0.33	1.18	0.30	0.1336	1.24	0.32	1.17	0.30	0.1725
	After 2 weeks (% change)	13.94	36.48	17.72	37.76	0.5302	13.84	34.70	20.93	41.87	0.3016
	After 5 weeks (% change)	3.03	33.97	4.59	36.14	0.3016	2.98	34.56	5.69	36.39	0.8388
RAP1GAP	Admission	1.12	2.06	1.56	7.54	0.7471	0.92	1.74	2.24	9.70	0.2255
	After 2 weeks (% change)	115.49	220.90	153.86	418.30	0.3219	156.49	353.31	97.11	335.46	0.5409
	After 5 weeks (% change)	157.89	314.37	86.77	165.87	0.1227	148.46	273.69	59.11	160.18	0.0596

Abbreviations: RORa, retinoid-related orphan receptor alpha; GCET, germinal center expressed transcript 2; CHl3L2, chitinase 3-like protein 2; LSP1, leukocyte-specific protein 1; SMARCC2, SWI/SNF-related, matrix-associated, actin-dependent regulator of chromatin, subfamily c, member 2; VPS13A, vacuolar protein sorting 13 homolog A; RECQL, RecQ protein-like; RAP1GAP, Rap1 GTPase-activating protein. Bold values signify P<0.05.

After 5 weeks.

GLM with age and sex as covariates.

Figure 1

Relative RORa expression in antidepressant responders and nonresponders. Note that RORa expression appears to be elevated in nonresponders throughout the observation period reaching statistical significance at admission and after 5 weeks (P=6.23 × 10−4, P=3.29 × 10−3, respectively; GLM with age and sex as covariates). ** P-value <0.01; *** P-value <0.001. Error bars are s.e.m. RORa, retinoid-related orphan receptor alpha.

Applying the conservative Bonferroni-based correction for multiple testing in our replication experiment (α=6.25 × 10−3), RORa and GCET2 remained significantly different at baseline, and LSP1 after 5 weeks. Given the suspected impact of RORa in psychiatric disorders, we focused on this gene in our further validation experiments.

RORa expression, sleep, diurnal depressive symptoms and hypothalamic–pituitary–adrenal axis

As RORa is a so-called clock gene involved in circadian rhythm generation, we further tested whether the expression of RORa was related to presence of sleep disturbances or diurnal variations of depressive symptoms at admission using the respective items on the HAMD scale (4, 5, 6 and 18, respectively). We could not find a correlation of these variables with RORa expression (P=5.56 × 10−2, P=8.35 × 10−2, P=3.46 × 10−2, P=7.35 × 10−2, respectively; partial correlation analysis with age and sex as covariates).

As functional RORA knockout mice show an altered ACTH and corticosterone response to stress,[23] we further tested the hypothesis whether RORa expression has any impact on the regulation of the hypothalamic–pituitary–adrenal-axis in the combined dexamethasone/CRH test. We found significant negative correlations (after controlling for age and sex) with RORa expression for Amax (−231; P=4.49 × 10−2; N=80) at admission and for Cbas (−0.366; P=2.58 × 10−2; N=39) after 6 weeks (Figure 2). No significant correlation was found with HAMD values at the time of testing (that is, at admission and at week 6; data not shown).

Figure 2

ACTH (a) and cortisol (b) response to the dex/CRH test at admission (left) and after 6 weeks (right), depending on the RORa expression level (upper vs lower 50% of expression values). Error bars are s.e.m. ACTH, adrenocorticotropic hormone; dex/CRH, dexamethasone/corticotropin-releasing hormone; RORa, retinoid-related orphan receptor alpha.

Further investigating the negative correlation of RORa and hypothalamic–pituitary–adrenal axis activity observed in our patients, we quantified the expression of RORa in the hippocampus of mice. Testing for a correlation of hippocampal RORa expression in intermediate reactivity animals of the ‘stress reactivity' mouse model[14] with plasma corticosterone concentrations under basal conditions, we indeed detected a significant negative correlation (r=−0.534; P=3.31 × 10−2; N=16).

Association of RORa mRNA with RORa polymorphisms (SNPs)

We tested the association of RORa expression with two recently described RORa SNPs implicated in the traits depression (rs12912233)[21] and antidepressant response (rs809736),[24] respectively. Using linear regression, we found a nominal association of RORa mRNA expression change after 5 weeks with the rs809736 genotype (β=11.98±5.61 (s.e.m.); P=3.46 × 10−2; genotype distribution: 85 (AA), 50 (AG) and 7 (GG); with the G allele being associated with lower RORa mRNA expression), but not with rs12912233 (β=−8.5160±5.15 (s.e.m.); P=1.01 × 10−1).

Interrogating the publicly available pharmacogenetic data set of the PCG project (http://www.broadinstitute.org/mpg/ricopili/; PhaCoGe Quan12 (Feb. 2013)), no association of RORa rs809736 with antidepressant response was found (OR=1.01; P=1.53 × 10−1), whereas data for rs12912233 were not available.

Discussion

We identified several gene transcripts emerging from an unbiased whole-genome expression approach that were differentially expressed in antidepressant remitters and nonresponders, and that replicated in a representative clinical sample.

Using whole-genome expression arrays for the identification of antidepressant predictors so far produced conflicting results.[4, 5, 6, 7] The fact that we were able to replicate our findings from a relatively small discovery sample may be related to the carefully selected and severity-matched clinical extremes of remitters and nonresponders in the microarray experiment.[8] Three of the five top ranked transcripts from the microarray experiment (GCET2, CHI3L2 and LSP1) could be replicated in a sample of 142 patients with the same direction of effect. CHI3L2 has recently been reported to activate the MAPK pathway through phosphorylation of ERK1/2 in glial tumor cells.[25, 26] GCET2 and LSP1 are expressed in leukocytes and have a role as markers in B-cell lymphoma (GCET2)[27] and cellular transmigration (LSP1),[28, 29] respectively. Interestingly, in this context, long-term treatment (⩾12 weeks) with selective serotonin-reuptake inhibitors modulates B-cell proliferation.[30, 31]

RORa, the most significantly associated transcript in our replication study, has very recently gained attention as a new candidate in stress-related disorders, especially depression. The clock gene RORa is a transcription factor belonging to the steroid hormone receptor superfamily and has been implicated in response to cellular stress, cellular differentiation and proliferation, as well as circadian rhythm generation in different tissues.[32, 33] RORa-responsive genes include genes involved in calcium second messenger signaling and glutamatergic signaling pathways,[33] as well as the control of metabolic pathways.[34, 35] RORa is rhythmically expressed via E-box-binding of CLOCK-BMAL heterodimers. RORa then drives Bmal1 expression as a part of the so-called stabilizing loop of the molecular clock of the hypothalamic suprachiasmatic nucleus.[36, 37]

For a long time, it has been speculated that an altered molecular clock is involved in the pathophysiology of depression as many patients show diurnal variations of depressive symptoms and disrupted sleep.[38, 39] Lithium, used as a treatment of mania and as prevention of bipolar disorder, inhibits glycogen synthase kinase 3b that, in turn, regulates protein stability of components of the molecular clock.[40, 41] Under treatment with lithium, increased amplitudes of clock gene oscillations have been observed on the cellular and tissue level,[42] as well as increased circadian period length on the behavioral level.[42, 43] Intriguingly, McCarthy et al.[44] showed in fibroblasts of bipolar patients and healthy controls that circadian period lengthening induced by lithium is dependent on the RORa rs1292233 genotype.

Rs12912233 was recently found to be associated with trait depression in a genome-wide association study in a Caucasian epidemiologic sample.[21] Another intronic RORa SNP (rs809736) ranked third in the genome-wide association study with antidepressant response to citalopram in the STAR*D sample.[24] RORa polymorphisms have further been implicated in posttraumatic stress disorder[22] and autism spectrum disorder.[45, 46] We found rs809736 to be nominally associated with RORa expression in the present sample with the same direction of effect, that is, the allele associated with response in STAR*D was associated with lower RORa expression in the present study. Nevertheless, rs809736 was not replicated in the large PGC pharmacogenetic sample, hence, the nominal significant association of rs809736 with RORa gene expression in our cohort needs further replication.

Functional RORa knockout Staggerer mice (RORa sg/sg) show an enhanced ACTH and corticosterone response to novelty stress and a flattened corticosterone circadian rhythm.[23] Consistent with this observation, we found in our patients a negative correlation of RORa expression with basal and stimulated stress hormones in the dexamethasone/CRH test, an effect that was independent of the depressive states on the hypothalamic–pituitary–adrenal-axis. As a validation of this observation and extension to central pathophysiology, we found a negative correlation of hippocampal RORa expression with plasma corticosterone levels in intermediate reactivity mice. Nevertheless, RORa expression in our patients was lower in responders compared with nonresponders in the microarray as well as in the replication study which is partly counter-intuitive regarding the negative correlation with the hypothalamic–pituitary–adrenal-axis.[2, 47] The physiological meaning of this finding in depression remains unclear, but may point to adaptive changes in nonresponders. The idea of adaptive processes may be indirectly supported by the recent finding of Ikeda et al.[48] who showed that RORa upregulates expression of dopamine D3 receptor (DRD3) in mice ventral striatum. This interpretation is in line with restored dopaminergic signaling that has been implicated in recovery from depression.[49]

The present study harbors several limitations. The study did not follow a strict hypothesis-free design. Although we chose an unbiased approach in the discovery sample, not only were the most significant hits included in the replication experiment, but also the ‘promising candidate' RORa that was then subject for further validation in functional experiments as it turned out to be the most significantly associated transcript in the replication study. Thus, our results can only be interpreted cautiously, and a further replication of RORa is crucial, as is the validation of the other transcripts, GCET2, CHI3L2 and LSP1.

RNA expression analysis is challenging in numerous methodological aspects that can significantly impact on quality and reproducibility of our data.[8] Although we cannot exclude environmental factors having affected gene expression more than antidepressant response in our patients, we undertook several measures to reduce biological variability accompanied with this method, for example, standardized pre-analytic conditions, consideration of batch effects, technical replication with qPCR and the use of four housekeeping genes for normalization of expression values yielding a high interclass reliability. By carefully selecting our patients, we further tried to reduce biological variability emerging from the subjects themselves. This approach allowed us to observe slight, but robust expression differences of RORa that were comparable to gene expression changes observed in similar studies[50] and independent of effects of sleep disturbances or circadian symptom patterns that have recently been reported to affect RORa expression in a transcriptomics study.[51] Further, although gender was included as a covariate in each analysis, the effect observed for RORa may be primarily mediated by male patients as the discovery sample included males only. Indeed, the association of RORa expression (at admission) with response was stronger when analyzing males only (N=75; P=3.17 × 10−4 vs P=6.23 × 10−4), but lost significance in women (N=67; P=2.19 × 10−2). This observation may be related to the reduced statistical power in this subanalysis, but needs further attention in future studies.

The current study does not have a placebo arm. Thus placebo-related and real pharmacological effects cannot be separated. Different dosages between responders and nonresponders may have affected gene expression levels. Nevertheless, analyzing our data in this respect for the most frequently prescribed antidepressants in our study (venlafaxine; escitalopram; citalopram; mirtazapine; amitriptyline, including amitriptylinoxide) at week 5, we could not find a significant difference in prescribed dosages between responders and nonresponders (data not shown). In addition, only a part of our patients (N=101) was treated with only one class of antidepressant agent and medication switches within the 5-week observation period have not been analyzed specifically. On the other hand, all available antidepressants act as monoaminergic enhancers and effects of antidepressants are not specific, particularly after prolonged administration. The current study design was well suited to identify possible new targets that make some patients respond more easily to antidepressant drugs than others, independent of the drug given. Due to the naturalistic design of the present study, our results are possibly more generalizable and more relevant for clinical practice because most patients are treated with more than one psychotropic agent drug.[2, 52, 53]

In summary, in this study, we identified several transcripts that have the potential to serve as biomarkers for antidepressant response. Among these biomarker candidates, the circadian gene transcript RORa may deserve special attention as a candidate conveying disposition for depression and as a predictor of antidepressant response.