Study population

The UK Biobank cohort is a population-based cohort consisting of 501,726 individuals recruited at 23 centres across the United Kingdom. Genotypic data was available for 488,380 individuals and was imputed with IMPUTE4 using the HRC reference panel31 to identify ~39 M variants for 487,409 individuals32. We excluded 79,990 individuals that were outliers based on heterozygosity, had a variant call rate <98%, or were not recorded as “white British”. We excluded a further 131,790 related individuals based on a shared relatedness of up to the third degree using kinship coefficients (>0.044) calculated using the KING toolset33; we then subsequently added back in one member of each group of related individuals by creating a genomic relationship matrix and selected individuals with a genetic relatedness less than 0.025 with any other participant (n = 55,745). We removed variants with a call rate < 98%, a minor allele frequency <0.01, deviation from Hardy–Weinberg equilibrium (P < 10−6) or an imputation accuracy (Info) score < 0.1, leaving a total of 7,666,894 variants for 331,374 individuals.

Extensive phenotypic data were collected for UK Biobank participants using health records, biological sampling, physical measures and touchscreen tests and questionnaires. We used three definitions of depression in the UK Biobank sample, which are explained in greater depth in Supplementary Note 2 and are described below. For the three UK Biobank depression phenotypes, we excluded participants who were identified with bipolar disorder, schizophrenia or personality disorder using self-declared data, touchscreen responses (per Smith, et al.10) or ICD codes from hospital admission records; and participants who reported having a prescription for an antipsychotic medication during a verbal interview. Further exclusions were applied to all control individuals if they had a diagnosis of a depressive mood disorder from hospital admission records, had reported having a prescription for antidepressants or had self-reported depression (see Supplementary Note 2 for full phenotype criteria).

This research has been conducted using the UK Biobank Resource—application number 4844. The UK Biobank study was conducted under generic approval from the NHS National Research Ethics Service (approval letter dated 17th June 2011, Ref 11/NW/0382). All participants gave full informed written consent.

Broad depression phenotype

The broadest phenotype (broad depression) was defined using self-reported help-seeking behaviour for mental health difficulties. Case and control status was determined by the touchscreen response to either of two questions: “Have you ever seen a general practitioner (GP) for nerves, anxiety, tension or depression?” (UK Biobank field: 2090) or “Have you ever seen a psychiatrist for nerves, anxiety, tension or depression?” (UK Biobank field 2010). Caseness for broad depression was determined by answering “Yes” to either question at either the initial assessment visit, at any repeat assessment visit, or if there was a primary or secondary diagnosis of a depressive mood disorder from linked hospital admission records (UK Biobank fields: 41202 and 41204; ICD codes: F32—Single Episode Depression, F33—Recurrent Depression, F34—Persistent mood disorders, F38—Other mood disorders and F39—Unspecified mood disorders). The remaining respondents were classed as controls if they provided “No” responses to both questions during all assessments that they participated in. This provided a total of 113,769 cases and 208,811 controls (n total = 322,580, prevalence = 35.27%) for the broad depression phenotype. Individuals classed as cases for the broad depression phenotype potentially included individuals seeking treatment for personality disorders.

Probable MDD phenotype

The second depression phenotype (probable MDD) was derived from touchscreen responses to questions about the presence and duration of low mood and anhedonia, following the definitions from Smith et al.10, whereby the participant had indicated that they were “depressed/down for a whole week (UK Biobank field: 4598); plus at least 2 weeks duration (UK Biobank field: 4609); plus ever seen a GP or psychiatrist for nerves, anxiety or depression” (UK Biobank fields: 2090 and 2010), or “ever anhedonia for a whole week (UK Biobank field: 4631); plus at least 2 weeks duration (UK Biobank field: 5375); plus ever seen a GP or psychiatrist for nerves, anxiety, or depression” (UK Biobank fields: 2090 and 2010). Cases for the probable MDD definition were supplemented by diagnoses of depressive mood disorder from linked hospital admission records (UK Biobank fields: 41202 and 41204) as per the broad depression phenotype. There were a total of 30,603 cases and 143,916 controls (n total = 174,519, prevalence = 17.54%) for the probable MDD phenotype.

There were 66,176 individuals that were classed as controls for probable MDD, who were classed as cases for broad depression. This was due to stricter conditions for classification as a case for probable MDD compared to the broad depression phenotype, and potentially removes individuals seeking treatment for personality disorders.

ICD-coded MDD phenotype

The ICD-coded MDD phenotype was derived from linked hospital admission records (UK Biobank fields: 41202 and 41204). Participants were classified as cases if they had either an ICD-9/10 primary or secondary diagnosis for a depressive unipolar mood disorder (ICD codes: F32—Single Episode Depression, F33—Recurrent Depression, F34—Persistent mood disorders, F38—Other mood disorders and F39—Unspecified mood disorders). ICD-coded MDD controls were participants who had linked hospital records, but who did not have any diagnosis of a mood disorder and were not probable MDD cases. There were 8,276 cases and 209,308 controls (n total = 217,584, prevalence = 3.80%) for the ICD-coded MDD phenotype.

There were no individuals classed as a case for ICD-coded MDD who were classed as a control for either the broad depression or probable MDD phenotypes. There were no individuals classed as controls in ICD-coded MDD who were classed as a case for probable MDD. However, there were 53,491 control individuals for ICD-coded MDD who were classed as a control for broad depression, which reflects the potentially looser definition stipulated under the broad depression phenotype. Of the 8,276 individuals classed as a case using the ICD-coded MDD definition, 7,471 (90.3%) also reported seeing either a GP or psychiatrist for nerves, anxiety, or depression (UK Biobank fields: 2090 and 2010), or had been depressed/down or suffered from anhedonia for a whole week (UK Biobank field: 4598).

A cross tabulation of case and control status between each pair of phenotypes is provided in Supplementary Table 7.

Association analysis

We performed a linear association test to assess the effect of each variant using BGENIE v1.132:

$$\begin{array}{l}{\mathbf{y}} = {\mathbf{{X}\beta }} + {\mathbf{\varepsilon }}_1\\ {\hat{\mathbf y}} = {\mathbf{{X}\beta }}\\ ({\mathbf{y}} - {\hat{\mathbf y}}) = {\mathbf{G}}{\mathrm{b}} + {\mathbf{\varepsilon }}_2\end{array}$$

where \({\mathbf{y}}\) was the vector of binary observations for each phenotype (controls coded as 0 and cases coded as 1). \({\mathbf{\beta }}\) was the matrix of fixed effects, including sex, age, genotyping array, and 8 principal components. \({\mathbf{X}}\) was the corresponding incidence matrices. \(({\mathbf{y}} - {\hat{\mathbf y}})\) was a vector of phenotypes residualized on the fixed effect predictors, \({\mathbf{G}}\) was a vector of expected genotype counts of the effect allele (dosages), \({\mathrm{b}}\) was the effect of the genotype on residualized phenotypes, and \({\mathbf{\varepsilon }}_1\) and \({\mathbf{\varepsilon }}_2\) were vectors of normally distributed errors.

Genome-wide statistical significance was determined by the conventional threshold of P < 5 × 10−8. To determine significant variants that were independent, the clump command in Plink 1.90b411 was applied using --clump-p1 1e-4 --clump-p2 1e-4 --clump-r2 0.1 --clump-kb 3000, mirroring the approach of Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium., et al.3. Therefore, variants that were within 3 Mb of each other and shared a linkage disequilibrium greater than 0.1 were clumped together, and only the most significant variant was reported.

Due to the complexity of major histocompatibility complex (MHC) region, an approach similar to that of The Schizophrenia Psychiatric Genome-Wide Association Study Consortium34 was taken, and only the most significant variant across that region was reported. To obtain odds ratios for those variants associated with depression, a logistic regression was conducted in Plink 1.90b411 using the same covariates as the linear model above. Regional visualisation plots were produced using LocusZoom35. Linkage Disequilibrium Score regression (LDSR)13 was used to determine whether there was elevation of the polygenic signal due to population stratification, by examining the intercept for evidence of significant deviation ( ± 1.96 standard error) from 1. The genomic inflation factor (\(\lambda _{\mathrm{GC}}\)) was also reported for each phenotype.

Replication cohort and meta-analysis

We sought to replicate those variants within UK Biobank that were identified as significantly associated (P < 5 × 10−8) with depression. To conduct the replication, we used the associated analysis results from the discovery sample from the 23andMe cohort4 (n = 307,354, cases = 75,607 and controls = 231,747). We firstly examined if the effect of the A1 allele was in the same direction across both UK Biobank and 23andMe. Secondly, we determined whether the variant within 23andMe was significant for the 17 variants (α = 0.05 / 17; P < 0.0029). Additionally, we used Metal36 to conduct an inverse variance-weighted meta-analysis of these 17 significant variants in UK Biobank and 23andMe.

Heritability of depression

We used GCTA-GREML12 to estimate the SNP-based h2 for each depression phenotype within each recruitment centre and each region. The specified population prevalence was matched to that observed across the whole cohort for each phenotype to allow transformation from the observed to the underlying liability scale. The recruitment centres within each region is provided in Supplementary Table 4. LDSR13 was used to provide a whole sample SNP-based estimate of the h2 on liability scale of the phenotypes using the whole-genome summary statistics obtained by the association analyses. The summary statistics for each phenotype were filtered using the default file, w_hm3.snplist, with the default LD Scores computed using 1000 Genomes European data (eur_w_ld_chr) used as a reference.

Genetic correlations with depression

To calculate the genetic correlations between the three depression-related phenotypes in UK Biobank LDSR13 was used. We also tested that hypothesis that the phenotypes were perfectly correlated with one another (r g = 1):

$$P-{\mathrm{value}} = 2 \times \mathrm{pnorm}\,\left( { - \mathrm{abs}\left( {\frac{{1 - r_{\mathrm{g}}}}{{\mathrm{standard}}\,{{\rm error}}}} \right)} \right)$$

Genetic correlations were also calculated between each of the three phenotypes and 235 other behavioural and disease related traits using LD Hub14. P-values were false discovery rate (FDR) adjusted using the Benjamini and Hochberg37 approach.

Gene- and region-based analyses

Two downstream analyses of the results were conducted using MAGMA18 (Multi-marker Analysis of GenoMic Annotation) by applying a multiple regression model, which incorporates linkage disequilibrium between markers and accounts for multi-marker effects, to the results of our association analyses. Estimates of the linkage disequilibrium between markers were obtained from the European populations sequenced as part of the 1,000 Genomes Project (phase 1, release 3)38. In the first downstream analysis, a gene-based analysis was performed for each phenotype using the results from our GWAS. The NCBI 37.3 build was used to determine the genetic variants that were ascribed to each gene. A total of 18,033 genes were assessed for an association with each depression phenotype with Bonferroni correction used to determine significance (α = 0.05 / 18,033; P < 2.77 × 10−6).

In the second downstream analysis, a region-based analysis was performed for each phenotype. To determine the regions, haplotype blocks identified by recombination hotspots were used as described by Shirali, et al.39 and implemented in an analysis of MDD by Zeng, et al.40 for detecting causal regions. Block boundaries were defined by hotspots of at least 30 cM per Mb based on a European subset of the 1,000 genome project recombination rates. This resulted in a total of 8,308 regions being analysed using the European panel of the 1,000 Genomes data (phase 1, release 3) as a reference panel to account for linkage disequilibrium. A genome-wide significance threshold for region-based associations was calculated using the Bonferroni correction method (α = 0.05 / 8,308; P < 6.02 × 10−6).

Gene-set pathway analysis

We used the results obtained from our gene-based analysis to conduct a further gene-set pathway analysis to test for gene enrichment within 5,917 biological classes using MAGMA18. The gene-set pathways were obtained from the gene-annotation files provided by the Gene Ontology (GO) Consortium (http://geneontology.org/)19, which are derived from the Molecular Signatures Database (MSigDB) v5.220. To correct for multiple testing, the default of 10,000 permutations was applied. The estimate of the effect size (beta) reflects the difference in association between genes in the gene set and genes outside the gene set from fitting a regression model to the data.

Tissue enrichment analysis

The DEPICT package21 was used to investigate whether specific tissues were enriched using the genome-wide association study results obtained for each phenotype. The recommended P-value threshold of <10−5 was applied to the variants used to assess enrichment. DEPICT incorporates data from 37,427 human gene expression microarrays to determine whether genes in associated loci for each phenotype were significantly expressed in any of 209 tissue/cell type annotations.

eQTL identification

The online GTEx portal (https://www.gtexportal.org/home/) was used to determine whether any of the genome-wide significant variants for each phenotype were eQTL9.

Code availability

All code used to conduct the research detailed in this manuscript is available on request from the corresponding author.

Data availability

The summary statistics relating to each of the three depression phenotypes assessed in this study are available on the Edinburgh DataShare website, via the following link: https://doi.org/10.7488/ds/2314

The raw genetic and phenotypic data that support the findings of this study are available from UK Biobank but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are, however, available from the authors upon reasonable request and with permission of UK Biobank (http://www.ukbiobank.ac.uk/).