|Home | About | Journals | Submit | Contact Us | Français|
Antidepressant response is likely influenced by genetic constitution, but the actual genes involved have yet to be determined. We have carried out a genome-wide association study to determine if common DNA variation influences antidepressant response.
Our sample is derived from Level 1 participants in the Sequenced Treatment Alternatives to Relieve Depression (STAR*D) study, all treated with citalopram. Association for the response phenotype included 883 responders and 608 non- responders. For the remission phenotype, 743 subjects that achieved remission were compared to 608 non-responders. We used a subset of SNPs (n = 430,198) from the Affymetrix 500K and 5.0 Human SNP Arrays, and association analysis was carried out after correcting for population stratification.
We identified three SNPs associated with response with p-values less than 1 × 10−5 near the UBE3C gene (rs6966038, p = 4.65 × 10−7), another 100kb away from BMP7 (rs6127921, p = 3.45 × 10−6), and a third that is intronic in the RORA gene (rs809736, p = 8.19 × 10−6). These same SNPs were also associated with remission. Thirty-nine additional SNPs are of interest with p-values ≤ 0.0001 for the response and remission phenotypes.
Although the findings reported here do not meet a genome-wide threshold for significance, the regions identified from this study provide targets for independent replication and novel pathways to investigate mechanisms of antidepressant response. This study was not placebo controlled, making it possible that we are also observing associations to non-specific aspects of drug treatment of depression.
Depression is a serious health concern with a lifetime prevalence of ~16% in the United States (1). Symptoms of Major Depressive Disorder (MDD) are often alleviated with use of antidepressant medication, but about 65% of those afflicted with MDD will not achieve remission after their first medication trial (2). Identification of patient factors that increase the likelihood of successful treatment would benefit patients, clinicians, and society. Such factors may include inter-individual genetic differences that contribute to variability in antidepressant response. It has been observed that antidepressant response clusters in families, supporting a role for genetic variation in antidepressant response, though the degree of heritability is still unknown (3, 4). Previously, we reported candidate gene association studies for pharmacokinetic genes (5) and serotonin-pathway genes (6, 7) in the Sequenced Treatment Alternatives to Relieve Depression (STAR*D) sample, which is comprised of 1,953 subjects treated with citalopram. Although we failed to find statistically significant associations between antidepressant response and these candidate genes, other interesting candidate genes have been reported for response (8, 9) and remission (10, 11) in the same sample. Given the lack of placebo in the STAR*D trial, it is possible that previous associations may be related to non-specific aspects of the treatment protocol.
An alternative to candidate gene studies are genome-wide association studies (GWAS), which have successfully identified genetic variants for numerous diseases, including type 2 diabetes (12) and Crohn’s Disease (13). Several pharmacogenetic GWASs have also identified novel candidate regions for phenotypes of interest including hydrochlorothiazide response for treatment of hypertension (14), anti-tumor necrosis factor response for rheumatoid arthritis (15), interferon beta therapy for multiple sclerosis (16) and iloperidone response for schizophrenia (17). Using the STAR*D sample, we carried out a GWAS to determine if common DNA variation influences citalopram response and remission.
The design of the Sequenced Treatments Alternatives to Relieve Depression (STAR*D) study has been previously described (2, 18, 19). To maximize power to detect associations, we focused on Level 1 data collected in the STAR*D trial, where all participants were initially treated with citalopram. Our DNA samples and phenotypes have been described in previous reports (5, 20). For this report, our outcome phenotypes were antidepressant response and remission, the definitions of which have been previously described (5). Briefly, the first comparison is between responders and non-responders. Both groups required a 42-day enrollment in the study with responders demonstrating a ≥ 50% reduction in the 16-item Quick Inventory of Depressive Symptomatology, self-report version (QIDS-SR) (21) from baseline to final visit. Non-responders did not meet this criterion. We use “response” to denote drug efficacy. For remission, we compared non-responders as described above to individuals scoring ≤ 5 on the QIDS-SR at follow-up, which closely corresponds with the conventional definition of remission by the Hamilton Rating Scale for Depression (HAM-D) ≤ 7. Figure S1 details the inclusion and exclusion of subjects for the association study. Tolerance was based on exit data for the 1st level of the STAR*D study. If a subject continued on citalopram or entered follow-up at the end of level 1, they were considered tolerant to citalopram. If the level-exit or study-exit forms indicated side-effects, they were considered intolerant. Figure S1 details the inclusion and exclusion of subjects for the association study. For secondary analysis, response was considered as a quantitative variable and we measured the percent change in QIDS-SR from baseline.
Genotyping of the 1,948 STAR*D samples was conducted at two locations and on two different platforms. A total of 969 subjects were genotyped at Affymetrix, Inc. (South San Francisco) on the Human Mapping 500K Array Set. We genotyped the remaining 979 samples using the Affymetrix Genome-Wide Human SNP Array 5.0. The two groups were balanced by ethnic grouping, gender and proportions of responders and non-responders. Samples run on Affymetrix 500K Array were called using the BRLMM algorithm, and samples analyzed on Affymetrix Array 5.0 were called using the BRLMM-P algorithm. Twelve samples were genotyped on both the 500K and 5.0 Arrays, and we found > 99% concordance across these platforms.
Hardy-Weinberg Equilibrium was evaluated using either Fisher’s exact test or chi-square test within each ethnic group, as appropriate. SNPs with P < 1.0 × 10−5 in all three ethnic categories (African-American, Non-Hispanic Caucasian, and Hispanic Caucasian) were excluded (n= 1,233 SNPs). This was done to remove the poorly performing SNPs that could be a result of genotyping error. The requirement for the SNP to violate HWE in all three ethnic groups derived from the assumption that genotyping error would not be specific to one ethnic group. SNPs with minor allele frequencies < 1.0% (n=10,792) or call rate < 95% (n=42,908) were excluded. The final number of SNPs for autosomal analysis was 430,198. Differential missingness between responders and non-responders was assessed, and considered to not influence the association statistics. Table S13 lists the 63 SNPs with differential missingness p-values < 1.0 × 10−3. Only one SNP (rs2443527, rank 8,431) fell within the top 10,000 association hits for remission, therefore, no SNPs were dropped due to differential missingness. SNPs on the X-chromosome (N=9,697) that passed the same quality control measures as the autosomal SNPs were also tested for association with response and remission. The estimated average rate of coverage of common variants across the genome is 46%. This is based on an r2 ≥ 0.8 between tagging SNPs and non-panel HapMap SNPs (release 21) in the 10 ENCODE regions (22). We estimated coverage using of these regions because they are neutral with respect to choice of location and density of the SNPs on the Affymetrix array we used for our association study.
Individuals were excluded if they were of self-reported Native American ancestry (n=16), Asian ancestry (n=21), Pacific Island ancestry (n=8), mixed ancestry (n=68), or unknown (n=1) due to small numbers in these groups. Consequently, for the response phenotype, we analyzed 1,067 Non-Hispanic Caucasians (71.5%), 241 African-Americans (16.2%) and 183 Hispanic Caucasians (12.3%). Further, we excluded individuals if they had more than 5% of their genotypes missing (n=71), did not have follow-up phenotypic data (n=16), or were cryptically related to others in the sample (n=32), which was determined using the 1-IBS metric tabulated in PLINK v1.04 (22). The final number of individuals in the analysis was 1,491 for response vs. non-responders, and 1,351 for remitters vs. non-responders.
Ancestry was controlled for using multidimensional scaling (MDS), which is a multivariate method that summarizes all SNPs into linear combinations of the count of rare alleles of SNPs. These MDS dimensions are linearly independent of each other, and are given consecutive numbering (i.e. MDS #1, MDS #2, etc.) based on how much variation in SNP rare allele count they explain. We used 140,701 independent SNPs in the MDS analysis. Independence was determined using the –indep function in PLINK for the following conditions: 50 SNPs in the window, five SNP window shift, and a variance inflation factor of two. The first ten MDS dimensions were then used as covariates in association analyses to adjust for population stratification, as previously established (23). Q-Q plots of expected and observed p-values are available for response and remission in Figures S2 and S3, and were created using publicly available scripts (24, 25). The inflation factor λ was computed based on the Cochran-Armitage trend tests.
All SNPs were tested for association with response and remission using logistic regression. The minor allele homozygous genotype was used as the reference group, and each SNP was modeled in individuals as having a log-additive effect, after adjusting for the ten MDS dimensions. Odds ratios (ORs) and 95% confidence intervals (CIs) are reported. Statistical analysis was done using PLINK v1.04 (26). Cluster plots for all SNPs contributing to results reported in Tables 1 and and22 were examined for artifacts (see Figure S5).
We estimated statistical power using Quanto (27). We used the following as parameters: α =1.16 × 10−7 (Bonferroni correction for 430,198 markers), risk of response (0.59) and remission (0.55), 1−β = 0.8, log additive model, and minor allele frequency of 0.25. For the response phenotype (883 responders and 608 non-responders), the genotypic relative risk that could be detected at the set power level was 1.74. For the remission phenotype (743 remitters and 608 non-responders), the equivalent detectable genotypic relative risk was 1.77. A range of power estimations for both phenotypes is graphically depicted in Figure S6.
Genotype and phenotype data are available from the NIMH Center for Collaborative Genetic Studies on Mental Disorders.
We evaluated 430,198 SNPs for association with antidepressant response and remission. For response, 1,491 (608 non-responders; 883 responders) subjects were analyzed, and for remission, 1,351 (608 non-responders; 743 responders) subjects were analyzed. The gender distribution between the groups within each analysis was not significantly different (Table S1). There were significant differences by ethnicity between responders and non-responders (p = 0.003) and between remitters and non-responders (p = 0.0003). African Americans were over-represented among non-responders (19.7%), compared to responders (13.7%) and remitters (12.5%). Significant differences among responders and non-responders also existed for a number of demographic and clinical measures (e.g., insurance and employment status, clinical co-morbidities, etc), as shown in Table S1. Since these variables are correlated with the response phenotype, they are considered potential confounding variables. Therefore, we tested association of the top 100 SNPs associated with response for association with insurance status, employment status, years of completed schooling, history of attempted suicide, and co-morbid disorders (generalized anxiety, hypochondriasis, panic disorder, somatoform disorder, social phobia, and obsessive compulsive disorder). None of these variables are associated with the top reported SNPs, and therefore are not considered confounding variables for the response and remission association results, as confounding variables must be associated with outcome (e.g. response) and the genetic variant (e.g. our list of top 100 SNPs).
We used the same approach for assessment of tolerance as a confounding variable. The top 100 SNPs from the response and remission association tests were tested for association with tolerance. This resulted in 145 SNPs be tested for association with tolerance, as there is overlap in the association results between response and remission. Seventeen SNPs were found to have association p-values < 0.05 with tolerance, with the smallest p-value = 0.0003, and are presented in Table S2. Only one SNP (rs11128623) met a Bonferroni corrected p-value (p= 0.00034) for significance given that 145 SNPs were evaluated, and thus could be considered as confounded by tolerance.
Q-Q plots and full genome plots for the association tests of both the response and remission phenotypes are found in supplemental data (Figures S2, S3 and S7). The genomic inflation factors for both phenotype association tests are 1.01, suggesting that our association results are unlikely to be confounded by population stratification. Figure S7 shows regional plots for the two strongest association results for response and remission. We found 39 SNPs and 41 SNPs with p-values < 1.0 × 10−4 for the response and remission phenotypes, respectively. Table 1 lists the top 25 findings for the response phenotype and Table 2 lists the top findings for the remission phenotype. A more comprehensive list of the top 100 results for both phenotypes is presented in Tables S3 and S4. Thirteen of the top 25 SNPs from results for the response phenotype analysis are also found in the top 25 findings for the remission phenotype. This is not surprising given the correlation between the two phenotypes. The top two results for response and remission are the same SNPs: rs6966038 and rs6127921, with p-values for response of 4.65 × 10−7 and 3.45 × 10−6, and for remission, 3.63 × 10−7 and 1.07 × 10−6. We sought to determine if additional SNPs in these regions of interest also showed evidence for association (Tables S5 and S6). Our top-ranked finding for both phenotypes does not have any additional SNPs with a p-value < 0.05 and an r2 ≥ 0.50 within 250kb. However, rs6127921 (rank 2) has a nearby SNP with a p-value < 0.001. We analyzed the surrounding regions for the most associated SNPs for proxies in the region, using the --proxy function in PLINK. For the top SNP finding (for both response and remission), we found no proxies for this. However, the ombnibus haplotype test supported association in the region (p = 1.19 × 10−5 for response, and 3.92 × 10−7 for remission). For the second ranked SNP associated with response, we found an omnibus haplotype test p = 0.00065, with 2 of 6 subhaplotypes meeting proxy criteria with p-values 0.00036 and 0.00046. These analyses suggest although no single SNPs show similar associations to the most associated SNPs in these regions, we found haplotypes that do show association results of a similar magnitude and that are correlated with the original associated SNPs. The top 10 association results for response and remission on the X chromosome and pseudo-autosomal regions are in Table S7. No statistically significant association was detected with these phenotypes on the sex chromosomes.
The ancestry groups were also analyzed separately for the response phenotype. The top 25 hits for each ancestry group (Non-Hispanic Caucasians, Hispanic Caucasians, and Non-Hispanic African Americans) are presented in Table S8. There is no correlation of the SNP rankings between the ancestry groups, most likely due to low power to detect risk variants in the two smaller sub-samples, or differential risk loci related to ethnicity. In an additional exploratory analysis where response was considered as a quantitative trait (percent change in QIDS-SR from baseline), the lowest p-value was 3.20 × 10−6 for rs6508329. The top 100 association p-values for this analysis are presented in Table S9.
We found high correlation between SNP ranking in association results when analyzed with 1 multi-dimensional scaling (MDS) covariate vs. 10 MDS covariates (r2 = 0.974). Table S10 demonstrates that whether we correct for 1 or 10 covariates, the results are largely the same. Further, as an alternative to SNP ranking we also evaluated change in odds ratios (ORs) and standard deviations and found little change across these comparisons (e.g., mean difference in ORs = 0.00013), suggesting they are comparable metrics.
Many candidate gene studies of antidepressant response have preceded this current GWAS, some with positive and many with null associations. We assessed our data for how these past candidates performed using WGAviewer (28). A list of the investigated genes, and their corresponding ranks and p-values from the GWAS data are listed for the response and remission phenotypes (Table S11). Additionally, the gene coverage [coverage rate = (tagged SNPs in HapMap + genotyped SNPs in GWAS set)/total common SNPs in HapMap; an algorithm created in WGA viewer (28, 29)], based on our SNP set, and numbers of SNPs with p-value < 0.05 are listed in Table S12. The serotonin receptor 2A (HTR2A) has several nominally associated SNPs for response and remission. This supports previous findings of HTR2A being involved in antidepressant response using the same STAR*D sample used here (7, 9). The D2 dopamine receptor (DRD2) shows some evidence of association for one SNP (rs2234689, p-value = 0.0001 for remission). Notably, no support for the serotonin transporter (SLC6A4) was detected in remission or response analyses. This finding also supports our previous report of no association in this gene in the same sample (6).
Our findings represent a GWAS conducted to determine if genetic variation influences population variation in antidepressant response and remission in Major Depressive Disorder in the STAR*D sample. Our single locus findings for response and remission include two SNPs located between two sets of genes. The top finding (rs6966038; response p = 1.65 × 10−7) is 51kb from the ubiquitin protein ligase E3C gene (UBE3C) and 77kb from the motor neuron and pancreas homeobox 1 (MNX1 or HLXB9) gene. The UBE3C protein modifies proteins to signal them for degradation (30) and the MNX1 protein is a DNA binding protein that controls gene expression at the caudal end of the developing notochord of an embryo (31). The second finding (rs6127921; response p = 3.45 × 10−6) is closest to the bone morphogenic protein 7 (BMP7) gene, but at least 106kb away. Many of the SNPs that appear in the top tier of results in the response phenotype also appear in the remission phenotype. This is likely due to the overlap in samples and phenotypic definition. In our case, all remitters are also responders, however, not all responders achieve remission, and thus remitters are a subgroup of the responders. Furthermore, both phenotypes were compared to a common non-responder group.
When we examined the local regions around our highest ranked findings, we found four gene regions supporting association to response and remission: amiloride-sensitive cation channel 1, neuronal (ACCN1), PR domain-containing protein 2 (PRDM2), aryl hydrocarbon receptor translocator-like (ARNTL) and nucleolar protein 4 (NOL4) genes. More detailed information regarding these genes is available in the Supplemental Materials.
Although these data potentially implicate novel pathways for investigation of mechanisms of antidepressant response, we do not believe that these data can reliably be incorporated into genetic tests for prediction of antidepressant response. Rigorous replication in a larger, independent sample demonstrating significant signals of similar effect size and direction would be required before a solid correlation between these variants and antidepressant response could be surmised. The effect sizes for our highest ranked SNPs are quite low and, though interesting, do not justify a clinical application at this early stage of investigation. While this study represents the largest pharmacogenetic study of antidepressant response, our sample was only large enough to detect genotypic relative risks in the range of 1.75 or greater, and we observed several findings with odds ratios in that range. Given the limited power in our sample, meta-analysis of all available samples with similar pharmacogenetic phenotypes will be required before any definitive conclusions can be generated. Similarly, while our data could be construed as not supporting the role of common genetic variation in citalopram response, this conclusion may be premature from a single study. It is possible that larger samples may detect true risk alleles of small effect, but still of scientific interest. It is also possible that higher density genotyping of common alleles not covered by our genotyping platform would uncover missed associations.
Imputation of these data was not performed due to the admixed and multi-ethnic nature of our sample. The choice of a reference panel on which to base the imputation is not straightforward and complicates an imputation procedure. Imputation would be useful to show regional support for a SNP that exceeds typical genome-wide statistical significance values; however we are not reporting that any of our SNPs surpass this threshold. Likewise, this report does not include copy number variation (CNV) analyses due to the technical challenges of calling CNVs in Affymetrix data when using two separate genotyping platforms and genotyped at two locations. Furthermore, the power to detect an association of CNVs with antidepressant response would be marginal given their rarity.
The limitations to our study design have been previously reviewed (5, 6, 56). Briefly, they include: phenotypic definitions, STAR*D study design for placebo and adherence, lack of generalizability to the general affected population and to other antidepressant medications, and population stratification. Each is discussed in turn.
Our phenotypes are antidepressant response and remission, and require a six week exposure to citalopram to be considered for analysis. This should theoretically increase our chances of finding a true drug response if we are excluding subjects who responded to treatment before the six week timeline; though it is possible that we are missing non-responders that dropped out of the study early. In a post-hoc analysis, we removed the six week treatment criterion, and found that inclusion of additional individuals did not improve the strength of our findings (Table S14 and Figure S8). Another concern is that the STAR*D trial did not have a placebo arm in the clinical trial design, and did not measure serum levels or other levels of adherence to treatment. Under these conditions we cannot thoroughly assess true response to treatment versus placebo or be sure that all responders are responding due to the medication. This is a situation shared with all treatment studies in depression which attempt to discover predictors of response are faced with, namely there is still no agreed upon method to separate true drug response from placebo response. Thus, in the absence of placebo or comparator in Level 1, it is possible that we are not studying citalopram response, per se, but that our observations may actually relate to the speed of recovery or the placebo effect. This intriguing possibility provides grounds for further research. Likewise, and more relevant is that we cannot know for sure that non-responders were actually taking their medication, and might be lacking a response due to non-compliance. Furthermore, we lack information on concomitant medications that could interact with citalopram and alter its effects (positively or negatively) and lead to an artificial inflation of the phenotype (response or non-response). Although this is a real issue in this design, STAR*D response rates are similar to clinical trials in which adherence was measured, suggesting that the rates observed in the STAR*D sample may be fairly accurate.
Another limitation of the STAR*D study design for genome wide association studies is the initiation of DNA collection after the start of the study. As a result, not all individuals in the STAR*D study donated DNA, and there are some differences between individuals who donated DNA and those who did not, potentially reducing the generalizability of our results to the entire STAR*D sample. These differences are outlined in previous reports (6) and include variables such as education, marital status, age, ethnicity and severity of disease.
The STAR*D study design was intended to mimic naturalistic clinical practices as closely as possible and allowed dosage adjustments as necessary, presenting a possible confound for a GWAS on antidepressant response phenotypes. We found this not to be the case, since there was no significant difference between responders and non-responders for final citalopram dose. Our findings may be specific to citalopram treatment, and may not be the same for treatment with other selective serotonin reuptake inhibitors (SSRI’s) or other classes of antidepressant medication. Notably, this is the largest collection of samples designed to study antidepressant response to date. Until another study of similar magnitude is conducted with a different SSRI or another class of treatment, we are unable to determine if the findings reported here are relevant for response to medications other than citalopram.
It is possible that for very few of the SNPs, the response or remission association results could be confounded by tolerance. However, only one SNP showed possible evidence of confounding, suggesting that confounding from tolerance is not influencing our results.
Population stratification is a concern in genetic association studies, regardless of phenotype. We controlled for population structure by using MDS vectors as covariates in our logistic regression. This approach has been shown to be effective in accounting for confounding due to population structure (23), and is additionally supported by our λ estimates of 1.01 (for both response and remission), which also strongly suggests limited confounding due to hidden population structure in our data.
In summary, we report results of the first genome wide association study on antidepressant response phenotypes in the STAR*D sample. These data will prove useful for meta-analysis of similar antidepressant GWAS datasets. Additionally, they provide a basis for several regions or specific genes to be pursued in fine-mapping and independent replication studies. Furthermore, these data invite investigation of the associated proteins in novel mechanisms of antidepressant response.
This work was supported by the National Institute of Mental Health (NIMH) to S.P.H. (MH072802); by NIMH training funds to H.A.G. (F32 MH082562-01A1 & T32 MH19552-11), and to S.I.S. (R25 MH060482-08 & T32 MH19126-19); by NARSAD New Investigator Award A109584 (H.A.G.); by a Howard Hughes Medical Institute Pre-Doctoral Fellowship to E.J.P.; and by the State of New York for partial support to P.J.M. for this work. The authors appreciate the efforts of the STAR*D Investigator Team for acquiring, compiling, and sharing the STAR*D clinical dataset. STAR*D was funded by NIMH via contract (N01MH90003) to the University of Texas Southwestern Medical Center at Dallas (A. John Rush, P.I.). We would also like to thank Stephen Wisniewski, Ph.D., Director, STAR*D Data Coordinating Center, University of Pittsburgh for providing demographic and clinical response data.
Dr. McGrath reports receiving research support from the National Institute of Mental Health, the New York State Department of Mental Hygiene, and Roche Pharmaceuticals. The other authors reported no biomedical financial interests or potential conflicts of interest.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.