Introduction

Major depressive disorder (MDD) is a common and debilitating disorder with pervasive impact on the quality of life for both sufferers and their families. Lifetime prevalence is estimated to be 15%,1, 2 and is consistently estimated to be twice as common in women as men.3, 4 MDD is associated with high morbidity, reflected in estimates of burden of disease and years lost in productivity,1, 5 and excess mortality from suicide6 and other causes. Our understanding of the etiology of MDD remains fragmented, despite wide-ranging research, but is the key to effective prevention and treatment.

MDD is familial, with heritability estimated to be 0.37 (95% confidence interval (CI) 0.31–0.42) and both early age of onset and recurrence of depression are associated with higher familial aggregation.7, 8 This implied genetic etiology has motivated studies designed to identify specific genetic variants associated with MDD. Results from genome-wide linkage, reviewed in Boomsma et al.,9 and candidate gene association studies10 have shown little consistency and hopes for new progress have spurred on a generation of genome-wide association studies (GWAS). Five GWAS for MDD have been published to date,11, 12, 13, 14, 15 each using control samples screened negative for MDD. None of these studies has identified variants that achieve genome-wide significance. Taken together, the results of these studies imply that specific genetic variants individually make very small contributions to the etiology of MDD. In this study, we present the largest GWAS for MDD to date, the MDD2000+ study comprising 2431 cases and 3673 screened controls. We compare our results with reports of the other published MDD GWAS and present a formal meta-analysis of our results with the two other largest studies (5763 cases and 6901 controls).

Materials and methods

Overview

The MDD2000+ project comprises a total of 2431 cases with MDD and 3673 screened controls from different sources and genotyped on different platforms (Tables 1 and 2). Samples were provided by the Queensland Institute of Medical Research (QIMR, Australia), The Netherlands Study of Anxiety and Depression (NESDA), The Netherlands Twin Registry (NTR), the University of Edinburgh (UK), and the Molecular Genetics of Schizophrenia study (controls only, US). Genotyping was conducted on different Illumina and Affymetrix platforms and because the overlap in genotyped single-nucleotide polymorphisms (SNPs) is limited, association analysis is based on a set of >1 M imputed SNPs. The number of SNPs for each analysis set (Table 1), represents genotyped SNPs surviving all quality control (QC) criteria that were used for imputation.

Table 1 Analysis set statistics: numbers of samples, numbers of SNPs, ages at interview and MDD onset
Table 2 Descriptive statistics by analysis set and sample source

Subject recruitment

QIMR

Study participants were adult twins and their families recruited through the Australian Twin Registry (http://www.twins.org.au). Only unrelated individuals were included in MDD2000+. All participants provided written informed consent under study protocols approved by the QIMR Human Research Ethics Committee. MDD cases were identified through psychiatric questionnaires, either the shortened Composite International Diagnostic Interview16 or the SSAGA-OZ interview instrument (a version of the Semi-Structured Assessment for the Genetics of Alcoholism17 modified for use in Australia), a comprehensive psychiatric interview designed to assess MDD and other psychiatric disorders17 according to DSM-IIIR18 and DSM-IV19 criteria. Structured interviews were administered by trained telephone interviewers, closely supervised by a clinical psychologist. Briefly, from 1988 to 1990, study participants were mailed an extensive health and lifestyle questionnaire, which included the shortened revised Eysenck personality questionnaire.20 Sum scores of 12 item responses in each personality domain resulted in quantitative scores for neuroticism. Between 1992–2000, an unselected subset of these participants were interviewed by telephone using the SSAGA-OZ.21 Over the period 1996–9 sibling pairs that were either concordant or discordant for extreme neuroticism scores (one sibling in the top or bottom decile, the other sibling in the top or bottom quintile) were recruited to complete the Composite International Diagnostic Interview, which provides DSM-IV19 lifetime diagnoses of MDD.22 Finally, some study participants completed the SSAGA-OZ telephone questionnaire in 2003–2007 as part of alcohol and nicotine dependence studies (the nicotine addiction genetics and Inter-related Project Grant studies described in Table 2 of Hansell et al.23). The Inter-related Project Grant/nicotine addiction genetics studies captured 28% of families who had already participated in earlier studies ascertaining those with either (a) large sibship size or containing a proband with either (b) nicotine dependence or (c) alcohol dependence. For this study, all cases met DSM-IV lifetime criteria for MDD. Screening items for mania were not consistent across interviews and screening items for psychosis were not included; the ability to assess accurately these less common criteria is difficult in large-scale community settings. Therefore, it is possible that a small number of individuals with a primary diagnosis of bipolar disorder or schizophrenia are included in the case group. If multiple cases were present within families then one case was selected in the following order of preference: age of onset <31 years, multiple episodes of depression, co-morbid anxiety disorders and high neuroticism score. Unrelated controls were selected as genotyped individuals from families in which no individuals qualified for diagnoses of MDD or anxiety disorders. If multiple controls were available from a family, the individual with the lowest neuroticism score was preferentially selected, otherwise an individual was selected at random.

NESDA and NTR

Additional MDD cases were selected from two parallel studies, NESDA (http://www.nesda.nl) and NTR (http://www.tweelingenregister.org); NTR also provided screened controls. These samples do not overlap with those included in a prior MDD GWAS12 but are drawn from the same parent studies (although a small number of related individuals are included, see meta-analysis section below). Details of the data collection methods are described elsewhere.9 Similar inclusion and exclusion criteria were used to select MDD cases from both the NESDA and NTR studies. Inclusion criteria were a lifetime diagnosis of DSM-IV MDD as determined by the Composite International Diagnostic Interview,16 age 18–65 years, and self-reported western European ancestry. Those not fluent in Dutch or with a primary diagnosis of schizophrenia or schizoaffective disorder, obsessive–compulsive disorder, bipolar disorder or severe substance use or dependence were excluded. All subjects provided written informed consent after study approval by the appropriate Institutional Review Boards.

University of Edinburgh

MDD cases were recruited through in- and out-patient services of psychiatric hospitals in Scotland and were tertiary referrals from primary care. All patients were interviewed by an experienced psychiatrist using the Schedule for Affective Disorders and Schizophrenia-Lifetime,24 supplemented by hospital case note review and information from informants. Final determination of MDD as the primary DSM-IV diagnosis was made by consensus of two psychiatrists. All cases had a lifetime history of recurrent MDD and IQ >70. The study was approved by the Central Office of Research Ethics Committees in Scotland and all subjects gave informed written consent for the collection of DNA samples for use in genetic studies.

Molecular Genetics of Schizophrenia controls

Briefly, random digit dialing was used to achieve a representative sample of individuals from the United States. Participants completed an online questionnaire including the short form Composite International Diagnostic Interview, supplemented by questions about schizophrenia, psychosis and bipolar disorder. Controls were required to never have met criteria for MDD. All subjects provided informed consent. These controls have been used in GWAS for multiple psychiatric disorders including MDD.13, 14

Genotyping and QC

Full details are given in the Supplementary File 1. All analysis sets were put through a common QC pipeline based on that used by the Psychiatric GWAS Consortium.25 SNPs were removed based on the following criteria: minor allele frequency (MAF) <0.01, Hardy–Weinberg equilibrium test P<1 × 10−6 in controls, missingness >0.01 or 0.02, missingness difference between cases and controls of 0.01 or 0.02, difference in frequency with HapMap of 0.07 or 0.15 (more stringent thresholds were applied when cases and controls were not genotyped concurrently, and were chosen to balance SNPs lost versus genomic control inflation factors). Samples were removed with missingness >0.02, if related to other samples or if identified as European ancestry outliers. Imputation was conducted in four analysis sets (I317, I370, I610, A6.0) in batches of 300 mixed cases and controls to a common set of SNPs present in HapMap3 CEU/TSI, consisting of 410 haplotypes, using Beagle 3.04.26, 27

Statistical analyses

Association analysis was conducted on 2431 MDD cases and 3673 screened controls using allelic dosages for imputed SNPs passing QC. The test of association was logistic regression, including the first three ancestry principal components (PC) and analysis set as covariates. The INFO score (ratio of observed to expected dosage variance), a measure of the quality of imputation of SNPs was used to interpret results. A second association analysis was conducted restricting cases to those with recurrent early onset MDD (REO, age of MDD onset <31 years with multiple episodes of MDD). Separate analyses were conducted by sex, as different prevalences of MDD between males and females could imply existence of sex-specific genetic risk variants.

Statistical power

Detailed power calculations are provided in the Supplementary File 1. The MDD2000+ sample affords >90% power to detect an associated variant with MAF 0.36 and genotype relative risk (GRR) 1.33 (the median values across all published complex disease/trait GWAS, with P<5 × 10−8).28

Comparison of association results with other studies including meta-analysis

We compared our results with those published for other GWAS of MDD. In particular, we conducted a formal genome-wide meta-analysis for autosomal SNPs from MDD2000+, GAIN–MDD12 and UK15 studies (that is, the three largest MDD GWAS). The GAIN–MDD sample, after removing individuals related to Dutch cases or controls of MDD2000+, comprised 1696 MDD cases and 1634 screened controls and the UK sample comprised 1636 MDD cases and 1594 screened controls. Meta-analysis29 on the logistic regression results from each study, restricting the imputed SNPs to those with INFO scores between 0.8 and 1.1.

Gene-based tests

To determine whether any genes harbored an excess of SNPs with small P-values, we undertook genome-wide gene-based tests which account for both gene length and linkage disequilibrium between SNPs using VEGAS.30 SNPs were allocated to one or more autosomal genes using gene boundaries ±50 kb. We investigated 183 candidate genes for MDD.12, 13, 14, 15, 31

Results

Of 1,251,157 imputed SNPs, 1,079,979 (86%) had MAF>0.01 and INFO>0.8. The quantile–quantile plots of the observed versus expected −log(P) from the four (total, male, female and REO) association analyses are presented in Figure 1 and Supplementary File 1. The genomic control λ (the median χ2 association statistic divided by the median expected under the null) shows no evidence for inflation of the test statistics. For example, for the full analysis, λ=1.05 (and standardized to a sample size of 1000, λ1000=1.02). There were no genome-wide significant results with P<5 × 10−8. The Manhattan plot of association P-values for all subjects is presented in the Supplementary File 1 and Table 3 shows associations with <1 × 10−5 in any analysis. Regional association and forest plots are provided in Supplementary File 2. SNPs associated with P<0.001 in any analysis are provided in Supplementary File 3.

Figure 1
figure 1

Quantile–quantile plots for the association analyses of (i) all cases and controls and (ii) recurrent early onset (iii) meta-analysis.

PowerPoint slide

Table 3 Regions containing at least one SNP with a P-value of P<10−5 in any one of the four association analyses—total sample (all), males only, females only and recurrent early onset (REO) MDD ordered by chromosomal position

Comparison of association results with other studies including meta-analysis

We compared the MDD2000+ results with other published studies. Muglia et al.11 listed 27 SNPs with P<1 × 10−5; 25 SNPs overlapped with our analysis but none had P<0.05 and association with the same allele. Of the top 200 SNPs in one or both of their study samples, 90% were analysed in our study but 5% had P<0.05, consistent with chance expectations. As the MDD2000+ controls were also used in GenRed13 and STARD,14 we did not make formal comparisons of our results, but we note that 5% of their top SNPs (P<1 × 10−5) had P<0.05 in MDD2000+. The SNP rs2251219 identified as a mood disorder risk factor in a GWAS meta-analysis (three bipolar studies plus the GAIN–MDD study)32 but was not associated in MDD2000+ (P=0.51, MAF=0.40, odds ratio (OR)=0.97, 95% CI 0.90–1.05). PCLO was the top finding in the GAIN–MDD study,12 replicated in Australian (QIMR)12 and Dutch33 samples, but we found no evidence for association with PCLO variants (for example, P=0.51 for rs2522833), despite partial overlap (562 cases and 264 controls) in the QIMR samples used here and in the original replication sample (see Supplementary File 1 for details). Meta-analysis results of the three largest MDD GWAS—427362 SNPs in MDD2000+, GAIN–MDD and the UK studies—yielded λ1000=1.01 (Figure 1iii). In all, 19 SNPs gave P<5 × 10−5 compared with 21 expected under the null hypothesis (ignoring linkage disequilibrium) and are listed in the Supplementary File 1.

Gene-based tests

MDD2000+ SNPs mapped to 18 454 genes. We list all genes with P<1 × 10−4 and all lie within linkage regions identified in the SLEP34 database of research findings for psychiatric disorders. Gene ontology annotations and names are provided for these genes (Table 4). Three are plausible candidate genes for MDD (GAL, ADCY3, PDK4) and others have a role in ion binding (ZDHHC19, ZNF83). None of the genes listed in Table 4 shows any evidence for association in the gene-based test applied to GAIN–MDD or the UK MDD GWAS results, although association at P<0.001 is retained for GAL in the gene-based test applied to the meta-analysis results.

Table 4 Top ten ranked genes from gene-based test for genes ±50 kb from start and stop positions

Of the 183 candidate genes investigated in other GWAS for MDD, we could test for association with 180 (Supplementary File 4). None were associated with P<0.00028 (0.05/180). The six most associated were IL10, OPRM1, HTT, HTR1B, GRIN1 and CACNA1C. Of these, OPRM1 was reported to have P<0.05 in the GAIN–MDD gene-based test and CACNA1C was top ranked in the GenRed13 study; the other studies all used a gene test based on best single associated SNP corrected for gene length. In this study, CACNA1C ranked 808 out of 18 454 (4th percentile). Our top associated SNP in CACNA1C is rs98545 (P=0.0019, OR=0.83, MAF=0.15) which is 387 kb from, and in near linkage equilibrium (r2=0.06) with, rs1006737 (P=0.020, OR=1.10, MAF=0.35), the SNP first identified in genome-wide association studies of bipolar disorder;35 an interaction between these SNPs was not significant (P=0.21). The 3rd ranked of the candidate genes (Huntingtin, HTT) was included in the candidate list because mood disorders are characteristic of presymptomatic carriers of the Huntington's disease trinucleotide repeat polymorphisms.31 HTT ranked in the top 0.7% of all genes; the SNP showing highest individual association within the gene is rs363099 (also known as rs4690074, P=4.7 E−04, OR=0.86, MAF=0.31) located in exon 29, 85 kb from the exon 1 CAG repeat.

Discussion

Study results

The MDD2000+ study is the largest GWAS for MDD reported to date, comprising 2431 cases with MDD and 3673 screened controls. Our analysis was of >1 M SNPs but none achieved genome-wide significance. Indeed, the quantile–quantile plots (Figure 1) show that the distribution of observed associations closely follows that expected under the null hypothesis. These results are consistent with other GWAS of MDD.11, 12, 13, 14, 15 Furthermore, comparison of our results with these other studies including formal meta-analysis of results from MDD2000+ and the next two largest studies12, 15 (a total of 5763 cases and 6901 controls) showed little evidence for replication between studies.

To determine whether any genes harbored an excess of low-associated variants, we undertook gene-based tests which account for both gene length and linkage disequilibrium between SNPs. Our top associated genes were not replicated in the gene-based test when applied to the GAIN–MDD or UK MDD results, although some genes retained low association P-values when applied to the meta-analysis results (Table 4). Despite this there are some suggestive results worthy of note. Among the top ranked genes is GAL, which encodes the neuropeptide galanin, proposed by Weiss et al.36 as having an important role in MDD. Their hypothesis, based on animal models is that GAL released in the ventral tegmentum inhibits the activity of dopaminergic cells resulting in decreased motor activity and anhedonia. GAL is a regulator of brain serotonin and 5-HT1A receptor-mediated transmission,37 agonists of GAL receptors have been proposed as potential drug targets for MDD.38 GAL also has an important role in the hippocampal processing of cognition.39 Our most associated SNP in GAL (rs2156464, P=2.7 × 10−5, OR=1.24, CI=1.12–1.37, MAF=0.19) lies in the same haplotype block as an association reported for panic disorder.40 Another top ranked gene (ADCY3) encodes the enzyme adenylate cyclase 3 that catalyses synthesis of cyclic adenosine monophosphate, the association underpinned by SNP rs2384061 (P=6.3 × 10−6, OR=1.20, CI=1.11–1.29, MAF=0.42, Table 3). ADCY3 is an plausible candidate gene because depressed patients display reduction in platelet ADCY activity.41 Perhaps most interesting is that ADCY and GAL activity are inter-related42 and ADCY3 binds CACNA1C (calcium channel, voltage dependent, L type, α-1C subunit).43 A SNP (rs1006737, intron 3) within the gene CACNA1C is associated with bipolar disorder44 schizophrenia and MDD.45 In this study, CACNA1C ranked in the top 4% of genes and we replicated the association with rs1006737 (P=0.020, OR=1.10, CI=1.01–1.19, MAF=0.34), but the top associated SNP in CACNA1C was rs98545 (P=0.0019, OR=0.83, CI 0.75–0.94, MAF=0.15, intron 29). rs98545 lies 387 kb from rs1006737; they are in near linkage equilibrium (r2=0.06) and in our data their effects are additive in a logistic regression. Neither rs1006737 nor rs98545 were genotyped in the UK MDD or GAIN–MDD studies; imputation of GAIN–MDD genotypes showed weak association with rs1006737 (P=0.070, OR=1.10, CI=0.99–1.21, MAF=0.32) but not with rs98545 (P=0.52). We cannot exclude that the association of rs98545 has occurred by chance in our sample. Regional association plots for ADCY3, GAL and CACNA1C can be found in Supplementary File 2.

Despite the large sample size of MDD2000+, a number of limitations may reduce the realized power for detection of association. First, our study combines data sets genotyped on different primary platforms, and the power of our sample will vary between SNPs depending on the linkage disequilibrium structure between genotyped and imputed SNPs. Second, some of our cases and controls were genotyped separately because of economic constraints. We were fortunate that we had >300 QIMR cases genotyped on both the Affymetrix and Illumina platforms, which we used extensively for QC. Third, when combining cases and controls genotyped as separate sets, it was necessary to impose more stringent QC constraints than for cases and controls genotyped on the same platform. These constraints were successful in removing artefactual differences between the case–control sets, but could also have eliminated true associations. Fourth, in MDD2000+ with the aim of maximizing sample size we have combined the clinical cases from the University of Edinburgh cohort with the community samples from QIMR and NTR. The extent to which this is important depends on the unknown genetic etiology of MDD. Lastly, 26% of both cases and controls contributed by the QIMR sample (which makes the largest contribution in this study) were from families ascertained on the basis of a sibling with nicotine dependence, as reflected by the high proportion of smokers in the QIMR samples. Smoking is associated with depression and if some genes are involved in smoking behavior and MDD, then the power to detect genetic variants for MDD may be reduced.

Implications

Our analysis and meta-analysis represent the largest and most powerful investigation into the genetic architecture of MDD to date. The lack of clear-cut evidence for association allows us to exclude genetic main effects that are common (MAF>0.1) with moderate to strong effect sizes (GRR>1.4). As noted earlier, the MDD2000+ sample alone possessed >90% power to detect the median MAF (0.36) and GRR (1.33) of GWAS associations published to date with P<5 × 10−8. Moreover, as we would expect to detect >20% of variants with OR 1.2 and MAF 0.2–0.7 (Supplementary File 1), and having detected none that are genome-wide significant, we can probably draw an even more stringent boundary on the risk allele architecture under a main effect model. These results are informative for researchers planning association studies of MDD.

These results are unsurprising in the light of results for a range of complex genetic disorders published since this study was conceived in 2007.46 Results from GWAS of psychiatric disorders47 have, on the whole, yielded replicated associations only for low prevalence/high heritability disorders, that is, autism, bipolar disorder and schizophrenia, with tobacco dependence being the exception.48 As a high prevalence/low heritability disorder, the challenge for MDD was always going to be harder. From Yang et al.,49 we estimate that sample sizes 2.4-fold greater are needed for GWAS of MDD (prevalence 0.15) compared with schizophrenia (prevalence 0.007) to identify a variant that explains the same proportion of phenotypic variance to liability of these disorders (Supplementary File 1). This result reflects the smaller mean difference in phenotypic liability between cases and controls for MDD compared with schizophrenia. Hospital-based MDD cohorts may represent a more extreme phenotype, with lower prevalence and higher heritability.50 Using a prevalence for such clinical samples to be 0.06 (the average across sexes)51 still requires a sample size 1.8 times greater for a case–control study of MDD compared with one for schizophrenia.

Further, for main effect models, as the heritability on the liability scale (h2) is a function of the number of variants (n), their frequencies (qi) and their effect sizes (ai),

(with total variance of liability of 1), then the lower heritability of MDD (0.37, although this may be higher in clinical samples)51 compared to schizophrenia (0.81), must be explained by fewer risk alleles, lower risk allele frequencies and/or smaller effect sizes. It seems most plausible that MDD might have smaller effect sizes compared with schizophrenia (except in the unlikely event that the number of loci for MDD is substantially fewer than for schizophrenia). We estimate that sample sizes 4–5 times those needed for schizophrenia studies may be needed for MDD to detect variants that explain an equal proportion of the known genetic variance (Supplementary File 1).

Study designs

It is clear that larger sample sizes are needed to achieve the power required to detect common variants of smaller effect, and that larger sample sizes will be needed for MDD compared with less prevalent but more heritable disorders. The prize is major: any replicated hit provides a priceless window into the etiology of this idiopathic disorder. However, alternative and complementary strategies are possible, particularly those that consider a broader range of genetic models and relax further assumptions of etiological homogeneity.

First, recognizing that effect size is the average effect of a variant across all genetic and environmental backgrounds,52 it may be possible to forgo larger sample sizes to concentrate on more homogenous subsets in which larger effect sizes may be detectable. There has been much debate about the genetic heterogeneity of MDD. Some genetic heterogeneity is consistent with the polygenic model; each affected individual could potentially harbor a different combination of risk variants which could imply a spectrum of phenotypic symptoms, with those presenting similar symptom profiles carrying more similar profiles of risk alleles, consistent with symptom sharing of closely related individuals.53 For this reason, most GWAS for MDD have prioritized genotyping of the less prevalent and more heritable recurrent early onset MDD, but even so genome-wide significant associations have been elusive. Stratification by sex has also not yielded consistent results. Other possibilities certainly exist but phenotypes may be unavailable (for example, accurate delineation of typical versus atypical MDD)54 or prohibitively expensive (for example, magnetic resonance imaging). Use of quantitative scores of severity and reliability may be the best strategy to balance sample size with phenotype definition,55 as MDD symptom profiles are consistent with a quantitative rather than dichotomous liability.56

Second, in addition to achieving a more homogenous genetic background, larger estimated effect sizes may come from selecting more homogenous environmental backgrounds. The low heritability of liability for MDD implies an important role for environmental risk factors. Although genotype × environment interaction cannot explain the so-called ‘missing heritability’,52 it can contribute to small effect sizes. Although genotype × environment studies are conceptually attractive, the lessons learned from the most studied genotype × environment hypothesis for MDD (5HTTLPR and stressful life event) are sobering.57 However, this broad conclusion has been challenged58, 59 with Uher and McGuffin58 arguing that replication has been achieved when stressful life event were recorded objectively, temporally and before the onset to MDD. More mileage might be gained by considering MDD in the context of more clearly defined exposures, such as pregnancy, in which women with perinatal and post-partum MDD have been exposed to a similar event.

Whichever way we look at it, and whether risk variants are common or rare, it seems that the challenge for MDD will be much harder than for the less prevalent more heritable psychiatric disorders. Larger samples are required whether we attempt to identify associated variants with small effect across average backgrounds or attempt to enhance detectable effects sizes by selection of homogeneity of genetic or environmental background. In the long-term, a greater understanding of the etiology of MDD will require large prospective, longitudinal, uniformly and broadly phenotyped and genotyped cohorts that allow the joint dissection of the genetic and environmental factors underlying MDD.