Genome-Wide Pharmacogenetics of Antidepressant Response in the GENDEP Project

After 50 years of experience with antidepressant drugs, the large individual variation in antidepressant treatment outcome remains poorly understood. Investigations of biologically related individuals point to a genetic contribution to treatment responsiveness (1), and several polymorphisms in candidate genes have been associated with antidepressant response (2, 3). However, few pharmacogenetic associations have been replicated, and they explain only a small fraction of individual differences in antidepressant response. Since our understanding of the mechanisms of action of antidepressants is incomplete, important genes and regulatory intergenic sequences may have been missed in candidate gene studies. Therefore, a systematic exploration of variation across the genome has the potential to detect further variants that may help in understanding the biology of antidepressant action and predict response to a specific antidepressant in an individual patient. In the present study, we report on a genome-wide pharmacogenetic analysis of more than 500,000 common genetic variants, tested for association with change in depression severity over a 12-week period after commencing treatment with a serotonergic or noradrenergic antidepressant among subjects with moderate to severe depression. We conducted four primary analyses to address general and drug-specific effects of genes on antidepressant response and gene-drug interactions. We used secondary analyses to explore the influence of variants in 72 candidate genes on antidepressant action.

Method

Study Design and Sample

The Genome-Based Therapeutic Drugs for Depression (GENDEP) project (3, 4) is a 12-week partially randomized open-label pharmacogenetic study with two active treatment arms. A total of 811 adult patients (men: N=297; women: N=514) with unipolar depression of at least moderate severity according to ICD-10 or DSM-IV criteria were recruited at nine European centers. Patients were aged 19–72 years and of Caucasian European parentage. Diagnoses were established using the semistructured Schedules for Clinical Assessment in Neuropsychiatry interview (5). Exclusion criteria were personal and family history of schizophrenia or bipolar disorder and current substance dependence. Eligible participants were allocated to treatment with either escitalopram or nortriptyline, which represent different mechanisms of antidepressant action. Escitalopram is a selective inhibitor of the serotonin transporter, with no effect on norepinephrine reuptake (6). Nortriptyline is a tricyclic antidepressant, with an affinity for the norepinephrine transporter that is 100 times higher than for the serotonin transporter (7). Patients with no contraindications were randomly allocated to flexible-dosage nortriptyline (50–150 mg daily) or escitalopram (10–30 mg daily) for 12 weeks. Patients with contraindications for one drug were offered the other. Severity of depression was assessed weekly by three established rating scales (8). A total of 628 participants (77%) completed at least 8 weeks of treatment with the originally allocated antidepressant (4). The GENDEP project was approved by ethics boards of participating centers, and all participants provided written informed consent. Previous investigations of the GENDEP project sample included three candidate gene pharmacogenetic studies (3, 9, 10).

Definition of Phenotype

Of the three scales used in the GENDEP project, the clinician-rated 10-item Montgomery-Åsberg Depression Rating Scale (MADRS) (11) was shown to be the most accurate measure of depression severity (8) and to have the closest agreement with assessment by experienced clinicians (12), and thus it was used as the primary outcome measure. Since improvement in depression severity with treatment is a matter of degree rather than a yes or no phenomenon and there is no indication of bimodal endpoint outcomes (13), continuous measures, such as percentage improvement, are preferable to cutoff-based dichotomous measures of response or remission (14). Percentage change in MADRS scores from baseline to week 12 was selected instead of absolute change, since it shows high correlation with the end-of-treatment score (r=0.84) and low correlation with baseline severity (r=–0.06) and better reflects clinicians' impression of improvement (15). Since the commonly used last-observation-carried-forward procedure introduces bias in the presence of missing data (16), we imputed missing week-12 measurements using the best unbiased estimate from a mixed-effect linear regression model, with fixed linear and quadratic effects of time and random effects of individual and center of recruitment, which are relatively free of bias under a wider range of assumptions (4, 16). Finally, we adjusted for age and recruitment center because outcome was associated with age and varied across centers of recruitment. This adjustment minimized the risk that minor genetic differences as a result of stratification of European populations could be spuriously associated with outcome.

Genotyping

DNA was extracted from blood samples provided by 795 participants and collected in ethylenediaminetetraacetic acid (17). A total of 727 samples that were available in sufficient quantity and quality were sent to the Centre National de Genotypage (Evry Cedex, France) and genotyped using the Illumina Human610-quad bead chip (Illumina, Inc., San Diego). This chip assays more than 610,000 single nucleotide polymorphisms (SNPs) and copy number variant markers selected to provide a comprehensive coverage across populations, and it captures the majority of known common variations in the human genome, based on HapMap (release 23).

Statistical Analysis

After quality control (see the data supplement accompanying the online version of this article), the association between genotypes and adjusted percentage improvement in MADRS scores was performed using linear regression analyses under an additive genetic model, with the first four principal components from population structure analysis used as covariates to control for population stratification (see the data supplement), implemented in PLINK (18). The following four tests were performed: 1) association between genotype and outcome across the whole sample to identify genetic variants associated with improvement in depression severity irrespective of the type of antidepressant treatment; 2) association between genotype and outcome among subjects treated with escitalopram to identify genetic variants associated with the outcome of treatment with a serotonin reuptake inhibitor; 3) association between genotype and outcome among subjects treated with nortriptyline to identify genetic variants associated with the outcome of treatment with a noradrenaline reuptake inhibitor; and 4) interaction between genotype and drug to identify genetic variants that differentially predict the outcome of treatment with escitalopram and nortriptyline. The quality of control for population stratification effects was checked using quantile-quantile plots and genomic control lambda (19).

Frequentist and Bayesian approaches were used to assess the likelihood of results being true or false. To account for multiple testing, genome-wide significance was set at p<5×10^–8 (20, 21). However, since many markers associated with important individual differences in outcome may be identified at more modest levels of significance, all associations detected at p<5×10^–6 are reported as suggestive findings of interest. To estimate the posterior probability of true positive findings in the context of multiple nonindependent tests, false discovery rate q values were calculated, using the Benjamini and Hochberg step-up procedure (22), which have been shown to retain desirable properties in genome-wide studies and can be interpreted as posterior probabilities of no association at a given locus (23). We included q values alongside p values, and all associations with a q value ≤0.1 in the primary analyses are reported.

The pharmacogenetic analyses had statistical power to detect clinically significant pharmacogenetic effects explaining ≥5% of variance but were expected to miss a large proportion of moderate or weak effects (see the online data supplement).

Exploration of Candidate Genes

We extracted results for 2,801 markers in the coding regions and 20 kilo-base-pair upstream and downstream sequences of 72 candidate genes identified in a literature review (2) as encoding proteins involved in antidepressant action, previously reported to be associated with the outcome of antidepressant treatment or with mental illness in a genome-wide study (see Table 1 in the online data supplement for a list of candidate genes and rationale for their inclusion). Results from these candidate genes that remained significant after correction for the number of markers within a gene are reported.

Results

Genotyping Quality Control and Sample

Illumina Human610-quad genotyping was obtained for 727 subjects. Of the 550,337 SNPs with a minor allele frequency >0.01, a total of 539,391 (98%) were at least 99% complete and retained for analyses. Individual quality control is summarized in Figure 1. There were no sex mismatches, but three individuals with ambiguous genotypic sex were excluded. Three individuals representing outliers on heterozygosity were also excluded. One individual in each of six pairs of related individuals (three first- and three second-degree pairs of relatives) was retained for further analyses. Four individuals representing outliers with genotyping completeness <95% were excluded. After excluding five population structure outliers, 706 of the 811 participants were retained for the main analyses, with a mean genotyping completeness of 99.82%.

The clean sample of unrelated subjects of European ancestry passing quality control and included in genome-wide pharmacogenetic analyses is described in Table 1. These 706 included subjects did not differ from the other 105 GENDEP participants on sex, medication history, depression severity, or outcome. There was a difference in age, with included individuals being younger (mean age: 42.0 years [SD=11.7] versus 45.0 years [SD=11.6]; t=2.50, p=0.01), which may have been the result of chance, given the number of comparisons, and should not have influenced results because the outcome was adjusted for age.

Enlarge table

Genome-Wide Pharmacogenetic Analysis

Whole sample

The first linear regression analysis tested the association between 539,391 SNP markers and the continuous outcome of age- and recruitment center-adjusted percentage change in MADRS scores under an additive genetic model in the 706 individuals treated with either escitalopram or nortriptyline. A quantile-quantile plot showed a uniform distribution of p values (see the online data supplement), and a genomic control lambda coefficient of 1.007 confirmed no inflation of test statistics. No association was identified at a genome-wide level of significance. The seven strongest associations had suggestive levels of significance and a false-discovery q value <0.1 (Table 2). As seen in Figure 2, these associations were grouped in two intergenic regions on chromosomes 1 and 10. Both associated regions were contained within previously described common copy number variants (24), and the minor alleles were associated with worse response to both antidepressants.

Enlarge table

The three associated markers on chromosome 1 lie within a 163-kilo-base-pair duplication, which was reported for 98 of 270 HapMap subjects (24). Another six SNPs within this region, all with a minor allele frequency of approximately 0.3, could be imputed with high certainty and were associated with treatment response at similar levels of significance as the genotyped markers (see the online data supplement). Although no known genes overlap these markers, there are three expressed sequence tags that could represent two genes, which are potentially affected. None of the expressed sequence tags was homologous to a known gene, and thus it was not possible to infer their potential function.

The associated markers on chromosome 10 map within several structural variants, including two overlapping copy number variants of >200 kilo base pairs identified in 126/270 HapMap samples. Another 14 imputed SNPs with a minor allele frequency of 6%–7% were associated with treatment response (see the online data supplement). Again, no known genes overlapped these markers, but a single spliced expressed sequence tag (DA737884) was identified, homologous with the mixed-lineage leukemia-translocated MLLT10 gene (a chromatin modifier), located approximately 1 mega base pair downstream.

Escitalopram-treated subjects

A second linear regression analysis tested the association with the continuous outcome of age- and recruitment center-adjusted percentage change in MADRS scores for 394 individuals treated with escitalopram. The quantile-quantile plot and genomic control lambda coefficient of 1.011 showed no inflation of test statistics (see the online data supplement). No marker was associated at a genome-wide level of significance. A synonymous SNP in the coding region of the interleukin-11 gene, IL11, was associated with outcome at a suggestive level of significance (Table 2, Figure 2). This SNP creates new consensus sequences for the exonic splice enhancer factor 2/alternative splicing factor, which could favor an alternative IL11 splice variant (25). It was not in strong linkage disequilibrium with other genotyped or imputed markers.

Nortriptyline-treated subjects

A third linear regression analysis tested the association with the continuous outcome of age- and recruitment center-adjusted percentage change in MADRS scores under an additive genetic model in 312 individuals treated with nortriptyline. A quantile-quantile plot and genomic control lambda coefficient of 0.989 confirmed no inflation of test statistics (see the online data supplement). The intronic rs2500535 SNP within the uronyl 2-sulphotransferase (UST) gene on chromosome 6 was associated with outcome at a genome-wide level of significance and with a very small posterior probability of being a false positive (Table 2, Figure 2). Carriers of the minor A allele at rs2500535 experienced less improvement than carriers of the GG homozygote from week 4 onward. Imputation identified another 11 SNPs, with a minor allele frequency of 4%–6%, associated at genome-wide or suggestive levels of significance. Two of these imputed SNPs (rs2486404 and rs2500525) were predicted to be associated more strongly than the genotyped SNP, after taking into account the uncertainty of imputation. None of the associated markers were located in exons or regions predicted to influence gene function.

Another two genotyped markers associated at suggestive levels of significance (rs4651156 and rs9425322) were contained within an intron of the Ral guanine nucleotide dissociation stimulator-like 1 (RGL1) gene on chromosome 1 (Table 2, Figure 2), encoding a protein involved in G-protein signaling. Imputation identified another five SNPs in this region, with minor allele frequencies between 0.17 and 0.28, associated with response to nortriptyline (see the online data supplement). The genotyped marker rs9425322 causes a guanine-to-adenine transition and is located in an intronic region with strong mammalian conservation. Although this level of intronic conservation suggests regulatory function, the SNP does not directly affect known transcription factor binding sites, and the functional effect of this variant is unknown.

Genotype-drug interaction

A fourth linear regression analysis tested interactions between the type of antidepressant drug and genotypes in their effects on age- and recruitment center-adjusted percentage change in MADRS scores under an additive genetic model. The genomic control lambda coefficient was 0.991. The following two markers interacted with treatment at a suggestive level of significance: 1) rs11666579 in the solute carrier family 27 member 1 (SLC27A1) gene on chromosome 19, encoding a fatty acid transporter, and 2) rs1013696 in an intergenic region 11 kilo base pairs downstream from the nucleolar protein 4 (NOL4) gene on chromosome 18 (Table 2, Figure 2). One additional marker in strong linkage disequilibrium with rs11666579 was identified through imputation (see the online data supplement). The rs2500535 marker in the UST gene, associated with nortriptyline response, was also among the strongest interactions.

Pharmacogenetic Analysis

Candidate genes

Among the 2,801 markers in 72 candidate genes, none was significantly associated with outcome after correction for the number of comparisons (Bonferroni-corrected p value=1.8×10^–5). Pharmacogenetic associations that remained significant after correction for the number of markers within a gene are reported (Table 3). The strongest associations were 1) a marker in the norepinephrine transporter gene (SLC6A2) with overall response, 2) markers in a glutamatergic receptor gene (GRIA4) and the interleukin-6 (IL6) gene with response to escitalopram, and 3) markers in the leptin gene (LEP) and another glutamatergic receptor gene (GRIK1) with response to nortriptyline. A marker downstream of the HTR3A gene, encoding a serotonergic receptor, showed a gene-drug interaction because its minor C allele was associated with a poor response to escitalopram and a good response to nortriptyline.

Enlarge table

Discussion

This genome-wide analysis of the GENDEP project has shown that previously unexpected genomic regions may be more potent predictors of antidepressant response than functional candidate genes. Several genomic loci emerged that were more likely to be associated with response than not, after taking into account the extensive multiple testing performed. The genes encoding uronyl 2-sulphotransferase and IL11 as well as two intergenic regions containing common structural variants were revealed as potentially important for antidepressant action. To our knowledge, none of the markers associated with response in the GENDEP project has previously been reported in genome-wide studies of depression susceptibility (26, 27) or antidepressant response (28). Since the statistical power of the present investigation was limited, it is likely that a number of moderately strong associations have been missed, and integration of data from multiple large samples may reveal further important pharmacogenetic associations.

Irrespective of which antidepressant was used, the outcome of treatment was associated with polymorphisms in two regions on chromosomes 1 and 10. The nearest known gene to the associated locus on chromosome 1 (zinc finger protein 326 [ZNF326]) was 380 kilo base pairs away from the associated markers, and the nearest gene to the locus on chromosome 10 (plexin domain containing 2 [PLXDC2]) was 241 kilo base pairs away from the associated markers. Common copy number polymorphisms that may influence the expression of genes at longer distances through structural changes of the chromatin have been reported in both regions. This finding demonstrates the importance of exploring the whole genome rather than concentrating on coding regions of known genes and their flanking sequences.

Antidepressant-specific analyses revealed a strong genome-wide significant association between a marker in the UST gene and response to nortriptyline. Although the UST gene did not appear on any previous list of candidate genes, the function of uronyl 2-sulphotransferase in the nervous system suggests a mechanism for its involvement in antidepressant action. This enzyme is responsible for the production of oversulfated proteoglycans that are essential for neurogenesis and neuronal migration (29). A knockdown of UST led to severe disturbance of neuronal migration that was rescued by a UST cDNA (30). The apparent delayed onset of the pharmacogenetic effect on response after 4 weeks of treatment was consistent with a neurogenesis-related mechanism. It is presently unclear why such an effect should be specific to nortriptyline, since various types of antidepressants increase neurogenesis (31). Response to nortriptyline was also associated with markers in the RGL1 gene, involved in Ras and Ral signaling pathways of neuronal differentiation and outgrowth (32). Nearby markers in this gene have been associated with conduct disorder and hyperactivity (33). These findings require replication and investigation of the role of neurogenesis-related genes in antidepressant action.

Escitalopram response was predicted by a marker in the gene encoding interleukin-11. Although statistically less robust, this finding was further supported by an association of escitalopram response with a marker in the IL6 gene, an established candidate gene for depression (34) and a close homologue of IL11. Interleukins 6 and 11 are members of a family of neuropoietic cytokines, which are widely expressed in the brain (35) and can potently inhibit serotonin signaling by inducing raphe neurons to produce acetylcholine instead of serotonin (36). A differential expression of these cytokines, or alternative splicing suggested by an exonic splicing enhancer consensus sequence (25), may make some individuals vulnerable to a subtype of depression that may be related to environmental exposures (37) and poor response to serotonergic antidepressants (38). This molecular mechanism is consistent with the role of inflammation in a subtype of depression (39) and could explain the specific moderation by IL6 and IL11 of response to the serotonergic antidepressant escitalopram.

A comprehensive list of 72 candidate genes, including those encoding key proteins in monoaminergic and glutamatergic signaling, have produced findings consistent with chance. However, markers in several candidate genes, including IL6, HTR3A, and GRIA4, showed gene-wide significant pharmacogenetic associations that would have been reported as positive results in a candidate gene study. The implication of IL6 in response to serotonergic antidepressants was supported by the pharmacogenetic association with a functionally related, but genomically independent, IL11 gene. The HTR3A gene, encoding the serotonergic 5-HT_3A receptor, is an important candidate as effective antidepressants are antagonists of this receptor (40) and colocalize with the receptor on cell surfaces (41). These genes should be considered for exploration in future studies.

A genome-wide pharmacogenetic study faces a number of challenges. The strengths of the GENDEP project include a homogenous sample specifically recruited for a pharmacogenetic investigation, detailed and complete data characterizing response to antidepressants, tight control of population stratification with no overall inflation of test statistics, and high-quality genotyping enabling stringent quality control while retaining most genotyped individuals. On the other hand, our findings are limited by the absence of a placebo-comparison group, which makes it impossible to distinguish between genetic predictors of nonspecific improvement and effects that are common to both antidepressants used in the GENDEP project. Across the four analyses, only one pharmacogenetic association was established at a genome-wide level of significance. The predominance of negative findings could be the result of limited power to detect weak to moderate associations or of the absence of even moderately strong common genetic determinants of antidepressant drug response. A meta-analysis of multiple large samples will be needed to establish the generalizability of findings reported in the present study.

The pharmacogenetic analysis of GENDEP has brought up a number of promising results, including two regions containing structural variants, a genome-wide significant association with nortriptyline response, and a functionally plausible marker associated with response to escitalopram. These results demonstrate the feasibility of genome-wide pharmacogenetic analyses with well-characterized moderately large samples. As in all genetic studies, replications will be required before the value of any particular finding can be fully assessed, but genome-wide studies may thus be important steps in the elucidation of the genetic basis of pharmacological response.