Genetic association of cyclic AMP signaling genes with bipolar disorder.

The genetic basis for bipolar disorder (BPD) is complex with the involvement of multiple genes. As it is well established that cyclic adenosine monophosphate (cAMP) signaling regulates behavior, we tested variants in 29 genes that encode components of this signaling pathway for associations with BPD type I (BPD I) and BPD type II (BPD II). A total of 1172 individuals with BPD I, 516 individuals with BPD II and 1728 controls were analyzed. Single SNP (single-nucleotide polymorphism), haplotype and SNP × SNP interactions were examined for association with BPD. Several statistically significant single-SNP associations were observed between BPD I and variants in the PDE10A gene and between BPD II and variants in the DISC1 and GNAS genes. Haplotype analysis supported the conclusion that variation in these genes is associated with BPD. We followed-up PDE10A's association with BPD I by sequencing a 23-kb region in 30 subjects homozygous for seven minor allele risk SNPs and discovered eight additional rare variants (minor allele frequency < 1%). These single-nucleotide variants were genotyped in 999 BPD cases and 801 controls. We obtained a significant association for these variants in the combined sample using multiple methods for rare variant analysis. After using newly developed methods to account for potential bias from sequencing BPD cases only, the results remained significant. In addition, SNP × SNP interaction studies suggested that variants in several cAMP signaling pathway genes interact to increase the risk of BPD. This report is among the first to use multiple rare variant analysis methods following common tagSNPs associations with BPD.

Introduction

Bipolar disorder (BPD) affects ∼2.5% of the United States population aged 18 and older based on data from the population-based National Comorbidity Survey Replication.[1] The risk of BPD for a first-degree relative of an affected individual is 9%, and the concordance rate for monozygotic twins is 40–45%.[2] Relatives of BPD probands are also at an increased risk of other related psychiatric disorders, such as psychotic, anxiety, substance abuse and impulse control disorders.[1, 2, 3, 4] With an overall heritability estimated at 80–90%,[2] these and other findings make clear that susceptibility to BPD has a very strong genetic basis that overlaps with the genetic susceptibility for other neuropsychiatric disorders.

The two most common subtypes of BPD are bipolar disorder I (BPD I) and bipolar disorder II (BPD II). The lifetime risks in the US population is 1% for BPD I and 1.1% for BPD II.[1] The clinical course of BPD I and BPD II differ in that BPD I patients exhibit one or more manic episodes while BPD II patients are characterized by recurrent depressive episodes.[2, 3] Recurrence of manic/depressive episodes and associated disabilities is common5, 6, 7 and there exists no curative treatment. Although BPD I and BPD II have been shown to co-segregate within families,[2, 8] the diverse clinical symptoms indicate that subtypes of the BPD phenotype exist and/or that phenotypic modifiers are involved.[9] The identification of BPD I and BPD II susceptibility genes may lead to more effective targets for therapy.

Several recent genome-wide genetic association studies and meta-analyses[2, 10] have jointly analyzed BPD I and BPD II subjects,[10, 11, 12, 13, 14, 15, 16] or have focused only on BPD I subjects.[17, 18, 19] One motivation for grouping the subtypes is that it remains possible that BPD II is a milder form of BPD I as opposed to being a distinct disorder. However, genetic heterogeneity can reduce the power to detect susceptibility variants when samples of BPD I and BPD II are combined for analysis. For this study, we analyzed BPD I and BPD II separately to increase the likelihood of genetic homogeneity among cases.

A substantial body of evidence indicates a relationship between BPD and cyclic adenosine monophosphate (cAMP) signaling although the specific nature of the relationship remains obscure (reviewed in Gould and Manji[20] and Dwivedi and Pandy[21]). Several studies on peripheral cells or post-mortem brain have reported that basal and receptor-activated adenylyl cyclase (AC) activities are increased in BPD patients.[22, 23] Other studies have reported that the mood stabilizers lithium and carbamazepine interact directly or indirectly to attenuate the receptor-activation of ACs,[24, 25, 26, 27, 28, 29] consistent with the possibility that an increased AC activity or activation is a biochemical phenotype associated with BPD. Several studies have shown that mRNA levels for Gαs, the stimulatory G-protein that couples with AC, are increased in peripheral cells of BPD patients.[30, 31, 32, 33, 34] Post-mortem analyses on brain have demonstrated an increased abundance of the Gαs protein, but not Gαi or Gαo,[32, 30] and lithium has been found to decrease the abundance of mRNAs for several G-proteins.[35, 36] Lithium may also decrease the coupling of receptors to G-proteins.[22, 26, 37, 38, 39] Post-mortem analyses have demonstrated a decreased level of the regulatory subunits for protein kinase A (PKA) in BPD,[40] with increased PKA activity observed in lymphocytes and platelets from BPD subjects.[41, 42] Some substrates of PKA have been found to be hyperphosphorylated and increased PKA activity and immunoreactivity have been observed in post-mortem brain[43, 44] and in peripheral cells[41, 45, 46, 47] of BPD patients. Finally, the antidepressant rolipram inhibits cAMP phosphodiesterases (PDE) of the dunce class of PDEs (PDE4),[48] implicating this class of enzymes in the general state of depression. PDE4 expression is also altered in both BPD and schizophrenia.[49]

Overall, the evidence strongly suggests that aberrant cAMP signaling is part of the spectrum of biochemical phenotypes associated with BPD. Because of the well-established roles of cAMP in regulating behavioral processes, we hypothesize that genetic variation in the cAMP pathway genes is associated BPD etiology. To test this hypothesis, and to identify the specific mutations within the signaling system that confer susceptibility to BPD, we focused our association study on the collection of genes comprising the cAMP signaling system, as opposed to testing a variety of genes representing many different signaling systems with unrelated functions in the central nervous system. These analyses detected associations of PDE10A with BPD I and of DISC1 and GNAS with BPD II. In addition, significant associations were observed for interaction of variants in DISC1 and PDE4 for both BPD I and BPD II.

Materials and methods

Study samples

DNA samples were obtained from the Systemic Treatment Enhancement Program for BPD study (STEP-BPD). STEP-BPD was a national, longitudinal study designed to examine the effectiveness of BPD treatments in over 4000 patients with BPD.[50] Assignment of diagnosis was based on information derived from clinical assessment and psychiatric diagnostic interviews using the Affective Disorders Evaluation[50] and confirmed by the Mini International Neuropsychiatric Interview.[51] The affective disorders evaluation includes a modified version of the mood and psychosis modules of the Structured Clinical Interview from the American Psychiatric Association's, Diagnostic and Statistical Manual of Mental Disorders Fourth Edition (DSM-IV; 1994). The interviews were administered by trained clinical specialists using DSM-IV criteria to establish diagnoses of BPD I, BPD II, BPD NOS, schizoaffective manic or BP type or cyclothymic disorder. Only subjects with BPD I or II diagnoses were included in the analysis. The National Institute of Mental Health-funded genetics repository (STEP-Genetics Repository) that created anonymous immortalized cell lines from participants in the STEP-BPD study provided genomic DNA samples for this study. Genomic DNA from controls selected to be of European-American descent and free of signs of depression was provided by the National Institute of Mental Health Genetic Repository. Controls in the repository were screened with the Composite International Diagnostic Interview to reduce the chances of inadvertently introducing susceptibility alleles for neuropsychiatric disorders into the control population. All subjects involved in both studies gave written, informed consent. As the subject and control DNA samples were both pre-existing and de-identified, this study was considered exempt human research by the Scripps Institutional and Baylor College of Medicine IRBs.

Single-nucleotide polymorphism (SNP) selection and genotyping

We selected 29 candidate genes representing biochemical functions that have been implicated in cAMP signaling. Cross-talk among signaling pathways is extensive in cellular physiology, which blurs their boundaries, so we selected genes that are widely recognized as comprising the conical cAMP signaling pathway so as to limit our set (Table 1). For these genes, tagSNPs were selected within each gene and flanking 10 kb, using HAPMAP data (European-American Utah residents with ancestry from Northern and Western Europe, which were collected by the Center d'Etude du Polymorphisme Humaine (CEU)) using an r2⩾0.8, as a measure of linkage disequilibrium (LD) to capture variation within each gene and a minor allele frequency (MAF) ⩾0.05. A subset of 238 SNPs were selected to control for European substructure from a panel previously shown to distinguish different European populations.[52] One-thousand and fifty tagSNPs and 238 admixture SNPs were genotyped in 1742 cases and 1757 controls in addition to 58 duplicate samples using the Illumina Golden Gate array (Illumina, Inc., San Diego, CA, USA). Beadstudio software (Illumina) was used for allele detection and genotype calling.

Table 1

Cyclic AMP signaling pathway genes for which tagSNPs were genotyped

Gene symbol	Gene name	Location	tagSNPs pre-QC	tagSNPs post-QC	Imputed SNPs	Gene size (kb)
GNAI3	Guanine nucleotide-binding protein (G-protein), alpha inhibiting activity polypeptide 3	1p13	14	12	27	47
PDE4B	Phosphodiesterase 4B, cAMP-specific (phosphodiesterase E4 dunce homolog, Drosophila)	1p31	100	93	361	582
PRKACB	Protein kinase, cAMP-dependent, catalytic, beta	1p36.1	21	21	66	160
DISC1	Disrupted in schizophrenia 1	1q42.1	148	137	335	339
PPP1CB	Protein phosphatase 1, catalytic subunit, beta isoform	2p23	9	8	30	51
CREB1	cAMP responsive element binding protein 1	2q34	6	5	28	69
GNAI2	Guanine nucleotide-binding protein (G-protein), alpha inhibiting activity polypeptide 2	3p21	4	4	56	33
PRKAR2A	Protein kinase, cAMP-dependent, regulatory, type II, alpha	3p21.3-p21.2	3	3	14	97
PDE4D	Phosphodiesterase 4D, cAMP-specific (phosphodiesterase E3 dunce homolog, Drosophila)	5q12	226	219	560	1519
PPP2CA	Protein phosphatase 2 (formerly 2A), catalytic subunit, alpha isoform	5q31.1	2	2	10	30
PDE10A	Phosphodiesterase 10A	6q26	149	140	333	331
ADCY1	Adenylate cyclase 1 (brain)	7p13-p12	33	30	59	149
PRKAR1B	Protein kinase, cAMP-dependent, regulatory, type I, beta	7p22	1	1	NA	22
GNAI1	Guanine nucleotide-binding protein (G-protein), alpha inhibiting activity polypeptide 1	7q21	30	29	56	85
PRKAR2B	Protein kinase, cAMP-dependent, regulatory, type II, beta	7q22	20	20	80	117
ADCY8	Adenylyl cyclase 8 (brain)	8q24	116	108	281	260
PPP1CA	Protein phosphatase 1, catalytic subunit, alpha isoform	11q13	3	3	NA	23
PPP1CC	Protein phosphatase 1, catalytic subunit, gamma isoform	12q24.1-q24.2	4	3	8	23
GNAO1	Guanine nucleotide10-binding protein, alpha activating activity polypeptide O	16q13	55	49	125	166
AKAP10	Protein kinase A anchoring protein 10	17p11.1	6	6	22	72
NF1	Neurofibromin, type1	17q11.2	8	7	94	283
PPP1R1B	Protein phosphatase 1, regulatory (inhibitor) subunit 1B	17q12	2	2	2	10
PRKAR1A	Protein kinase, cAMP-dependent, regulatory, type I, alpha	17q23-q24	8	8	25	39
PPP4R1	Protein phosphatase 4, regulatory subunit 1	18p11.22	27	26	41	68
PRKACA	Protein kinase, cAMP-dependent, catalytic, alpha	19p13.1	1	1	1	26
PDE4C	Phosphodiesterase 4C, cAMP-specific (phosphodiesterase E1 dunce homolog, Drosophila)	19p13.11	8	7	12	47
PDE4A	Phosphodiesterase 4A, cAMP-specific (phosphodiesterase E2 dunce homolog, Drosophila)	19p13.2	3	3	1	49
GNAS	GNAS complex locus	20q13.3	31	31	18	71
GNAZ	Guanine nucleotide-binding protein(G-protein), alpha Z polypeptide	22q11.22	12	10	40	55

Abbreviation: cAMP, cyclic adenosine monophosphate; NA, not applicable; SNP, single-nucleotide polymorphism.

Data quality control

Data cleaning and analysis were performed using PLINK software.[53] Samples missing >2% of their genotypes were excluded because missingness can be indicative of low-quality DNA and remaining genotypes may have high error rates. SNP markers were excluded if they were not successfully genotyped in >95% of the samples, had a MAF <1% or deviated from Hardy–Weinberg equilibrium with a P-value <1 × 10−6 in controls. SNPs with deviation from Hardy–Weinberg equilibrium with a P-value of <1 × 10−3 were flagged for follow-up and were also removed before SNP marker loci were imputed. Samples were also removed if they exhibited excess homozygosity or excess heterozygosity (>4 s.d. from the median). In addition, identity-by-state statistics were calculated between samples to assess duplicates and to evaluate relatedness among individuals included in the study. Multidimensional scaling (MDS) was used to identify outliers using the population substructure SNPs. Samples were clustered with data from HAPMAP,[54] which was collected from individuals from CEU, Yoruba in Ibadan, Nigeria (YRI) and Japanese in Tokyo, Japan (JPT) and Han Chinese in Beijing, China (CBT). Samples were removed from further analysis if they deviated by <3 s.d. from the median of the MDS component 1 and MDS component 2 within the CEU cluster. After data cleaning, 1728 controls, 1172 BPD I and 516 BPD II remained available for analysis.

To evaluate potential population substructure/admixture, λ, the median of the distribution of test statistics over the median expected under the null of no association[55] was calculated before and after sex, age and the first MDS component were included in the logistic regression model to control for potential confounding.

Association testing

PLINK was used to perform logistic regression analyses under an additive model and was used to test for associations between single tagSNPs and BPD I and BPD II.[53] Haplotype analysis was performed using GenABEL[56] and Haplo Stats[57] using a sliding window of three SNPs. GenABEL was used to evaluate two-way interactions between tagSNPs in the cAMP signaling pathway. Findings were assessed for confounding with age and sex. The first MDS component calculated from the admixture SNPs for BPD I and for the analysis of BPD II sex was also included in the model to control for population substructure.

Imputation

MACH was used to impute additional SNPs that were not genotyped in genes within the cAMP pathway.[58, 59] It was not possible to impute additional SNPs for PRKAR2B and PPP1CA (Table 1) because of the limited number of genotyped SNPs. CEU phase haplotypes from HAPMAP phases 1 and 2[60] were used as templates to infer SNP markers that were not genotyped. The imputed SNPs with Rsq scores >0.8 were analyzed in PLINK using the probability of each genotype in the analysis, which were analyzed as a dosage.

Single-nucelotide variant discovery

Thirty BPD I patients homozygous for the PDE10A region 1 risk haplotype were randomly selected for single-nucleotide variant (SNV) variant discovery. An ∼23 kb region (chr6: 166063714–166040354) was amplified as three long-range PCR products using the Takara LA polymerase (Fisher Scientific, Pittsburgh, PA, USA). Nested PCR products ranging in size from 2–3 kb were amplified from the long-range PCR products using the Roche high fidelity and GC rich polymerase systems (Branchburg, NJ, USA). The PCR products were then sequenced using Sanger methodology (Genewiz, South Plainfield, NJ, USA), analyzed for variants using the Sequencher v.5 software (Gene Codes Corporation, Ann Arbor, MI, USA), and compared with a reference sequence (UCSC genome browser, Feb 2009 GRCh37/hg19 build). All coding exons of PDE10A isoforms 1 and 2 (Genbank NM001130690.1 and NM 006661.2) were also sequenced and analyzed to rule out any obvious non-synonomous amino acid or splice site changes that could be disease associated.

SNV genotyping and analysis

All variants from the reference sequence were investigated using dbSNP build 135. Any novel variants, not reported in dbSNP, were analyzed using Taqman SNP genotyping assays (Applied Biosystems, Carlsbad, CA, USA). The following numbers of randomly chosen individuals were analyzed for SNV allele frequencies: 807 BD I, 192 BD II and 801 controls. Taqman assays were performed in a 384-well plate format according to manufacturer's instructions on a 7900HT Fast real-time PCR instrument (Applied Biosystems).

Rare variant analysis was restricted to those SNVs with MAF <1%. We previously demonstrated that directly combining sequence and genotype data can lead to inflated type I errors in two-stage association studies.[61] This experimental design involves sequencing a portion of the cases to discover variants followed by genotyping the identified variants in the remaining sample. We therefore developed a likelihood-based method named SEQCHIP, for integrating sequencing and genotype data from two-stage association studies. It corrects for the genotypes of sequenced BPD I and BPD II cases, so that the corrected sequenced genotypes follow the same distribution as that of the genotyped samples.[62] The corrected genotypes can then be analyzed by existing rare variant tests, which were shown to have correct type I errors.[62] More specifically, the association analyses described here utilized the ANRV (aggregated number of rare variants),[63] the WSS (weighted sum statistics)[64] and the VT (variable threshold)[65] tests.

Results

Nine hundred and eighty-eight tagSNPs that met quality control criteria within 29 candidate genes (Table 1) were tested for associations using 1172 BPD I cases, 516 BPD II cases and 1728 controls. There was evidence of population substructure in the sample comparing the controls to the BPD cases (λ=1.11). The addition of one MDS component calculated from the admixture SNPs decreased λ to 1.03.

None of the association tests of the tagSNPs with either BPD I or BPD II reached genome-wide significance (P<5 × 10−8). Association results for tagSNPs with P<0.01 with BPD I and BPD II are shown in Table 2. Those SNPs that displayed a nominal significance level of P<0.05 are shown in Supplementary Table 1. An additional 2675 tagSNPs were imputed in the cAMP signaling pathway genes using the 988 tagSNPs that passed quality control (Table 1). All imputed SNPs were analyzed and results of imputed SNPs with P<0.05 are included in Supplementary Table 2.

Table 2

TagSNP associations with BPD I and BPD II at P<0.01

Chromosome	SNP	Position (GRCh37/hg19)	MAF	OR	95% CI	P	Gene
BPD I
1	rs4384209	66279125	0.23	1.2	(1.1–1.4)	1.2 × 10−3	PDE4B
1	rs600674	84701251	0.20	0.8	(0.7–0.9)	6.5 × 10−3	PRKACB
1	rs10495310	231918776	0.11	1.3	(1.1–1.5)	1.8 × 10−3	DISC1
3	rs2282751	50291785	0.14	1.2	(1.1–1.4)	6.1 × 10−3	GNAI2
5	rs17741863	58858302	0.10	1.3	(1.1–1.5)	4.0 × 10−3	PDE4D
6	rs9459417	165837992	0.13	0.8	(0.7–0.9)	8.4 × 10−3	PDE10A
6	rs470922	165844311	0.31	0.8	(0.7–0.9)	3.1 × 10−3	PDE10A
6	rs4709082	166040621	0.39	1.2	(1.1–1.3)	8.3 × 10−3	PDE10A
6	rs1033701	166040859	0.28	1.2	(1.1–1.4)	5.4 × 10−4	PDE10A
6	rs9459464	166042705	0.42	1.2	(1.1–1.3)	1.7 × 10−3	PDE10A
6	rs1079418	166047034	0.31	1.2	(1.1–1.3)	2.7 × 10−3	PDE10A
6	rs9459465	166048868	0.23	0.8	(0.7–0.9)	2.7 × 10−3	PDE10A
6	rs2983517	166051731	0.27	1.2	(1.1–1.4)	1.0 × 10−3	PDE10A
6	rs2983521	166057203	0.24	1.3	(1.1–1.5)	9.3 × 10−5	PDE10A
6	rs9365899	166062379	0.31	0.8	(0.7–0.9)	2.3 × 10−3	PDE10A
6	rs3008049	166063617	0.32	1.2	(1.1–1.4)	1.1 × 10−3	PDE10A
6	rs10214709	166066893	0.24	0.8	(0.7–1.0)	8.2 × 10−3	PDE10A
8	rs13278912	132002770	0.16	1.2	(1.1–1.4)	8.1 × 10−3	ADCY8
20	rs6064714	57414140	0.15	0.8	(0.7–0.9)	1.5 × 10−3	GNAS
20	rs6026565	57439308	0.10	1.3	(1.1–1.5)	1.7 × 10−3	GNAS
BPD II
1	rs2812391	231911976	0.27	1.2	(1.1–1.5)	5.3 × 10−3	DISC1
1	rs10495310	231918776	0.11	1.4	(1.1–1.7)	3.2 × 10−3	DISC1
1	rs4658954	231979357	0.37	0.8	(0.7–0.9)	4.0 × 10−3	DISC1
20	rs6064714	57414140	0.15	0.7	(0.6–0.9)	1.1 × 10−3	GNAS
20	rs35113254	57435532	0.23	1.4	(1.2–1.6)	5.8 × 10−5	GNAS
20	rs6026565	57439308	0.10	1.4	(1.1–1.7)	6.9 × 10−3	GNAS

Abbreviations: BPD I, bipolar disorder I; BPD II, bipolar disorder II; CI, confidence interval; MAF, minor allele frequency; OR, odds ratio; SNP, single-nucleotide polymorphism.

Twelve different tagSNPs in the PDE10A gene were associated with BPD I at P<0.01, including rs2983521 with an odds ratio (OR) of 1.3, a 95% confidence interval ranging from 1.1 to 1.5, and a P=9.3 × 10−5 (Table 2). Seven of the PDE10A tagSNPs, such as rs1033701 (OR=1.2, P=5.4 × 10−4), that associated with BPD I with ORs ranging from 1.2 to 1.3 (P<0.01) were found to be in moderate LD (r2>0.3). An additional 333 SNPs were imputed in PDE10A (Figure 1, panel a). Our analyses of these imputed SNPs indicated that an additional seven SNPs were associated with BPD I at P<0.01 (Supplementary Table 2). One example is rs2983516, with OR of 1.3 and P=5.4 × 10−5.

Figure 1

PDE10A association with bipolar disorder I (BPD I). The points in panel (a) represent genotyped (black circles) and imputed (open triangles) single-nucleotide polymorphisms (SNPs) and their probability of association with BPD I. Panel (b) illustrates probability as a sequential sliding window with three SNPs per window for the genotyped SNPs (empirical haplotype) and all SNPs (haplotype). The results were adjusted for age and population substructure. Panel (c) shows the structure of the PDE10A gene aligned with panels A and B with exons represented by bars. The direction of transcription, 5′→3′ is from right to left in this figure.

The analysis of PDE10A haplotype windows revealed a series of haplotypes that were associated with BPD I (Figure 1, panel b). Two regions in PDE10A (region 1 at 166.1 Mb and region 2 at 165.8 Mb) were associated with BPD I. Region 1 is large non-coding region within intron 1 that may be important for the regulation of PDE10A expression. Because of this possibility and the large number of associated tagSNPs that fell within this region, SNV identification was focused to a 23-kb region on 30 BPD I patients homozygous for the region 1, minor allele risk SNPs (Figure 2).

Figure 2

PDE10A region 1 risk linkage disequilibrium (LD) plot. Shading and values in cells correspond to r2 values between single-nucleotide polymorphism (SNP) pairs. The SNPs within region 1 and with odds ratio (OR) >1.0 were used to define the boundaries for subsequent sequencing.

Thirteen SNVs were identified among these BPD I cases. These SNVs were genotyped in 999 cases (807 BPD I and 192 BPD II) in addition to 801 controls. Subsequent analysis was restricted to the eight novel SNVs with a MAF <1% to exclude common polymorphisms (Table 3). Interestingly, the SNVs chr6:166054721 insA and chr6: 166046780_166046778 delATT were found simultaneously in 2 of the 807 BD I cases genotyped. These SNVs are in complete LD in all individuals genotyped. As the newly identified variants are individually rare, they were aggregated and jointly analyzed to determine their collective effect. When only genotyped samples were analyzed, the P-values were P(ANRV)=3.9 × 10−2, P(WSS)=7.1 × 10−2 and P(VT)=8.0 × 10−2. When SEQCHIP was used to integrate both sequence and genotype data, the P-values were P(ANRV)=2.2 × 10−2, P(WSS)=2.0 × 10−3 and P(VT)=3.0 × 10−2. Using SEQCHIP, genotype and sequence data can be more efficiently utilized and a stronger association signal was identified.

Table 3

PDE10A Region 1 SNVs with MAF <0.01

Observed SNP	Chromosome	Position	A1	A2	Controls (N=801)	BPD cases (N=999)	BPD I cases (N=807)
					MAF	MAF	MAF
C>T	6	166042004	T	C	0	0.000502	0.0006219
G>A	6	166045644	A	G	0.000624	0.001505	0.001863
T>G	6	166046710	G	T	0.000625	0.001503	0.001861
delATT	6	166046780_166046778	A	T	0	0.001001	0.001239
C>G	6	166052752	G	C	0.000624	0.001001	0.001239
insA	6	166054721	A	T	0	0.001001	0.001239
A>T	6	166056798	T	A	0.000624	0.003504	0.003717
C>T	6	166058713	T	C	0	0.000502	0.0006203

Abbreviations: BPD, bipolar disorder; MAF, minor allele frequency; SNP, single-nucleotide polymorphism; SNV, single-nucleotide variant.

Position: February 2009 (GRCh37/hg19) build physical position. A1: minor allele. A2: major allele.

One tagSNP in DISC1, rs10495310, was associated with BPD I (Table 2). However, three tagSNPs in DISC1 were associated with BPD II with P<0.01 (Table 2). Two of the DISC1 tagSNPs, rs2812391 and rs10495310, conferred increased risk of BPD II with ORs of 1.2 and 1.4 for each SNP, whereas the third SNP, rs4658954, conferred decreased risk with an OR of 0.8. The tagSNPs rs2812391 and rs10495310 are in mild LD (r2=0.28), whereas both are not in LD with rs4658954 (r2=0.01 and r2=0.02, respectively). Furthermore, an additional 335 SNPs were imputed in DISC1 (Figure 2, panel a) and one SNP, rs3082, was associated with BPD II with an OR of 1.3 (P=5.5 × 10−3). DISC1 haplotype analysis (Figure 2, panel b) indicated that two regions (region 1 at 231.9 Mb and region 2 at 232.1 Mb) were associated with BPD II. The DISC1 region 1 haplotype consisting of rs10495310, rs4658945 and rs16854967 had an empirical global P=1.7 × 10−3 for its association with BPD II. The DISC1 region 2 haplotype consisting of rs12137417, rs9431736 and rs7548491 had an empirical global P=3.9 × 10−2.

Two tagSNPs in GNAS, rs6064714 and rs6026565, were associated with BPD I at P<0.01 (Table 2). However, GNAS appeared to have a stronger association with BPD II in this study with three tagSNPs associated at P<0.01 (Table 2). Two of tagSNPs conferred increased risk (rs35113254 and rs6026565) for BPD II, whereas the third (rs6064714) conferred decreased risk. The SNPs rs35113254 and rs6026565 are in moderate LD (r2=0.35) but are both not in LD with rs6064714 (r2=0.04 and r2=0.01, respectively). The rs35113254 SNP in GNAS was the most significantly associated tagSNP representing the cAMP signaling pathway with P=5.8 × 10−5 with an OR of 1.4 for BPD II. Similarly, an additional 18 SNPs were imputed in GNAS (Figure 3, panel a) and an additional SNP, rs6026560, was associated with BPD II with an OR of 0.7 (P=1.4 × 10−4). Haplotype analysis in GNAS indicated a region around 5.74 Mb that was associated with BPD II (Figure 4, panel b). This region overlaps with the results from imputed and genotyped SNPs.

Figure 3

DISC1 association with bipolar disorder II (BPD II). The points in panel (a) represent genotyped. (black circles) and imputed (open triangles) single-nucleotide polymorphisms (SNPs) and their probability of association with BPD II. Panel (b) illustrates probability as a sequential sliding window with three SNPs per window for the genotyped SNPs (empirical haplotype) and all SNPs (haplotype). The results were adjusted for age and population substructure. Panel (c) shows the structure of the DISC1 gene aligned with panels (a) and (b) with exons represented by bars. The direction of transcription, 5′→3′ is from left to right in this figure.

Figure 4

GNAS association with bipolar disorder II (BPD II). The points in panel (a) represent genotyped. (black circles) and imputed (open triangles) single-nucleotide polymorphisms (SNPs) and their probability of association with BPD II. Panel (b) illustrates probability as a sequential sliding window with three SNPs per window for the genotyped SNPs (empirical haplotype) and all SNPs (haplotype). Panel (c) shows the structure of the GNAS gene aligned with panels (a) and (b) with exons represented by bars. The direction of transcription, 5′→3′ is from right to left in this figure.

Single tagSNPs in PDE4B (rs4384209), PRKACB (rs600674), GNAI2 (rs2282751), PDE4D (rs17741863) and ADCY8 (rs13278912) were also associated with BPD I with P <0.01. Additionally, imputed SNPs in PRKACB (four imputed SNPs), PDE4D (nine imputed SNPs) and ADCY8 (one imputed SNP) were associated with BPD I at P<0.01 (Supplementary Table 2). Numerous haplotype windows had P<0.05, however, only PDE10A, DISC1 and GNAS had notable peaks like those presented in Figures 1, 2, 3, 4 panel (b) for these genes.

Results for all two-way interactions between SNPs in the cAMP signaling pathway with P<0.0001 are presented in Supplementary Table 4. A total of 29 two-way tagSNP interactions with BPD I and 34 two-way tagSNP interactions with BPD II had P<0.0001. For both BPD I and BPD II, several interactions between PDE4 and DISC1 SNPs were significant. PDE4D tagSNP rs4699931 and ADCY8 tagSNP rs928136 interaction was associated with an increased ORinteraction of 3.65 for BPD I (P=1.06 × 10−5), while PDE4B tagSNP rs12404118 and DISC1 tagSNP rs823162 interaction was associated with an increased ORinteraction of 3.53 for BPD II (P=5.47 × 10−5).

Discussion

This study represents the first genetic association study of BPD I and BPD II using dozens of genes comprising the cAMP signaling pathway. Collectively, our results indicate that some of the analyzed SNPs in several genes of the pathway are in LD with causal variants for BPD I and BPD II. Furthermore, the SNV discovery aspect of our study has identified new variants in the PDE10A gene of BPD subjects that may be causally related to the disorder. Our genetic studies are consistent with prior biochemical studies implicating cAMP signaling in BPD, and highlight specific genes and SNVs that could be of major importance for these disorders.

We employed three noteworthy strategies for elucidating the genetic basis for susceptibility to BPD. First, our association study utilized 29 candidate genes representing a single-signaling pathway. The hypothesis underlying this strategy was that if one gene in the signaling pathway is involved, then casual variants in other genes of the same pathway are likely involved as well. We recognize that there exist known and obvious intersections and cross-talk between the cAMP signaling pathway and other signaling pathways, but we selected genes that are widely recognized as comprising the canonical cAMP signaling pathway in order to sample a defined and limited set of genes. A ‘molecular pathway hypothesis' as contrasted with a ‘single gene hypothesis,' may account for the genetic heterogeneity of neuropsychiatric disorders in general and the marked difficulty in discovering genetic variants that have prominent roles in these disorders. It can also explain heterogeneity in clinical phenotypes or severity within a disorder, since a variant in an AC gene, for instance, might have largely overlapping but some distinct phenotypic effects compared with a variant in cAMP PDE. The interaction analyses were motivated by the ‘molecular pathway hypothesis' as functional variants in two molecular functions of a pathway could produce a more pronounced effect than either single functional variant. For instance, the relatively large OR produced by considering an interaction between PDE4D and ADCY8 for BPD I (Supplementary Table 4) could be due to a functional variant that decreases PDE4D activity synergistically interacting with one that increases ADCY8 activity. Similarly, an interaction between variants in PDE4B and DISC1 for BPD II could be because of two separate functional variants that have synergistic effects on signaling. These ideas remain speculative and depend on the sign and magnitude that the functional variants have on the molecular factors involved.

Second, the choice of candidate molecular pathways considered the known and major role of the signaling pathway in the regulation of behavior in model systems. Genetic studies of learning and memory, circadian rhythms and other behaviors in model systems like Drosophila and the mouse have identified hundreds of genes that have behavioral roles. A good example of this is the cAMP PDE encoded by the Drosophila dunce gene, which is involved in learning and memory in Drosophila,[66] but whose human counterparts are involved in mood regulation[48] and schizoaffective disorders.[67] Candidate genes selected from the human homologs of genes with known behavioral roles in model systems should provide an extremely effective filter in gene discovery for neuropsychiatric disorders.

Third, we followed-up the initial association study with SNV discovery in the vicinity of the associated PDE10A tagSNPs that identified additional variants that may be causally involved in BPD. This required the development of a new method for integrating information from sequencing projects in a cost effective way. Using this method, named SEQCHIP, we were able to combine the sequencing data from 30 BPD I patients used for discovery with the data from the follow-up genotyping in a larger subset while accounting for potential bias. This approach was used to identify potentially causal variants that are tagged by the genotyped SNPs. When the data are integrated by SEQCHIP, all rare variant tests have correct type I errors. The genetic information from both sequence and genotype data were efficiently utilized, and therefore, more significant association signals were identified.

Suggestive evidence from other studies is consistent with the possibility that these newly discovered SNVs are involved in BPD. A recent RNA sequencing study employing human prefrontal cortex suggested that two of these SNVs (chr6:166052752 C>G and chr6: 166046780_166046778 delATT) may be associated with novel PDE10A transcripts.[68] Another study performing chromatin immunoprecipitation sequencing on neuronal cells from human prefrontal cortex tissue showed that SNVs, chr6:166058713 C>T, chr6:166056798 A>T and chr6: 166046780_166046778 delATT are associated with an enriched H3K4me3 signal, suggesting that they participate in the transcriptional regulation of PDE10A.[69] Although this preliminary evidence must be confirmed in human tissue expressing PDE10A, it does suggest that these variants may have a role in PDE10A protein expression.

Our study highlights the association of PDE10A with BPD I. Of related and high interest is the observation that the PDE10A gene resides at 6q26, within one of the four leading linkage regions for BPD.[2] The PDEs degrade cAMP and/or cyclic guanosine monophosphate to their 5′-AMP and 5′-GMP counterparts and members of this large enzyme family are known to have broad roles in controlling behavioral output in diverse organisms from flies to humans.[66, 70] The family member PDE10A codes for an 89-kDA protein capable of hydrolyzing both cAMP and cyclic guanosine monophosphate that is expressed in the brain and preferentially in the basal ganglia.[71, 72, 73, 74] This expression pattern paired with basal ganglia dysfunction in schizophrenia have led to robust efforts in identifying drugs that target PDE10A for the treatment of schizophrenia.[75, 76, 77] The association of PDE10A with BPD I expands the potential range over which PDE10A drugs might exhibit efficacy.

The gene disrupted in schizophrenia 1 (DISC1) was first implicated in psychiatric disorders after a t(1;11) translocation was discovered in a large Scottish family cosegregating with schizophrenia, BPD and recurrent major depression.[78, 79] The DISC1 region has subsequently been implicated with bipolar spectrum traits in several candidate gene association studies.[80, 81, 82, 83, 84, 85] The selection of tagSNPs has varied from study to study so we were unable to compare the precise haplotype windows between this and other studies. In addition, the majority of studies failed to control for population admixture and in most cases termed a P-value of less than nominal significance (P-value <0.05) as significant. Relative to individual tagSNP associations that have been reported, Perlis et al.[86] observed an association between rs1934909 and BPD. We observed an association between this tagSNP and BPD II (OR=0.8, P=0.01) but not BPD I cases (Supplementary Table 1). Similarly, Schosser et al.[83] observed an increased prevalence of BPD for carriers of the minor allele of rs2492367. We observed an association between this SNP and BPD II but not BPD I (OR=0.8, P=0.03).

It seems highly likely given the complexity of the DISC1 protein and it's multiple molecular roles[87] that there exist multiple functional variants linked to the various tagSNPs used in this and other studies. DISC1 interacts with several key proteins in pathways that have important roles in neuronal migration, dendritic organization and neurogenesis in addition to cAMP signaling.[67, 87, 88, 89, 90, 91, 92] The binding sites for interacting proteins are plausible regions for multiple susceptibility regions as are regions of the protein that modify binding. This could underlie the DISC1 and PDE4 tagSNP interactions observed to be significantly associated with BPD I and BPD II in this study. Furthermore, DISC1 variants in intronic regions may influence DISC1 mRNA and subsequent protein levels that also could interfere with key pathways. Indeed, Maeda et al.[93] demonstrated reduced DISC1 mRNA levels in lymphoblasts from BPD subjects compared with unaffected family controls. Thus, it is plausible that multiple genetic variants in LD with DISC1 tagSNPs exist that alter protein function or mRNA levels that influence susceptibility to neuropsychiatric disorders. Identifying these variants and characterizing their effect on DISC1 molecular functions may provide useful information for the development of novel therapeutics for bipolar spectrum traits.

The guanine nucleotide-binding protein, α-stimulating polypeptide 1 (GNAS) encodes the Gαs subunit of heterotrimeric G-proteins with Gαs activating AC for increased synthesis of cAMP. Although the exact mode of action of lithium in treating BPD remains elusive, there exists evidence that lithium decreases G-protein subunit mRNA levels.[35] Moreover, it has been speculated that lithium functions to stabilize the inactive conformation of G-proteins.[94] GNAS variants, like rs6064714 and rs6026565, may be tagging functional variants in BPD I patients that interfere with Gαs' role in regulating cAMP signaling.

It is of interest that SNPs representing certain genes, such as PDE4B and PDE10A, were found as associated with BPD I but not BPD II. It is possible that this represents an authentic genetic difference between the two BPD subtypes. However, it could also be that the power to detect the association was stronger with the BPD I population, given that the BPD I population was twice the size of the BPD II population. If the data reflect authentic genetic differences between the two subtypes, they could be important to the pharmacogenetics of BPD treatment. For instance, perhaps an altered regulation of PDE10A, specific to BPD I, is most related to full-blown mania associated with BPD I. Drugs with PDE10A action may then have a primary effect on mania.

Overall, this study provides insight into the variation that exists in genes in the cAMP signaling pathway and the association of that variation with BPD I and BPD II. Key regions were identified in PDE10A, GNAS and DISC1 that are associated with BPD I and BPD II. The effect size for these associations was modest, as observed in most other genetic association studies for neuropsychiatric disorders. Our study also illustrates the feasibility of performing large candidate gene or GWAS (genome-wide associated studies) studies followed by targeted direct sequencing to uncover novel, rare SNVs that may be associated with disease.