Genome-wide association study identifies 30 loci associated with bipolar disorder.

Bipolar disorder is a highly heritable psychiatric disorder. We performed a genome-wide association study (GWAS) including 20,352 cases and 31,358 controls of European descent, with follow-up analysis of 822 variants with P < 1 × 10-4 in an additional 9,412 cases and 137,760 controls. Eight of the 19 variants that were genome-wide significant (P < 5 × 10-8) in the discovery GWAS were not genome-wide significant in the combined analysis, consistent with small effect sizes and limited power but also with genetic heterogeneity. In the combined analysis, 30 loci were genome-wide significant, including 20 newly identified loci. The significant loci contain genes encoding ion channels, neurotransmitter transporters and synaptic components. Pathway analysis revealed nine significantly enriched gene sets, including regulation of insulin secretion and endocannabinoid signaling. Bipolar I disorder is strongly genetically correlated with schizophrenia, driven by psychosis, whereas bipolar II disorder is more strongly correlated with major depressive disorder. These findings address key clinical questions and provide potential biological mechanisms for bipolar disorder.

Bipolar disorder (BD) is a severe neuropsychiatric disorder characterized by recurrent episodes of mania and depression that affect thought, perception, emotion, and social behaviour. A lifetime prevalence of 1-2%, elevated morbidity and mortality, onset in young adulthood, and a frequently chronic course make BD a major public health problem and a leading cause of the global burden of disease [1]. Clinical, twin and molecular genetic data all strongly suggest that BD is a multifactorial disorder [2]. Based on twin studies, the overall heritability of BD has been estimated to be more than 70% [3,4], suggesting a substantial involvement of genetic factors in the development of the disorder, although non-genetic factors also influence risk.

BD can be divided into two main clinical subtypes [5,6]: bipolar I disorder (BD1) and bipolar II disorder (BD2). In BD1, manic episodes typically alternate with depressive episodes during the course of illness. Diagnosis of BD2 is based on the lifetime occurrence of at least one depressive and one hypomanic (but no manic) episode. Although modern diagnostic systems retain the Kraepelinian dichotomy [7] between BD and schizophrenia (SCZ), the distinction between the two disorders is not always clear-cut, and patients who display clinical features of both disorders may receive a diagnosis of schizoaffective disorder-bipolar type (SAB). Likewise, in genetic studies BD and SCZ are usually treated separately, although recent epidemiological and molecular genetic studies provide strong evidence for some overlap between the genetic contributions to their etiology [2,8].

Recent genome-wide association studies (GWAS) in BD have identified a number of significant associations between disease status and common genetic variants [9-23]. The first large collaborative BD GWAS by the multinational Psychiatric Genomics Consortium (PGC) Bipolar Disorder Working Group comprised 7,481 BD patients and 9,250 controls and identified four genome-wide significant loci [9]. Three subsequent meta-analyses that included the PGC BD data [10,12,18] identified an additional 5 loci.

Estimates of the proportion of variance in liability attributable to common variants genome-wide (SNP-heritability) indicate that ~30% of the heritability for BD is due to common genetic variants [8]. To date, only a small fraction of this heritability is explained by associated loci, but results from other human complex traits suggest that many more will be identified by increasing the sample size of GWAS [24]. Here, we report the second GWAS of the PGC Bipolar Disorder Working Group, comprising 20,352 cases and 31,358 controls of European descent in a single, systematic analysis, with follow up of top findings in an independent sample of 9,412 cases and 137,760 controls. Some of our findings reinforce specific hypotheses regarding BD neurobiology; however, the majority of the findings suggest new biological insights.

RESULTS

GWAS of bipolar disorder (BD)

We performed a GWAS meta-analysis of 32 cohorts from 14 countries in Europe, North America and Australia (Supplementary Table 1A), totaling 20,352 cases and 31,358 controls of European descent (effective sample size 46,582). This is a large GWAS of BD, a 2.7-fold increase in the number of cases compared to our previous GWAS [9], and includes 6,328 case and 7,963 control samples not previously reported. We imputed variant dosages using the 1,000 Genomes reference panel, retaining association results for 9,372,253 autosomal variants with imputation quality score INFO > 0.3 and minor allele frequency ≥ 1% in both cases and controls. We performed logistic regression of case status on imputed variant dosage using genetic ancestry covariates. The resulting genomic inflation factor (λGC) was 1.23, 1.01 when scaled to 1,000 cases and 1,000 controls (λ1000) (Supplementary Figure 1). The LD Score regression intercept was 1.021 (se=0.010), and the attenuation ratio of 0.053 (se=0.027) was non-significant, indicating that the observed genomic inflation is indicative of polygenicity rather than stratification or cryptic population structure [25]. The LD-score regression SNP-heritability estimates for BD were 0.17-0.23 on the liability scale assuming population prevalence of 0.5-2%. See Supplementary Table 1A, Online Methods and Supplementary Note for sample and method details.

We find a marked increase in phenotypic variance explained by genome-wide polygenic risk scores (PRS) compared to previous publications (sample size weighted mean observed Nagelkerke’s R2 = 0.08 across datasets, liability scale R2=0.04, for p-threshold 0.01; Supplementary Figure 2 and Supplementary Table 2). Among the different datasets, we observed no association between the PRS R2 and: (i) the gender distribution of the BD cases (p=0.51); (ii) the proportion of cases with psychosis (p=0.61); (iii) the proportion with a family history of BD (p=0.82); or (iv) the median age of onset for BD (p=0.64). In our primary genome-wide analysis, we identified 19 loci exceeding genome-wide significance (P< 5×10−8; Table 1).

Table 1.

Genome-wide significant bipolar disorder risk loci

Locus Name[*1]	Lead SNP	CHR	BP	A1/A2	GWAS Meta-analysis			Follow-up samples		Combined
					Freq. A1	OR	P-value[*2]	OR	P-value[*3]	OR	P-value[*2]
A. Thirty loci with lead SNP P < 5×10−8 in combined GWAS+followup analysis
1,PLEKHO1	rs7544145	1	150,138,699	T/C	0.81	1.095	4.8E-07	1.064	0.010	1.085	4.8E-08
2,LMAN2L**	rs57195239	2	97,376,407	I/D	0.34	0.92	5.8E-09	0.96	0.030	0.93	3.8E-09
3,SCN2A	rs17183814	2	166,152,389	A/G	0.075	0.87	1.5E-07	0.89	0.0017	0.88	2.0E-09
4,[Intergenic]***	rs61332983	2	194,465,711	I/D	0.41	0.93	2.3E-08	0.95	0.0031	0.93	7.9E-10
5,TRANK1**	rs9834970	3	36,856,030	T/C	0.51	0.90	5.5E-14	0.98	0.15	0.93	5.7E-12
6,ITIH1**	rs2302417	3	52,814,256	A/T	0.49	0.92	4.9E-09	0.94	0.0012	0.93	6.6E-11
7,CD47	rs3804640	3	107,793,709	A/G	0.53	1.075	9.3E-08	1.044	0.016	1.065	2.0E-08
8,FSTL5	rs11724116	4	162,294,038	T/C	0.16	0.90	3.3E-08	0.95	0.031	0.92	2.4E-08
9,ADCY2**	rs200550695	5	7,587,236	I/D	0.82	0.91	1.2E-07	0.94	0.011	0.92	1.5E-08
10,SSBP2	rs10035291	5	80,796,368	T/C	0.68	1.081	1.1E-07	1.047	0.018	1.070	2.7E-08
11,RIMS1	rs57970360	6	72,519,394	D/I	0.44	1.066	3.1E-06	1.062	0.0016	1.064	3.5E-08
12,POU3F2**	rs2388334	6	98,591,622	A/G	0.52	0.93	8.6E-08	0.95	0.0051	0.94	4.0E-09
13,RPS6KA2	rs10455979	6	166,995,260	C/G	0.53	0.93	4.6E-08	0.97	0.046	0.94	4.3E-08
14,THSD7A	rs113779084	7	11,871,787	A/G	0.30	1.068	7.3E-06	1.095	2.9E-05	1.076	2.5E-09
15,SRPK2	rs73188321	7	105,048,158	T/C	0.33	0.92	7.0E-08	0.94	0.0015	0.92	1.1E-09
16,MRPS33	rs201231874	7	140,700,006	D/I	0.25	0.92	9.4E-08	0.93	0.0008	0.92	6.2E-10
17,ANK3**	rs10994318	10	62,125,856	C/G	0.057	1.151	4.5E-07	1.130	0.0021	1.145	6.8E-09
18,ADD3**	rs59134449	10	111,745,562	I/D	0.16	1.105	5.0E-08	1.059	0.017	1.090	1.2E-08
19,FADS2**	rs12226877	11	61,591,907	A/G	0.29	1.095	1.2E-08	1.062	0.0073	1.085	9.9E-10
20,PACS1	rs10896090	11	65,945,186	A/G	0.81	1.094	2.1E-07	1.062	0.0089	1.084	1.9E-08
21,PC	rs7122539	11	66,662,731	A/G	0.35	0.93	2.2E-07	0.96	0.015	0.94	3.8E-08
22,SHANK2	rs12575685	11	70,517,927	A/G	0.31	1.066	1.2E-05	1.088	5.7E-05	1.073	7.7E-09
23,CACNA1C**	rs10744560	12	2,387,099	T/C	0.34	1.087	2.9E-09	1.052	0.0086	1.076	3.6E-10
24,STARD9	rs4447398	15	42,904,904	A/C	0.12	1.112	1.1E-07	1.072	0.0079	1.099	9.4E-09
25,ZNF592	rs139221256	15	85,357,857	I/D	0.28	0.92	8.5E-09	0.97	0.082	0.93	2.7E-08
26,GRIN2A	rs11647445	16	9,926,966	T/G	0.65	0.93	1.2E-07	0.93	9.8E-05	0.93	1.1E-10
27,HDAC5	rs112114764	17	42,201,041	T/G	0.69	0.93	1.7E-06	0.94	0.0021	0.93	2.5E-08
28,ZCCHC2	rs11557713	18	60,243,876	A/G	0.29	1.074	1.2E-06	1.059	0.0038	1.069	3.6E-08
29,NCAN**	rs111444407	19	19,358,207	T/C	0.15	1.124	2.4E-10	1.040	0.075	1.097	1.3E-09
30,STK4	rs202012857	20	43,682,549	I/D	0.28	0.923	3.0E-07	0.942	0.0043	0.929	1.1E-08
B. Additional loci with lead SNP P < 5×10−8 in GWAS analysis
TFAP2B	rs55648125	6	50,816,718	A/G	0.90	0.89	4.9E-08	0.95	0.068	0.91	8.5E-08
DFNA5	rs17150022	7	24,771,777	T/C	0.88	0.89	2.7E-08	0.96	0.087	0.91	8.6E-08
SLC25A17	rs138321	22	41,209,304	A/G	0.50	1.083	4.7E-09	1.012	0.28	1.060	1.9E-07
HLF	rs884301	17	53,367,464	T/C	0.37	1.084	5.8E-09	1.013	0.26	1.061	2.1E-07
PHF15	rs329319	5	133,906,609	A/G	0.43	1.082	1.5E-08	1.019	0.18	1.061	2.1E-07
ODZ4**	rs73496688	11	79,156,748	A/T	0.14	1.11	1.0E-08	1.016	0.29	1.083	4.2E-07
[Intergenic]***	rs57681866	2	57,975,714	A/G	0.06	0.85	5.0E-08	0.97	0.23	0.89	1.2E-06
[Intergenic]***	rs13231398	7	110,197,412	C/G	0.11	0.89	3.4E-08	0.998	0.47	0.92	4.6E-06

Loci are numbered 1 to 30, ordered by genomic position, with previously reported gene name for published loci

P-values for GWAS and combined analyses are two-tailed, bold and underlined if p < 5×10−8.

P-values for follow-up are one-tailed based on the direction of effect in the discovery GWAS, bold and underlined if p < 0.05.

Previously published and named loci. (Locus 12 would be named as Intergenic, nearest gene is POU3F2 691Kb.)

Intergenic loci nearest genes: Locus 4 PCGEM1 824kb, Table 1B chr2 locus VRK2 298Kb, Table 1B chr7 IMMP2L 106Kb.

Follow-up of suggestive loci in additional samples

We tested lead variants that were significant at P<1×10−4 in our discovery GWAS meta-analysis, a total of 794 autosomal and 28 X chromosome variants, for association in follow-up samples totaling 9,412 cases and 137,760 controls of European ancestry (effective sample size 23,005; Supplementary Note and Supplementary Table 1B). We first compared discovery and follow-up sample summary statistics using LD score regression, and estimated their genetic correlation to be 0.98 (se=0.07), consistent with homogeneous genetic effects between the two samples. Discovery and follow-up samples also show similar patterns of significant genetic correlations with a range of other human diseases and traits in the LD Hub database [26] (Supplementary Table 3; correlation of 0.93, p = 8.3×10−14, Supplementary Figure 3).

Thirty autosomal loci achieved genome-wide significance (P< 5×10−8) in fixed-effect meta-analysis of our GWAS and follow-up samples (Figure 1, Table 1A, Supplementary Data 1–3, Supplementary Table 4). In Supplementary Table 5, we present detailed descriptions of the associated loci and genes, with bioinformatic and literature evidence for their potential roles in BD. Of the 30 genome-wide significant loci from our combined analysis, 20 are novel BD risk loci. These include 19 loci that were significant only in the combined analysis, of which three were reported to have genome-wide significant SNPs in previous studies (ADCY2 [18], POU3F2 [18], ANK3 [12,18]), and 11 that were significant in our primary GWAS. We refer to loci by the gene name attributed in previous BD GWAS publications, or by the name of the closest gene for novel loci, without implication that the named gene is causal. Results for all variants tested in the follow-up study are presented in Supplementary Table 4.

Figure 1.

Manhattan plot for our primary genomewide association analysis of 20,352 cases and 31,358 controls. GWAS -log10P-values are plotted for all SNPs across chromosomes 1-22 (diamonds, green for loci with lead SNP GWAS P < 10−6). Combined GWAS+followup −log10P-values for lead SNPs reaching genome-wide significance in either GWAS or combined analysis (triangles, inverted if GWAS+followup −log10P > GWAS -log10P). Labels correspond to gene symbols previously reported for published loci (black) and the nearest genes for novel loci (blue), at top if GWAS+followup P < 5×10−8. Loci with one-tailed follow-up p > 0.05 (Table 1) have dotted underlined locus names.

Of the 19 variants that were genome-wide significant in the discovery GWAS, 8 were not genome-wide significant in the combined analysis (Table 1B), and 11 were non-significant in one-tailed association tests in the follow-up samples (p>0.05 in Table 1). Still, the follow-up results for these 19 variants are clearly non-null in aggregate: all 19 had consistent directions of effect between discovery GWAS and follow-up (9.5 expected by chance, binomial test p=4×10−6), and eight of the 19 had follow-up 1-tailed p<0.05 (1 expected by chance, sign test p=2×10−6). Using effect sizes corrected for winner’s curse [27,28] for each of the 19 variants that were genome-wide significant in the GWAS, we calculated power to achieve significant results (1-tailed p<0.05) in the follow-up samples or genome-wide significance in combined analysis (Supplementary Note, Supplementary Table 6, Supplementary Figure 4). We found that the number of variants significant in follow-up is close to expectation (8 observed with follow-up p<0.05, 8.26 expected, Poisson binomial p = 0.57), and that 11 variants achieving genome-wide significance in the combined analysis is also within the expected range (p = 0.29). As an alternative to winner’s curse correction, we conducted a polygenic inference analysis using a mixture of Gaussian effect size distributions to model BD genetic architecture and estimate the variants’ true effect sizes [29] (Supplementary Note, Supplementary Figure 5). Under this model, we found that only two variants were nominally significantly weaker in follow-up than expected by chance (TRANK1 rs9834970 p = 0.012, and rs13821 p = 0.026; Supplementary Table 7), and none were Bonferroni significant (p>0.05/19=0.0026). Thus, the overall replication rate is within the expected range given the polygenic architecture of BD.

We next asked if the variants tested in the follow-up samples were, in aggregate, consistent with the presence of additional sub genome-wide significant BD association signals. After excluding 47 variants that were genome-wide significant in our GWAS, our combined analysis or previous BD GWAS, 775 variants remained in our follow-up experiment. 551 variants had the same direction of effect in the discovery GWAS and follow-up (71%, compared to a null expectation of 50%, sign test p = 1.3×10−32), and 110 variants had the same direction of effect and were nominally significant (p<0.05) in the follow-up (14%, compared to an expected value of 5%, binomial test p = 2.1×10−22). This consistency between our GWAS and follow-up results suggests that many more true BD associations exist among these variants.

To identify additional independent signals, we conducted conditional analyses across each of the 30 significant BD loci (Supplementary Table 8). We used the effective number of independent variants based on LD structure within loci [30] to calculate a multiple test-corrected significance threshold (p=1.01×10−5, see Supplementary Note). Only one locus showed evidence for an independent association signal (rs114534140 in locus #8, FSTL5; pconditional = 2×10−6). At another locus (#30,STK4 on chr 20), we found two SNPs with genome-wide significance in low LD (R2 < 0.1); however, conditional analysis showed that their associations were not independent.

Shared loci and genetic correlations with other traits

We next examined the genetic relationships of BD to other psychiatric disorders and traits. Of the 30 genome-wide significant BD loci, 8 also harbor schizophrenia (SCZ) associations [31,32]. Based on conditional analyses, the BD and SCZ associations appear to be independent at 3 of the 8 shared loci (NCAN, TRANK1 and chr7q22.3:105Mb) (Supplementary Table 9). No genome-wide significant BD locus overlapped with those identified for depression (DEPR), including 44 risk loci identified in the most recent PGC GWAS of major depression [33], and those reported in a large study of depressive symptoms or subjective well-being [34]. As previously reported [35], we found substantial and highly significant genetic correlations between BD and SCZ [31] (LD-score regression estimated genetic correlation rg = 0.70, se = 0.020) and between BD and DEPR [33] (rg = 0.35, se = 0.026). The BD and DEPR genetic correlation was similar to that observed for SCZ and DEPR (rg = 0.34, se = 0.025) (Supplementary Table 10A).

We found significant genetic correlations between BD and other psychiatric-relevant traits (Supplementary Table 10B), including autism spectrum disorder [8] (rg = 0.18, P=2×10−4), anorexia nervosa [36] (rg = 0.23, P=9×10−8), and subjective well-being [34] (rg = −0.22, P=4×10−7). There was suggestive positive overlap with anxiety disorders (rg=0.21, P=0.04) [37] and neuroticism (rg=0.12, P=0.002) [38]. Significant rgs were seen with measures of education: college attendance [39] (rg = 0.21, P=1=x10−7) and education years [40] (rg=0.20, P=6×10−14), but not with childhood IQ [41] (rg=0.05, P=0.5) or intelligence [42] (rg=−0.05, P=0.08). Among a large number of variants in BD risk loci that were associated with additional traits in the GWAS catalog [43], we found a handful of loci with non-independent associations (in one overlapping locus with each of educational attainment, biliary atresia, bone mineral density, lipid-related biomarkers) (Supplementary Table 9). Biliary atresia and lipid- related biomarkers, however, did not show significant genetic correlation with BD (Supplementary Table 10B).

BD subtypes

We performed a secondary GWAS focusing on three clinically recognized subtypes of bipolar disorder: BD1 (n=14,879 cases), BD2 (n=3,421 cases), and SAB (n=977 cases) (Supplementary Note, Supplementary Tables 1A & 11, Supplementary Figure 6). We observed variants in 14 loci with genome-wide significance for BD1, 10 of which were in genome-wide significant loci in the combined BD GWAS analysis. Not surprisingly given the sample overlap, 3 of the 4 remaining loci genome-wide significant for BD1 have P < 10−6 in either our discovery GWAS or combined analysis. The remaining locus (MAD1L1, chr7:1.9Mb, discovery GWAS p = 2.4×10−6) was recently published in two BD GWAS that included Asian samples [44,45]. We did not observe genome-wide significant results for the smaller BD2 and SAB analyses. BD1, BD2 and SAB all have significant common variant heritabilities (BD1 h2snp = 0.25, se = 0.014, p = 3.2×10−77; BD2 h2snp = 0.11, se = 0.028, p = 5.8×10−5; SAB h2snp = 0.25, se = 0.10, p = 0.0071). Genetic correlations among BD subtypes show that these represent closely related, yet partially distinct, phenotypes (Supplementary Table 12).

We conducted polygenic risk score (PRS) analyses to explore the relationship between genetic risk of SCZ and DEPR, and BD subtypes and psychosis (Figure 2, Supplementary Table 13). PRS calculated from SCZ [31] were significantly higher in BD1 cases than in BD2 cases (p=5.6×10−17, P threshold = 0.1) and in cases with psychosis compared to those without psychosis (p=2.12×10−6, P threshold =0.1). Conversely, PRS calculated from DEPR [33] were significantly higher in BD2 cases than in BD1 cases (P=8.5×10−10, P threshold = 0.01), independent of psychosis. Genetic correlations from LD-score regression support these results; genetic correlations were greater for SCZ with BD1 (rg = 0.71, se = 0.025) than with BD2 (rg = 0.51, se = 0.072), and were greater for DEPR with BD2 (rg = 0.69, se = 0.093) than with BD1 (rg = 0.30, se = 0.028) (Supplementary Table 12).

Figure 2.

Association of BD1 and BD2 subtypes with schizophrenia (SCZ) and major depression (DEPR) polygenic risk scores (PRS). Shown are mean PRS values (1 s.e. error bars), adjusted for study and ancestry covariates and scaled to the PRS mean and sd in control subjects, in BD1 (red) and BD2 (blue) cases, for increasing source GWAS P-value thresholds (increasing grey) as indicated. P-values (italics) test BD1 vs BD2 mean PRS, in logistic regression of case subtype on PRS with covariates. Results are detailed in Supplementary Table 13.

Systems biology and in silico functional analyses

We tested for functional genomic enrichment in our BD GWAS using partitioned LD-score regression and a range of functional annotations across tissues [46] (Supplementary Note, Supplementary Table 14). SNP-based BD heritability was most enriched in open chromatin annotations in the central nervous system (proportion SNPs = 0.14, proportion h2snp = 0.60, enrichment =3.8, p = 3 × 10−14). We also used DEPICT [47] to test for expression of BD-associated genes across tissues, and found significant enrichment of central nervous system (p < 1.4×10−3, FDR < 0.01) and neurosecretory system (p = 2.0×10−6, FDR < 0.01) genes (Supplementary Table 15).

To prioritize genes that may play a functional role in BD, we integrated BD GWAS association statistics with eQTL (SNP-gene expression association) and mQTL (SNP-DNA methylation association) data using summary Mendelian randomization (SMR) [48,49,50] (Supplementary Table 16; Supplementary Note). SMR identified 21 genes using eQTL data that were significant after multiple testing correction, without evidence of heterogeneity between GWAS and eQTL association signals. Association with GNL3 was observed in both brain and blood, highlighting the utility of using blood eQTL data as proxy for brain eQTLs [50]. Methylation profiles at 6 CpGs in brain and 10 CpGs in blood were associated with BD, four of which overlapped between brain and blood mQTL: MUSTN1, GLT8D1, HAPLN4 and FADS2.

Finally, we used MAGMA [51] to conduct a gene-wise BD GWAS and to test for enrichment of pathways curated from multiple sources (see Supplementary Note). We note that significance levels were assigned to genes by physical proximity of SNPs, and do not imply that significant genes are causal for BD. Genic association results included 154 Bonferroni significant genes (MAGMA pJOINT < 2.8×10−6), including 82 genes in 20 genome-wide significant loci, and 73 genes in 27 additional loci that did not reach genome-wide significance (Supplementary Table 17). Nine related pathways were significantly enriched for genes with BD associations (p < 7.0×10−5, FDR < 0.05), including abnormal motor coordination/balance pathways (from mice), regulation of insulin secretion, and endocannabinoid signaling (Supplementary Table 18, Supplementary Figure 7).

DISCUSSION

We carried out a large bipolar disorder (BD) GWAS and identified 30 genome-wide significant loci, including 20 that were novel. Previous BD GWAS have reported a total of 20 loci significantly associated with BD[9-23] ; twelve of these previously reported loci were not genome-wide significant in our GWAS meta analysis, but all had PGWAS ≤ 1.3×10−5 (Supplementary Table 4C). Our recent GWAS of BD and SCZ [52], which included our discovery GWAS data jointly analyzed with published SCZ data [31] (without overlapping control subjects), highlighted similarities and differences in BD and SCZ in terms of known associated SNPs and PRS-subphenotype associations; here we maximized power to identify BD associations. Phenotypic variance explained by polygenic risk scores (PRS) based on our BD GWAS data is ~8% (observed scale; 4% on the liability scale [53]), an increase from 2.8% (1.2% on the liability scale) in our previous study [9]. The results of our BD subtype PRS analyses support the nosological distinction between BD1 and BD2, but also highlight the importance of psychosis beyond DSM subtypes, corroborating and expanding evidence from previous clinical [54] and genetic studies [52,55,56]. The DEPR vs. BD PRS analyses provide further support for the distinction between BD1 and BD2, independent of the presence of psychosis.

Of the 19 loci identified in our discovery GWAS, only 11 were genome-wide significant in meta-analysis of our GWAS and follow-up samples. These results are not unexpected given small effect sizes, the winner’s curse [28,57] (Supplementary Note and Supplementary Figure 4); SNPs can teeter-totter around the genome-wide significance threshold even as sample sizes increase. Genetic heterogeneity observed among BD GWAS cohorts[8] could also contribute to inconsistent replication results; we observed variable polygenic effects between BD subtypes (Figure 2, Supplementary Table 13) as well as between cohorts in our study (Supplementary Figure 2, Supplementary Table 4) which used a diversity of criteria to define cases (Supplementary Note). Remarkably, the strongest association signal from the discovery GWAS, at the TRANK1 locus (rs9834970; pcombined = 5.7E-12, OR = 0.93), exhibited significant heterogeneity among discovery GWAS cohorts (Cochran’s Q p = 1.9×10−4), and did not replicate in the follow-up sample (1-tailed pfollowup = 0.3) (Supplementary Data 2 & 3). This locus has been significant in recent [11,12,17,18] but not earlier BD GWAS [9,13,20]. Thus, complex genetic architecture as well as phenotypic heterogeneity may contribute to the inconsistency of genome-wide significant findings within and across BD GWAS studies. The observed heterogeneity is a major challenge for GWAS of psychiatric disorders and calls for careful and systematic clinical assessment of cases and controls in parallel with continued efforts to collect larger sample sizes.

Of the 30 BD associated loci, 8 also harbor associations [31,32,58] with schizophrenia (SCZ); however, conditional analyses suggest that the BD and SCZ associations at 3 of the 8 shared loci (in the NCAN, TRANK1 and chr7q22.3:105 Mb loci) may be independent (Supplementary Table 9). Differential BD and SCZ associations may represent opportunities to understand the genetic distinctions between these closely related and sometimes clinically difficult to distinguish disorders. We did not find BD loci that overlap with those associated with major depression[33].

The confirmed association within loci containing CACNA1C and other voltage-gated calcium channel genes supports the rekindled interest in calcium channel antagonists as potential treatments for BD, with similar examination ongoing for other genes implicated in SCZ GWAS [59]. Other genes within novel BD-associated loci include those coding for other ion channels and transporters (SCN2A, SLC4A1), neurotransmitter receptors (GRIN2A) and synaptic components (RIMS1, ANK3). Further study will confirm whether or not these are the causal genes in the loci. These processes are important in neuronal hyperexcitability[60], an excess of which has been reported in iPSC derived neurons from BD patients, and which has been shown to be affected by the classic mood stabilizing drug lithium [61]. In addition, SMR eQTL and mQTL analyses implicate GLT8D1, which is involved in proliferation and differentiation of neural stem cells [62]. Pathway analyses reveal new genetic evidence for insulin secretion and endocannabinoid signaling in BD. There is evidence of insulin action in the brain [63] and in BD [64]. The endocannabinoid system has possible roles in schizophrenia [65,66] and depression [67]. Top genes appearing in these pathways include calcium and potassium channel subunit, MAP kinase and GABA-A receptor subunit genes (Supplementary Table 18).

We observe significant positive genetic correlations with educational attainment, but not with either adult or childhood IQ, suggesting that the role of BD genetics in educational attainment may be independent of general intelligence. This result is inconsistent with suggestions from epidemiological studies [68], but in agreement with a recent clinical study [69].

In summary, findings from the genome-wide analysis of BD reveal an extensive polygenic genetic architecture of the disease, implicate brain calcium channels and neurotransmitter function in BD etiology, and confirm that BD is part of a spectrum of highly correlated psychiatric and mood disorders.

Accession codes:

Accessioned data were part of the following datasets analyzed in this study. Fat2: phs000167.v1.p1 or PGC bundle phs001254.v1.p1 (MGS nonGAIN controls). Gain: dbGAP phs000017.v3.p1 (GAIN Bip cases/controls), phs000021.v3.p2 (GAIN scz controls); also part of the PGC dbGAP bundle phs001254.v1.p1. Jjst: phs000092.v1.p1 or PGC bundle phs001254.v1.p1 (SAGE controls). St2c, NIMH RGR Bipolar Study 19 (STEP-BD), dbGAP phs000294.v1.p1 (MIGEN controls). Mich: NIMH RGR Bipolar Study 2 (Pritzker). Wtcc: EGAD00000000002.

Data Availability.

The PGC’s policy is to make genome-wide summary results public. Summary statistics for our meta-analysis are available through the PGC (https://www.med.unc.edu/pgc/results-and-downloads). Data are accessible with collaborative analysis proposals through the Bipolar Disorder working group of the PGC (https://med.unc.edu/pgc).

ONLINE METHODS

GWAS and follow-up cohorts.

Our discovery GWAS sample was comprised of 32 cohorts from 14 countries in Europe, North America and Australia (Supplementary Table 1A), totaling 20,352 cases and 31,358 controls of European descent. A selected set of variants (see below) were tested in 7 follow-up cohorts of European descent (Supplementary Table 1B), totalling 9,025 cases and 142,824 controls (Neff = 23,991). The Supplementary Note summarizes the source and inclusion/exclusion criteria for cases and controls for each cohort. All cohorts in the initial PGC BD paper were included [9]. Cases were required to meet international consensus criteria (DSM-IV or ICD-10) for a lifetime diagnosis of BD established using structured diagnostic instruments from assessments by trained interviewers, clinician-administered checklists, or medical record review. In most cohorts, controls were screened for the absence of lifetime psychiatric disorders and randomly selected from the population.

GWAS cohort analysis

We tested 20 principal components for association with BD using logistic regression; seven were significantly associated with phenotype and used in GWAS association analysis (PCs 1-6, 19). In each cohort, we performed logistic regression association tests for BD with imputed marker dosages including 7 principal components to control for population stratification. For all GWAS cohorts, X-chromosome association analyses were conducted separately by sex, and then meta-analyzed across sexes. We also conducted BD1, BD2, and SAB GWAS, retaining only cohorts with at least 30 subtype cases and filtering SNPs for MAF > 0.02. Results were combined across cohorts using an inverse variance-weighted fixed effects meta-analysis [70]. We used Plink ‘clumping’ [71,72] to identify an LD-pruned set of discovery GWAS meta-analysis BD-associated variants (P < 0.0001, and distance >500 kb or LD r2 < 0.1, n variants =822) for analysis in the follow-up cohorts. Conditional analyses were conducted within each GWAS cohort and meta-analyzed as above.

Follow-up cohort analysis.

In each follow-up cohort we performed BD association analysis of the 822 selected GWAS variants (when available) including genetic ancestry covariates, following QC and analysis methods of the individual study contributors. We performed inverse variance-weighted fixed-effects meta-analyses of the association results from the follow-up cohorts, and of the discovery GWAS and follow-up analyses.

Polygenic risk score (PRS) analyses.

We tested PRS for our primary GWAS on each GWAS cohort as a target set, using a GWAS where the target cohort was left out of the meta-analysis (Supplementary Table 2). To test genetic overlaps with other psychiatric diseases, we calculated PRS for DEPR and SCZ in our GWAS cohort BD cases [73]. In pairwise case subtype or psychosis analyses (Figure 2, Supplementary Table 13), we regressed outcome on the PRS adjusting for ancestry principal components and a cohort indicator using logistic regression, and visualized covariate-adjusted PRS in BD1 and BD2 subtypes (Figure 2). Outcome sample sizes were BD1 n=8,044, BD2 n=3,365, SAB n=977; BD1 cases with and without psychosis n= 2175 and 798 respectively, BD2 cases with and without psychosis n= 146 and 660.

Linkage disequilibrium (LD) score regression.

LD score regression [25,26] was used to conduct SNP-heritability analyses from GWAS summary statistics. LD score regression bivariate genetic correlations attributable to genome-wide common variants were estimated between the full BD GWAS, BD subtype GWASs, and other traits and disorders in LD-Hub [26]. We also used LD score regression to partition heritability by genomic features [46].

Relation of BD GWA findings to tissue and cellular gene expression.

We used partitioned LD score [47,74] and DEPICT [47] regression to evaluate which somatic tissues and brain tissues were enriched in the BD GWAS. We used summary-data-based Mendelian randomization (SMR) [48,50] to identify SNPs with strong evidence of causality of brain or blood gene expression or methylation in BD risk (Supplementary Table 16), with a test for heterogeneity to exclude regions with LD between distinct causal SNPs (pHET < 0.01).

Gene-wise and pathway analysis.

Guided by rigorous method comparisons conducted by PGC members [51,75], p-values quantifying the degree of association of genes and gene sets with BD were generated using MAGMA (v1.06) [51]. We used ENSEMBL gene coordinates for 18,172 genes giving a Bonferroni corrected P-value threshold of 2.8×10−6. Joint multi-SNP LD-adjusted gene-level p-values were calculated using SNPs 35 kb upstream to 10 kb downstream, adjusting for LD using 1,000 Genomes Project (Phase 3 v5a, MAF ≥ 0.01, European-ancestry subjects) [76]. Gene sets were compiled from multiple sources. Competitive gene set tests were conducted correcting for gene size, variant density, and LD within and between genes. The pathway map (Supplementary Figure 5) was constructed using the kernel generative topographic mapping algorithm (k-GTM) as described by [77].

Genome build.

All genomic coordinates are given in NCBI Build 37/UCSC hg19.