Introduction
Major depressive disorder (MDD) is a highly recurrent and prevalent disorder with societal costs. Despite the increased amount of available pharmacological treatments, many patients do not respond with the first-line treatment, while others experience costly delays until they find which treatment is efficacious for them. Whereas it is estimated that about 50% of patients with uncomplicated MDD will respond to any single antidepressant (Papakostas & Fava, Reference Papakostas and Fava2008, Reference Papakostas and Fava2009), other reports show lower response rates (Souery et al. Reference Souery, Oswald, Massat, Bailer, Bollen, Demyttenaere, Kasper, Lecrubier, Montgomery, Serretti, Zohar and Mendlewicz2007). For example, in a large, multicentre trial, the Sequenced Treatment Alternatives to Relieve Depression (STAR*D), approximately one in three patients with MDD achieved remission after being treated with the selective serotonin reuptake inhibitor (SSRI) citalopram (Trivedi et al. Reference Trivedi, Rush, Wisniewski, Nierenberg, Warden, Ritz, Norquist, Howland, Lebowitz, McGrath, Shores-Wilson, Biggs, Balasubramani and Fava2006a ).
There is substantial evidence that antidepressant response is to a certain extent heritable (O'Reilly et al. Reference O'Reilly, Bogue and Singh1994; Franchini et al. Reference Franchini, Serretti, Gasperini and Smeraldi1998). Pharmacogenetic research seems to be a promising path toward personalized treatment; however, the link between genetic risk variants and response to treatment is far from direct. To date, genomewide association studies (GWASs) on treatment response of antidepressants (Ising et al. Reference Ising, Lucae, Binder, Bettecken, Uhr, Ripke, Kohli, Hennings, Horstmann, Kloiber, Menke, Bondy, Rupprecht, Domschke, Baune, Arolt, Rush, Holsboer and Muller-Myhsok2009; Garriock et al. Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010; Uher et al. Reference Uher, Perroud, Ng, Hauser, Henigsberg, Maier, Mors, Placentino, Rietschel, Souery, Zagar, Czerski, Jerman, Larsen, Schulze, Zobel, Cohen-Woods, Pirlo, Butler, Muglia, Barnes, Lathrop, Farmer, Breen, Aitchison, Craig, Lewis and McGuffin2010) provide modest results, which do not reach genomewide significance and are hard to replicate, thereby making any translation into potential clinical utility difficult. On the other hand, there is a substantial research showing that candidate polymorphisms of target proteins may be involved in antidepressant drug action (Serretti et al. Reference Serretti, Benedetti, Zanardi and Smeraldi2005; Kato & Serretti, Reference Kato and Serretti2010; Porcelli et al. Reference Porcelli, Drago, Fabbri, Gibiino, Calati and Serretti2011). Environmental factors, both early and late in life, also influence treatment response (Uher, Reference Uher2011). It has been recently observed that the same genes, environments, and gene–environment (G × E) interactions implicated in aetiological pathways to depression may also play a role in predicting treatment response (Keers & Uher, Reference Keers and Uher2012). Genetic variation, as on the serotonin-transporter-linked polymorphic region (5-HTTLPR) for example, may increase the risk for MDD (Uher & McGuffin, Reference Uher and McGuffin2010) or may lead to a worse response to the SSRI treatment (Mandelli et al. Reference Mandelli, Marino, Pirovano, Calati, Zanardi, Colombo and Serretti2009; Keers et al. Reference Keers, Uher, Huezo-Diaz, Smith, Jaffee, Rietschel, Henigsberg, Kozel, Mors, Maier, Zobel, Hauser, Souery, Placentino, Larsen, Dmitrzak-Weglarz, Gupta, Hoda, Craig, McGuffin, Farmer and Aitchison2011) under unfavourable psychosocial environmental conditions. Although there is accumulating evidence of G × E interactions on the aetiology of depression emerging from candidate gene studies, statistical and methodological challenges hinder GWASs to examine the potential moderating role of life stressors (Khoury & Wacholder, Reference Khoury and Wacholder2009). However, considering the complexity of the phenotype, it has been suggested that alternative pharmacogenetic study designs should be considered, in particular those that take into account the manifold sociodemographic and environmental variables that make an impact on treatment response in the clinical setting (Malhotra, Reference Malhotra2010). Such factors may interact with genetic variation to influence antidepressant response. For example, the patient's quality of life (QoL) has been previously found to be an independent predictor of acute antidepressant treatment response (Pyne et al. Reference Pyne, Bullock, Kaplan, Smith, Gillin, Golshan, Kelsoe and Williams2001) and low QoL was associated with lower remission scores in the STAR*D trial (Trivedi et al. Reference Trivedi, Rush, Wisniewski, Nierenberg, Warden, Ritz, Norquist, Howland, Lebowitz, McGrath, Shores-Wilson, Biggs, Balasubramani and Fava2006a ).
The aim of the present study was twofold. First, in an attempt to improve the predictive model of treatment response, we examined the possible moderating role of the patient's level of QoL, in a GWAS model of treatment response to citalopram, in a large group of depressed patients from the STAR*D trial. Other sociodemographic and clinical characteristics, such as age, sex and symptom severity, known to influence response to treatment (Drago & Serretti, Reference Drago and Serretti2011), were also included in the model. Second, since a number of candidate genes have emerged from prior G × E interaction studies on depression, we were interested in investigating whether this specific and possibly putative cluster of genes is likely to represent a causal pathway in predicting treatment response. We therefore constructed an enrichment analysis to examine whether those genes are associated with antidepressant treatment response.
Method
Participants
The STAR*D trial has been previously described in detail elsewhere (Fava et al. Reference Fava, Rush, Trivedi, Nierenberg, Thase, Sackeim, Quitkin, Wisniewski, Lavori, Rosenbaum and Kupfer2003; Rush et al. Reference Rush, Trivedi and Fava2003). In brief, we focused on data collected at level 1 of the trial, wherein all participants were treated with citalopram. Genotypic and phenotypic data are available from the National Institute of Mental Health (NIMH) Center for Collaborative Genetic Studies on Mental Disorders.
Genotypic information was available for 1939 STAR*D samples. Genotyping was performed for 964 samples from the STAR*D discovery project using the Human Mapping 500 K Array set (Affymetrix Inc., USA). An additional 975 samples was genotyped using the Genome-Wide Human SNP Array 5.0 (Affymetrix Inc., USA). This sample and genotyping techniques have been described in detail by Garriock et al. (Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010). The two groups were balanced in terms of variables such as ethnicity, sex and responders and non-responders (Garriock et al. Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010).
Quality control was performed using PLINK v1.07 (Purcell et al. Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007). A total of 496 415 markers were initially included and the following single nucleotide polymorphisms (SNPs) were subsequently removed: 10 144 SNPs with minor allele frequencies (MAFs) < 1%; 44 483 SNPs that failed the missingness per marker test (>0.05); 23 121 SNPs on the basis of the Hardy–Weinberg equilibrium test (10 × 10−5); and 74 SNPs that existed in duplicate. The remaining 418 865 SNPs were used for further analyses.
Population heterogeneity
The –genome function was used in PLINK and independence was computed by the –indep function on the basis of 50 SNPs in the window, five SNP window shift, and a variance inflation factor of 2. A total of 140 313 independent SNPs were yielded after frequency and genotype pruning. Clusters were computed using a forced entry of 10 clusters using the k function in PLINK. This was done to simulate the 10 clusters used by Garriock et al. (Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010) and as previously suggested (Price et al. Reference Price, Patterson, Plenge, Weinblatt, Shadick and Reich2006). Quantile–quantile (Q-Q) plots of expected and observed p values were constructed for the phenotype and adjustment for clusters significantly improved inflation (λ = 1.03; for Q-Q plot, see online Supplementary Fig. S1).
Phenotype
The phenotype was constructed on the basis of the Garriock et al. (Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010) study. In brief, the responders group was comprised of participants that showed a 50% or greater reduction, from baseline to final visit, in the self-report 16-item Quick Inventory of Depressive Symptomatology (QIDS-SR). Participants who did not meet this criterion were non-responders. Only participants who completed a 6-week window of treatment with citalopram were included. Due to different versions of the data files in the STAR*D, in our study there were available phenotypic data for 1947 participants; 969 showed non-response, whereas 978 responded to antidepressant treatment. Participants who belonged to ethnic groups that comprised only a small proportion of the sample were excluded: Asian ancestry (n = 21), Pacific Island ancestry (n = 10), Native American (n = 16) or mixed (n = 68) or unknown (n = 1). Subsequently, there were 1279 non-Hispanic Caucasians, 305 African-Americans and 247 Hispanic Caucasians available for analyses. Of those 18 31, 189 had not completed 42 days of treatment and were excluded. Another 88 participants were removed because they had more than 5% of their genotypes missing; 30 participants were shown to be cryptically related to others in the sample [assessed using an identity by descent (IBD) > 0.1875 in PLINK; Anderson et al. Reference Anderson, Pettersson, Clarke, Cardon, Morris and Zondervan2010]. Of the remaining participants, 40 participants were coded with ambiguous sex (genetic–clinical sex disagreement assessed in PLINK) and therefore excluded. Covariate data, including the QoL measurement, were not available for 64 participants. This yielded a final sample size of 1426 participants for analyses on the response phenotype.
Assessment of QoL and covariates
The Quality of Life Enjoyment and Satisfaction Questionnaire was used to have an estimate of a person's functioning and QoL upon entering the trial (Daly et al. Reference Daly, Trivedi, Wisniewski, Nierenberg, Gaynes, Warden, Morris, Luther, Farabaugh, Cook and Rush2010). The measure was administered using an interactive voice response telephone interview, within 72 h of the baseline visit (Daly et al. Reference Daly, Trivedi, Wisniewski, Nierenberg, Gaynes, Warden, Morris, Luther, Farabaugh, Cook and Rush2010). It is a self-report 16-item questionnaire, which measures specific domains such as quality of social relationships, living or housing situation and physical health (Endicott et al. Reference Endicott, Nee, Harrison and Blumenthal1993). Each item is scored on a five-point Likert scale (1 = very poor to 5 = very good) indicating the degree of satisfaction achieved during the past week.
Covariates
We included age and sex as covariates since sociodemographic variables were associated with antidepressant response in the STAR*D (Drago & Serretti, Reference Drago and Serretti2011) and it has been found that sex mediates the effects of experimental manipulation of serotoninergic levels on mood (Booij et al. Reference Booij, Van der Does, Benkelfat, Bremner, Cowen, Fava, Gillin, Leyton, Moore, Smith and Van der Kloot2002). Patients' symptom severity was previously associated with QoL levels, in the STAR*D trial (Trivedi et al. Reference Trivedi, Rush, Wisniewski, Warden, McKinney, Downing, Berman, Farabaugh, Luther, Nierenberg, Callan and Sackeim2006b ), so we included the severity of depressive symptoms at baseline as another covariate (continuous variable). We also controlled for genetic stratification (population clusters).
Statistical analysis
Part A, GWASs
-
(i) A genomewide interaction model was examined in PLINK, using logistic regression, including the –interaction, and –parameters functions, with QoL as the interaction factor and treatment response as the outcome, while controlling for covariates sex, age, depression severity and population clusters. Interaction was assessed using the additive model, and QoL was binary (median split: high > 43 v. low ⩽43). The log odds ratio (OR) (β Gadd × QoL) of the interaction effect was the main outcome. This model was adopted from Cornelis et al. (Reference Cornelis, Tchetgen, Liang, Qi, Chatterjee, Hu and Kraft2012) who suggested that dichotomizing the interaction variable results in a model that is less subject to inflation (λ = about 1.01). Although some information may be lost due to this dichotomization, this saturated model may eliminate some mis-specification under the null hypothesis and decrease type I error (Cornelis et al. Reference Cornelis, Tchetgen, Liang, Qi, Chatterjee, Hu and Kraft2012). We also run an exploratory analysis with response as a continuous trait (percentage change on the QIDS from baseline).
-
(ii) A genomewide logistic regression model on treatment response, including all covariates (and QoL) and omitting the interaction term, was tested. This analysis was conducted in order to control for a possible lack in reliability due to the inclusion of the interaction term in the first model (Purcell et al. Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007). The WGAViewer program (http://compute1.lsrc.duke.edu/softwares/WGAViewer/) was used for the inspection of the GWAS results (Ge et al. Reference Ge, Zhang, Need, Martin, Fellay, Urban, Telenti and Goldstein2008).
Part B, enrichment analysis
Genes were selected on the basis of their prior implication in G × E interaction studies. G × E studies were searched via published articles and reviews in the National Center for Biotechnology Information (NCBI) PubMed database. Screening of the resulting studies resulted in a list of 15 genes that have been previously examined together with environmental factors in predicting depression status or severity, or antidepressant response. The list constructed for the enrichment analysis involved the following 15 genes: SLC6A4 (serotonin transporter); TPH1 (tryptophan hydroxylase 1); TPH2 (tryptophan hydroxylase 2); MAOA (monoamine oxidase A); HTR1A (serotonin receptor 1A); HTR2A (serotonin receptor 2A); HTR2C (serotonin receptor 2C); CRHR1 (corticotropin-releasing hormone receptor 1); NR3C1 (glucocorticoid receptor subfamily 3, group C, member 1); FKBP5 (FK506 binding protein 5); BDNF (brain-derived neurotrophic factor); SLC6A3 (dopamine transporter); DRD2 (dopamine receptor D2); COMT (catechol-O-methyltransferase); and GABRA2 (GABA A receptor α 2). The enrichment analysis was conducted using the following steps. The genomic sections that contained the genes were identified (PLINK annotation); these sections were imputed (using the CEU HapMap 1000 genomes), checked for quality (imputation quality score > 0.9) and pruned (r 2 > 0.2); SNPs with low MAF were excluded. Then, the association between the phenotype and the list of SNPs (original and imputed) harboured by those genes was tested. From the remaining part of the genome a random selection of variations was taken 10−5 times (the number of total SNPs being the same as in the ‘gene list of interest’, and all of them belonging to genes randomly selected throughout the genome). The association of this random list of variations with the phenotype was tested. Then, the number of variations significantly (p < 0.05) associated with the outcome in the ‘gene list of interest’ and the ‘random list’ were tested for a significant differential distribution (Fisher's exact test). The permuted p was then extracted from the 105 Fisher tests, as suggested previously (Phipson & Smyth, Reference Phipson and Smyth2010). Subsequently, the GeneMANIA program (University of Toronto; http://www.genemania.org/) was used to produce networks of the top pathway list genes on the basis of their genetic interaction, physical interaction, co-expression or co-localization with the identified proteins, as used by other investigators (Siddiqui et al. Reference Siddiqui, Li, Luo, Karaiskakis, Hou, Kislinger, Westwood, Morris and Lipshitz2012).
Results
The GWAS analysis involved 418 865 SNPs for association with treatment response and moderation by QoL. Of the 1426 participants included in the analysis, 774 participants were responders and 652 were non-responders to treatment. Sociodemographic characteristics are shown in Table 1. Responders had a lower mean age compared with non-responders. There were no significant sex differences on the response outcome, although it can be observed that males were slightly over-represented in the non-responder group (χ 2 1 = 3.47, p = 0.06). QoL scores were normally distributed in the sample. In the whole sample, QoL scores correlated with severity of depression: the higher the symptoms, the lower the QoL (Pearson's r = −0.54, p < 0.001). Responders were more likely to report a higher QoL compared with non-responders (χ 2 1 = 33.26, p < 0.001). Responders also had lower symptom severity upon entering the trial (t 1424 = 4.51, p < 0.001). There were also differences in ethnicity between responders and non-responders (χ 2 2 = 23.92, p < 0.001); African-Americans were over-represented in the non-responder group compared with the responder group.
Sociodemographic characteristics of the sample (n = 1426)
s.d., Standard deviation; QIDS, Quick Inventory of Depressive Symptomatology; QoL, quality of life.
* p < 0.001.
Genomewide association analysis
Genomewide interaction
The results of this model are presented in Table 2, where p values of the interaction effect are reported. We found one SNP (rs520210) within the NEDD4L (neural precursor cell expressed, developmentally down-regulated 4-like E3) gene that was associated with treatment response. Responders were most likely to be carriers of the minor allele (A), compared with non-responders, especially those reporting low QoL upon entering the trial. Visual inspection of the results showed that for all top 10 resulting SNPs, responders had higher frequencies of the minor allele, compared with non-responders. Non-responders with high QoL were more likely to have low minor allele frequencies for the top SNPs (online Supplementary Fig. S2). The Manhattan plot of this analysis is reported in online Supplementary Fig. S3. For purposes of ethnic specificity, and since African-Americans were over-represented in the non-responder group, we repeated the model including only Caucasian individuals. Of the 1426 patients, 238 were excluded, who were of African-American ancestry; this group was too small for a separate GWAS analysis. The results in the Caucasian sample (n = 1188; 171 Hispanics included) correlate with those of the whole sample; however, lower significance values were noted for the top SNP rs520210 (p = 1.87 × 10−6). The rs520210 (NEDD4L), rs11775176 (tumour necrosis factor receptor superfamily, member 10b; TNFRSF10B), rs156466 (intergenic: LOC100505738; SEMA5A), rs408465 (SLC37A1; PDE9A) and rs11701162 (SLC37A1; PDE9A) were consistently among the top 20 SNPs in both analyses (Table 2; online Supplementary Table S1). Running the model in Caucasians with Hispanics excluded (n = 1017), we observed that the rs12427491 (LRCH1), rs156466 (intergenic: LOC100505738; SEMA5A) and rs408465 (SLC37A1; PDE9A) variants were the only ones resulting among the top 15 (online Supplementary Table S2), as in the original model with the whole sample. The rs520210 SNP did not reach genomewide significance (p = 6.98 × 10−5), but the effect was in the same direction with an OR of 0.79 for minor allele carriers. Finally, in an exploratory analysis we investigated response to treatment by using the continuous quantitative trait as an outcome, namely the percentage change in QIDS from baseline. The rs520210 SNP did not reach genomewide significance (p = 1.35 × 10−5) and it did not result among the top 15 SNPs, but the effect remained in the same direction (online Supplementary Table S3). The top SNPs of the latter analysis did not correlate with the ones with response using a dichotomous outcome.
Top 15 SNPs of the genomewide interaction model on citalopram treatment response
SNP, Single nucleotide polymorphism; CHR, chromosome; OR, odds ratio; NEDD4L, neural precursor cell expressed, developmentally down-regulated 4-like E3; ASTN1, astrotactin 1; LOC100505738, uncharacterized LOC100505738; SEMA5A, TNFRSF10B, tumour necrosis factor receptor superfamily, member 10b; SLC37A1, solute carrier family 37 (glycerol-3-phosphate transporter), member 1; PDE9A, phosphodiesterase 9A; FKBP1A, FK506 binding protein 1A; NSFL1C, NSFL1 (p97) cofactor (p47); LRCH1, leucine-rich repeats and calponin homology domain containing 1; NMD, nonsense-mediated decay.
a Minor allele.
b Odds ratio for minor allele.
c NMD transcript variant.
d Base pairs.
A logistic regression model on treatment response, including all covariates
We performed another logistic regression model on treatment response, entering the QoL factor simply as a covariate (no interaction factor), together with the rest of the covariates (age, sex, symptom severity, population clusters). Results were very similar to the interaction model: the top SNP was also rs520210 (p = 7.01 × 10−7). The top eight SNPs of this analysis were also found within the top SNPs of the interaction analysis. Top SNPs for this analysis are shown in online Supplementary Table S4, while a Manhattan plot is shown in online Supplementary Fig. S4. Results for the two models are very similar; however, higher p values result from the interaction model. From the latter regression model, we also examined each main effect of the covariates separately, including QoL. No SNP reached genomewide significance; however, QoL and population-cluster variables showed the most significant moderations; Manhattan plots for each covariate are presented in online Supplementary Fig. S5.
Enrichment analysis
10 × 105 Fisher's exact tests served for extrapolating the permutated p value describing the enrichment in the original lists of genes compared with 10 × 105 random ones, with respect to the association of single SNPs with the phenotype under analysis. Each of the 10 × 105 random groups was created as to meet two criteria: (1) being as large as the original one (same number of SNPs); and (2) enlisting SNPs that only belong to genes selected randomly in the genome. We found that the number of SNPs significantly (p < 0.05) associated with the phenotype in the original list was four times as large as what was expected by chance (21% compared with 5%, permutated (105) p value: 0.01 976). Further exploration of the original gene list showed that significant SNPs were mainly found within the serotonergic genes (Table 3). The MAOA, GABRA2, CREB1, CRHR1, NR3C1 and BDNF genes did not contain SNPs that were significantly associated with response. A full list of SNPs of the investigated genes and their level of significance is reported in online Supplementary Table S5.
Significant SNPs within significant genes of the pathway a
SNP, Single nucleotide polymorphism; CHR, chromosome; MAF, minor allele frequency; OR, odds ratio; CI, confidence interval; HTR2A, serotonin receptor 2A; SLC6A3, dopamine transporter; SLC6A4, serotonin transporter; TPH1, tryptophan hydroxylase 1; TPH2, tryptophan hydroxylase 2; COMT, catechol-O-methyltransferase; FKBP5, FK506 binding protein 5; DRD2, dopamine receptor D2.
a Total number of SNPs examined in: HTR2A, 39; SLC6A3, 18; SLC6A4, 9; TPH1, 9; TPH2, 32; COMT, 14; FKBP5, 9; DRD2, 19. For a full list of SNPs and genes within the pathway please refer to online Supplementary Table S5.
b OR for minor allele.
The association between the genes in the form of a biological network is presented in Fig. 1. According to functional enrichment analysis in the GeneMANIA program, most network genes are involved in the regulation of neurotransmission and drug binding, which appear as top functions of the network. Furthermore, the SLC6A4, SLC6A3 and DRD2 genes are involved in monoamine transport; the HTR2A and DRD2 genes are involved in G-protein coupled amine receptor activity, and the SLC6A3 and COMT genes are involved in the neurotransmitter biosynthetic process. A more detailed account of the functionality of the network of genes can be found in online Supplementary Table S6.
Network of the enriched genes. Network analyses showed that the genes involved in monoamine transport and regulation of neurotransmitter levels (DRD2, SLC6A3, SLC6A4; dopamine receptor D2, dopamine transporter, serotonin transporter, respectively) are interrelated; HTR2A (serotonin receptor 2A) is primarily involved in drug binding modulation and is interconnected with other serotonergic genes. MPD2, monophosphate deaminase 2; GRM2, glutamate receptor metabotropic 2; PICK1, protein interacting with PRKCA 1; AC133561.1, uncharacterized protein; YWHAE, tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, epsilon polypeptide; TPH2, tryptophan hydroxylase 2; TPH1, tryptophan hydroxylase 1; PAH, phenylalanine hydroxylase; TH, tyrosine hydroxylase; ACT-bd, amino acid binding ACT domain protein; FKBP5, FK506 binding protein 5; CHUK, conserved helix-loop-helix ubiquitous kinase; LRTOMT, leucine rich transmembrane and 0-methyltransferase domain containing; COMT, catechol-O-methyltransferase; LITAF, lipopolysaccharide-induced tumour necrosis factor (TNF) factor; COMTD1, catechol-O-methyltransferase domain containing 1. Further functions are shown in online Supplementary Table S6. The network was designed using GeneMANIA (University of Toronto; http://www.genemania.org/).
Discussion
The present study undertook two approaches in order to predict antidepressant treatment response. First, we examined whether genetic variations interact with the patients' levels of QoL to predict antidepressant response, after controlling for demographic characteristics, depression severity and population stratification. Second, we conducted an enrichment analysis, which explored whether candidate genes (that have emerged from prior gene–environment interaction studies) are associated with treatment response, in a genomewide context.
Genomewide interaction
Our model for predicting response to citalopram in the STAR*D sample yielded interesting candidate genetic loci. The top finding that survived genomewide significance involves a SNP (rs520210) of the NEDD4L gene, which has not been previously reported. We found that responders were most likely to be carriers of the minor allele (A), compared with non-responders, especially those reporting low QoL upon entering the trial. QoL probably explains a certain degree of variance in treatment response, since more pronounced genetic effects were observed in the interaction model and since QoL was shown to be among the most (statistically) relevant covariates (online Supplementary Fig. S5). However, a lower significance level (not reaching genomewide significance) was observed for the rs520210 SNP in the analysis with the Caucasian sample only, where it was among the top SNPs, when Hispanics were included, and it did not result among the top SNPs when Hispanics were excluded (n = 1017). Furthermore, few SNPs and genes were consistent across all models. Decreased statistical power due to the lower sample size could explain these discrepancies; however, other factors may also play a role. African-Americans were over-represented in the non-responder group, which may have influenced the results observed in the whole sample. Although statistically controlled for, effects of population stratification cannot be excluded. For this reason, we have also reported results in Caucasians separately (see online Supplementary Tables S1 and S2). Several SNPs showed a trend of association in the whole sample, but not in the smaller ethnic homogeneous samples, whereas a number of genes or intergenic regions (NEDD4L, LRCH1, SLC37A1, TNFRSF10B, SEMA5A) appeared among the top hits across models. The relative lack of consistency raises the possibility that some observations are false positive, largely dependent on sample size, or ethnic group, and easily disappear when group sizes change.
Nevertheless, our analyses yielded some interesting loci that can be further investigated in future studies. The NEDD4L (18q21.31) gene is the gene for encoding the ubiquitin protein ligase enzyme. Ubiquitins are small regulatory proteins that direct protein recycling and are then functional to the housekeeping activities of cells. Our result is consistent with the one reported by Garriock et al. (Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010) who found another ubiquitin E3 protein ligase (UBE3C) gene involved in antidepressant response in the same sample. Both genes are involved in the ubiquitin-dependent protein catabolic process, and inhibitor-SMAD (I-SMAD) binding. Ubiquitin proteasome pathway components (PSMD9, USP36) were recently observed in a pathway analysis exploring the risk for MDD (Wong et al. Reference Wong, Dong, Andreev, Arcos-Burgos and Licinio2012). Furthermore, the area where NEDD4L is located (namely, area 18q21–23) has been previously shown to contain significant loci in a GWAS of bipolar II disorder (Nwulia et al. Reference Nwulia, Miao, Zandi, Mackinnon, DePaulo and McInnis2007) and in a linkage study of bipolar I and bipolar II (Baron, Reference Baron2002). In addition, considering the role of the NEDD4L gene in processes involving sodium imbalances, polymorphisms in NEDD4L have been implicated extensively in hypertension (Kristjansson et al. Reference Kristjansson, Manolescu, Kristinsson, Hardarson, Knudsen, Ingason, Thorleifsson, Frigge, Kong, Gulcher and Stefansson2002; Araki et al. Reference Araki, Umemura, Miyagi, Yabana, Miki, Tamura, Uchino, Aoki, Goshima, Umemura and Ishigami2008), as well as in treatment efficacy of β-blockers in hypertensive patients (Svensson-Färbom et al. Reference Svensson-Färbom, Wahlstrand, Almgren, Dahlberg, Fava, Kjeldsen, Hedner and Melander2011). Association between this gene and hypertension is particularly interesting, since it has been recently shown that depressed patients under serotonergic antidepressant treatment have an increased risk of hypertension (Licht et al. Reference Licht, de Geus, Seldenrijk, van Hout, Zitman, van Dyck and Penninx2009). Future studies may want to further investigate this gene in G × E models of treatment response, with hypertension or related illness as the ‘E’ parameter.
From the top findings of the GWAS in the whole sample, several SNPs belong to the astrotactin (ASTN1) gene. The neuronal protein astrotactin (ASTN1) is a well-known receptor for glial-guided neuronal migration (Fishell & Hatten, Reference Fishell and Hatten1991; Zheng et al. Reference Zheng, Heintz and Hatten1996), a fundamental process in the development of the laminar structure of cortical regions of the brain (Wilson et al. Reference Wilson, Fryer, Fang and Hatten2010). Interestingly, a pharmacogenetic GWAS on antidepressant response in an independent sample (GENome-based therapeutic drugs for DEPression; GENDEP) located a gene [UST (uronyl-2-sulfotransferase) gene, 6q25.1)], which also plays a role in neuronal migration (Uher et al. Reference Uher, Perroud, Ng, Hauser, Henigsberg, Maier, Mors, Placentino, Rietschel, Souery, Zagar, Czerski, Jerman, Larsen, Schulze, Zobel, Cohen-Woods, Pirlo, Butler, Muglia, Barnes, Lathrop, Farmer, Breen, Aitchison, Craig, Lewis and McGuffin2010). The TNFRSF10B gene is another interesting gene that resulted among the top 10 SNPs in the whole sample and in Caucasians (including Hispanics). This gene binds to the TRAIL protein [tumour necrosis factor (TNF)-related apoptosis-inducing ligand, a ligand that induces the process of cell death], and plays a major role in the apoptotic process and activation of nuclear factor-κB (NF-κB)-inducing kinase activity. NF-κB has been reported as a mediator of stress-impaired neurogenesis and depressive-like phenotypes (Koo et al. Reference Koo, Russo, Ferguson, Nestler and Duman2010), and mediator of stress response (LaPlant et al. Reference LaPlant, Chakravarty, Vialou, Mukherjee, Koo, Kalahasti, Bradbury, Taylor, Maze, Kumar, Graham, Birnbaum, Krishnan, Truong, Neve, Nestler and Russo2009), in animal models.
Furthermore, other genes that were consistently observed in the Caucasian, as well as in the whole, sample were an intergenic region near SEMA5A, which has been previously associated with autism (Weiss et al. Reference Weiss, Arking, Daly and Chakravarti2009), intergenic regions near SLC37A1 (which is a glycerol-3-phosphate transporter), and the LRCH1 gene; the latter two genes have not been associated yet with psychiatric disorders. Finally, intergenic variations near FKBP12 arose as promising loci. FKBP12 is the prototype of the FKBP family and binds immunosuppressive drugs, such as FK506 and rapamycin (Harding et al. Reference Harding, Galat, Uehling and Schreiber1989). From the FKBP group, FKBP5 has been previously associated with the recurrence of depressive episodes and antidepressant response in an independent study (Binder et al. Reference Binder, Salyakina, Lichtner, Wochnik, Ising, Putz, Papiol, Seaman, Lucae, Kohli, Nickel, Kunzel, Fuchs, Majer, Pfennig, Kern, Brunner, Modell, Baghai, Deiml, Zill, Bondy, Rupprecht, Messer, Kohnlein, Dabitz, Bruckl, Muller, Pfister, Lieb, Mueller, Lohmussaar, Strom, Bettecken, Meitinger, Uhr, Rein, Holsboer and Muller-Myhsok2004), and in the STAR*D trial using a candidate gene approach (Lekman et al. Reference Lekman, Laje, Charney, Rush, Wilson, Sorant, Lipsky, Wisniewski, Manji, McMahon and Paddock2008). Furthermore, the FKBP5 modulated effects of childhood trauma in predicting major depression (Appel et al. Reference Appel, Schwahn, Mahler, Schulz, Spitzer, Fenske, Stender, Barnow, John, Teumer, Biffar, Nauck, Volzke, Freyberger and Grabe2011; Zimmermann et al. Reference Zimmermann, Bruckl, Nocon, Pfister, Binder, Uhr, Lieb, Moffitt, Caspi, Holsboer and Ising2011). The FKBP5 gene was also one of the significant genes implicated in our enrichment analysis.
Enrichment
We found that serotonergic genes implicated in previously published G × E interaction studies are part of a network, which is likely to affect antidepressant treatment response. Variation in the HTR2A gene has been associated with antidepressant response or remission in a number of studies (Uher et al. Reference Uher, Huezo-Diaz, Perroud, Smith, Rietschel, Mors, Hauser, Maier, Kozel, Henigsberg, Barreto, Placentino, Dernovsek, Schulze, Kalember, Zobel, Czerski, Larsen, Souery, Giovannini, Gray, Lewis, Farmer, Aitchison, McGuffin and Craig2009; Horstmann et al. Reference Horstmann, Lucae, Menke, Hennings, Ising, Roeske, Muller-Myhsok, Holsboer and Binder2010; Kishi et al. Reference Kishi, Yoshimura, Kitajima, Okochi, Okumura, Tsunoka, Yamanouchi, Kinoshita, Kawashima, Naitoh, Nakamura, Ozaki and Iwata2010; Lucae et al. Reference Lucae, Ising, Horstmann, Baune, Arolt, Muller-Myhsok, Holsboer and Domschke2010). The HTR2A gene was previously studied in the STAR*D sample, using a candidate gene approach, by McMahon et al. (Reference McMahon, Buervenich, Charney, Lipsky, Rush, Wilson, Sorant, Papanicolaou, Laje, Fava, Trivedi, Wisniewski and Manji2006), who examined 68 candidate genes and found significant associations with only two intronic variants in the HTR2A and treatment response. HTR2A influence on citalopram response was reproduced by Paddock et al. (Reference Paddock, Laje, Charney, Rush, Wilson, Sorant, Lipsky, Wisniewski, Manji and McMahon2007) who used a larger STAR*D sample and examined another large number of candidate markers. To date, several variants within HTR2A have been proposed to affect antidepressant response, but failures to replicate findings also exist (Drago et al. Reference Drago, De Ronchi and Serretti2009). By using a different methodological approach from the above studies, our observations from the enrichment analysis further support the role of the serotonergic genes in influencing antidepressant response to citalopram. A variant in the TPH1 gene also showed a strong effect (OR = 3.2), which could be worthy of further investigation. Furthermore, several significant SNPs were found within the FKBP5 gene, which has been previously associated with treatment outcome in depression in two independent studies (Binder et al. Reference Binder, Salyakina, Lichtner, Wochnik, Ising, Putz, Papiol, Seaman, Lucae, Kohli, Nickel, Kunzel, Fuchs, Majer, Pfennig, Kern, Brunner, Modell, Baghai, Deiml, Zill, Bondy, Rupprecht, Messer, Kohnlein, Dabitz, Bruckl, Muller, Pfister, Lieb, Mueller, Lohmussaar, Strom, Bettecken, Meitinger, Uhr, Rein, Holsboer and Muller-Myhsok2004; Kirchheiner et al. Reference Kirchheiner, Lorch, Lebedeva, Seeringer, Roots, Sasse and Brockmöller2008). In the STAR*D trial, it was previously examined using a candidate gene approach (only two SNPs) and an association was found with remission (Lekman et al. Reference Lekman, Laje, Charney, Rush, Wilson, Sorant, Lipsky, Wisniewski, Manji, McMahon and Paddock2008).
Challenges and limitations
Our sample size has insufficient power to detect interaction effects (Murcray et al. Reference Murcray, Lewinger, Conti, Thomas and Gauderman2011). Nevertheless, we hypothesize that by considering a multifactorial, or more individualized, model of antidepressant response, genetic effects may become more detectable. There is increasing evidence emerging from GWASs that there are no genetic variants that contribute with a large effect (Flint & Munafo, Reference Flint and Munafo2013). Hence, apart from the large sample sizes that would provide sufficient power to observe small effects, designing multivariate models of prediction with the tools available in each study could assist in unravelling variants of risk (Cohen-Woods et al. Reference Cohen-Woods, Craig and McGuffin2013).
Several limitations need to be acknowledged. A major limitation with regard to the QoL measure is that the level of depression severity of the patient could have biased responses and might be intertwined with the measure itself. Exploring the relationship between depression severity and QoL in our sample, we found a moderate correlation (r = 0.54); hence, approximately 29% (r 2 = 0.29) of a patient's QoL can be explained by his/her levels of depressive symptoms upon entering the trial. Consequently, we statistically controlled for the severity of depression in our analyses. Depression severity is an inherent part of the disorder; in fact, we can diagnose the disorder due to the level of severity. QoL, however, also contains aspects on physical health and quality of social relationships that could be independent of (or not affected by) the individual's symptomatology. This is supported by a factor analysis on the QoL measure showing distinct factors ‘physical health’, ‘leisure time activities’, ‘social relationships’ and ‘subjective feelings’, the latter being more likely to overlap with variance explained by depressive symptoms (Ritsner et al. Reference Ritsner, Kurs, Gibel, Ratner and Endicott2005). Interestingly, a recent study in the STAR*D sample showed that QoL was a factor that contributed significantly to the patients’ Individual Burden of Illness Index, and the variance explained by QoL was additional to that explained by depressive symptoms in a principal component analysis (Cohen et al. Reference Cohen, Greenberg and Ishak2013). These additional aspects (satisfaction with work, health issues, social relationships, family relationships, economic condition, and leisure time) represented by the QoL measure could be considered as ‘external’ conditions, possibly approaching a proxy of an environmental contributor. Apart from depressive symptoms, several other factors that may have affected reports on QoL could be the duration of the episode, education, family status and other variables. Within the STAR*D trial, factors related to QoL have been previously examined and only age and depression severity were shown to have an important contribution (Daly et al. Reference Daly, Trivedi, Wisniewski, Nierenberg, Gaynes, Warden, Morris, Luther, Farabaugh, Cook and Rush2010), both factors being taken into account in the present analysis. Nevertheless, future studies would benefit from more robust and independent assessments. Ultimately, the assessment of important environmental adversities, such as childhood maltreatment, would be an important step in genomewide environment interactions (Karg & Sen, Reference Karg and Sen2012), since such events play a significant role in the onset of depression and response to treatment (Uher, Reference Uher2011). It would be interesting, therefore, to explore these G × E interactions in large samples, such as the GENDEP or MARS (Munich Antidepressant Response Signature), where environmental variables may be available. Finally, other factors that are likely to moderate the relationship between genes and antidepressant response could not be reported. For example, epigenetic factors possibly conceal some important effects in such investigations and are worthy of future research. Our study also suffers from other limitations of the STAR*D trial, as already reported previously (Laje et al. Reference Laje, Perlis, Rush and McMahon2009; Garriock et al. Reference Garriock, Kraft, Shyn, Peters, Yokoyama, Jenkins, Reinalda, Slager, McGrath and Hamilton2010).
Our study provides new insights for interesting genetic markers and mechanistic pathways; however, the variance explained is still modest. In addition, our enrichment analysis revealed that certain genes implicated in prior G × E studies are likely to affect antidepressant treatment response. These findings provide further incentives for research into the mechanisms of the proposed gene-related proteins in relation to antidepressant action.
Supplementary material
For supplementary material accompanying this paper visit http://dx.doi.org/10.1017/S0033291713001554.
Acknowledgements
The authors thank the STAR*D study team for conducting the study, including the University of Texas Southwestern Medical Center (National Coordinating Center) Dallas; the Massachusetts General Hospital, Boston; the University of Pittsburgh Medical Center, Pittsburgh; the Columbia College of Physicians and Surgeons, New York; the Laureate Healthcare System, Tulsa; the University of North Carolina, Chapel Hill; the University of California at Los Angeles; the University of California at San Diego; the Northwestern University Medical School, Chicago; the Psychiatric Research Institute, Wichita, Kansas; the University of Michigan, Ann Arbor; the Vanderbilt University Medical Center, Nashville; the Tuscaloosa Veterans Affairs Medical Center; and the Virginia Commonwealth University, Richmond (for further details, please visit http://www.edc.pitt.edu/stard/). We thank the NIMH for placing the data in the public domain, and we also thank the patients and their families for having agreed to be part of the STAR*D project.
The work of N.A. is supported by the ‘Rubicon’ Grant (no. 446–11–004), awarded by the Netherlands Organization for Scientific Research and the Marie Curie Cofund Action. The funding bodies had no involvement in the research.
Declaration of Interest
A.S. is or has been a consultant/speaker for Abbott, Angelini, AstraZeneca, Clinical Data, Boheringer, Bristol–Myers Squibb, Eli Lilly, GlaxoSmithKline, Italfarmaco, Janssen, Lundbeck, Pfizer, Sanofi and Servier.
