Assessment of Whole-Exome Sequence Data in Attempted Suicide within a Bipolar Disorder Cohort

Despite continuing work to improve the diagnosis and care for patients suffering from severe psychiatric disorders with substantially increased risk for suicidal behavior, such as bipolar disorder (BP) [1], rates of attempted and completed suicide have not fallen. Death by suicide now accounts for over 800,000 deaths each year around the world and is the second leading cause of death for individuals aged 15-29 [2]. In addition, current worldwide estimates suggest that suicide accounts for a total of 50% of all violent deaths in men and 71% of violent deaths in women, with an estimated 20 attempts for every death by suicide, imposing a terrible and ongoing societal cost [2].

Evidence from twin-based studies has demonstrated that suicidal behavior has a genetic component, with an estimated heritability of 30-50% [3]. Two sources are suspected to drive this heritability: psychiatric disorders such as mood and alcohol/substance use disorders and independent heritable factors such as impulsive aggression [4,5]. Discovering the underlying genetic basis of these factors could offer critical insight into biological pathways and mechanisms that contribute to suicidal behavior risk and provide new targets for patient assessment and treatment.

Toward this end, numerous candidate gene and pathway-driven studies have been undertaken that have identified many genetic associations with suicidal behavior (online suppl. Table S1; for all online suppl. material, see www.karger.com/doi/10.1159/000454773). The diversity and number of these findings underscore the complexity of suicidal behavior genetic risk and make follow-up efforts in the framework of individual candidate gene studies inefficient. This has led to the use of broad hypothesis-free investigation methods in many recent studies, including gene expression [6,7,8,9,10,11,12], linkage [13,14,15,16], and genome-wide association studies (GWAS) of common variants [17,18,19,20,21,22,23,24]. These studies have identified additional targets that may contribute to the risk of suicidal behavior, but very few implicated sites have been replicated [25,26]. Further, few implicated candidate genes have been deeply assessed for the possibility that rare functional variation within them might contribute to the phenotype [27,28,29].

The primary aim of our study was to broadly examine rare functional variation throughout the human exome in order to identify individual variants, genes, and/or pathways with significant variation differences between individuals who have attempted suicide and those who have not. To do this, we took advantage of a large BP whole-exome sequencing project that has generated data on the coding exons of approximately 19,000 genes in 1,018 BP subjects with an available suicide attempt history. The resulting attempted-suicide whole-exome sequencing project is the first such effort in the field of suicide genetics.

Our BP sample consists of 1,018 age-matched (p = 0.45) unrelated individuals of European-American ancestry. These subjects overlapped with our prior GWAS study (91% overlap) [19]. Briefly, subjects were diagnosed as BP type 1 (942 subjects), schizoaffective disorder bipolar type (75 subjects), or BP not otherwise specified (1 subject) in accordance with research diagnostic criteria, Diagnostic Statistical Manual (DSM)-III-R, or DSM-IV criteria. All individuals were interviewed using either the Diagnostic Interview for Genetic Studies (DIGS) [30] or the Schedule for Affective Disorders and Schizophrenia (SADS) [31], which both include self-reported suicide attempt histories. Subjects were also asked about previous suicide attempts and their intent to die. Subjects were included in the study that either had no self-reported suicide attempts (nonattempters; 631 subjects) or had at least 1 self-reported suicide attempt with definite or serious intent to die (attempters; 387 subjects). All included subjects, following a complete explanation of the parent study, supplied institutional review board-approved written consent. Additional sample set details, including demographic details (online suppl. Table S2), can be found in the online supplementary materials and methods.

Sequencing was performed according to the NimbleGen SeqCap EZ Exome protocol (NimbleGen, Madison, WI, USA). Target capture used NimbleGen SeqCap EZ arrays v1, v2, or v2 plus added targets for the promoter and untranslated regions of 1,422 neuronal postsynaptic density genes [32] and 57 genes suspected to be of importance in BP that were obtained from a number of individual candidate gene studies, including those cited here [33,34,35,36] (online suppl. Table S3). The additional targets in the v2+ array represent potential regulatory sites that have been shown to be particularly enriched for associations with complex disease in past studies [37]. The inclusion of these regulatory sites within genes of suspected importance in psychiatric disease allows us to economically increase the potential yield of the data set. Paired-end sequencing was performed on an Illumina GA-IIx or HiSeq 2000 (Illumina, San Diego, CA, USA).

Individual sample data were processed via a Burrows Wheeler Aligner [38], SAMtools [39], BAMtools [40], Picard (http://picard.sourceforge.net), and the Genome Analysis Tool Kit (GATK) [41] pipeline. GATK UnifiedGenotyper processing of Reducereads [42] files allowed final genotyping of all subjects at once. Variants were predicted to be “functional” based on annotation data and were classified into 2 levels of evidence for both coding and noncoding (“regulatory”) sites. Coding variants were annotated by ANNOVAR [43]. Stopgain and essential splice site variants were classified as “coding disruptive.” All nonsynonymous variants predicted to be damaging by at least 1 of 6 bioinformatic packages were combined with all coding disruptive variants under the classification “coding broad.” Both coding classifications were modeled after a recent schizophrenia exome study [44]. For consistency, regulatory variants were annotated by RegulomeDB [45] with scores of 1-2 being classified as “regulatory narrow” and 1-6 as “regulatory broad.”

The final variant set included only calls/sites that passed GATK variant recalibration, Hardy-Weinberg (p > 1 × 10^-6), depth ≥10, and genotyping quality ≥20. Detected insertion/deletions (indels) and tri/quad allelic sites were removed from the data set due to the technological limitations of accurately calling these complex alleles. Subjects were assessed via principal component analyses (PCA) and a sex-check algorithm to remove mislabeled or outlier individuals. PCA components showed no clustering based on phenotype or platform/assay (online suppl. Fig. S1). In addition, subjects were required to have depth ≥20× in ≥70% of the targeted sites. Finally, singleton variant distribution was assessed in attempters and nonattempters, demonstrating no systematic genotyping bias between groups (p = 0.29). Additional details regarding the sample prep, platforms used, assay versions, genotyping pipeline, and quality control measures can be found within the online supplementary materials and methods.

All single-variant tests and collapsed-variant tests (“gene burden”; using the Combined Multivariate Collapsing method [46]) were performed using the R-package “logistf” [47]. Additional gene tests were performed via the Sequence Kernel Association Test (SKAT) [48] for autosomal genes only. Pathway analyses were divided into 2 categories. In a primary pathway analysis, 33 pathways hypothesized to be potentially important within suicidal behavior risk were assessed for genetic association via PLINK/SEQ v0.10 (https://atgu.mgh.harvard.edu/plinkseq/) SMP tests (online suppl. Table S4), which allows for correction of potential systematic variant calling biases between comparison groups. Additionally, a secondary pathway analysis was performed within a more comprehensive set of 3,621 pathways derived primarily from the molecular signature database (MSigDB) [49] which includes Gene Ontology (GO) [50], Kyoto Encyclopedia of Genes and Genomes (KEGG) [51], and many other annotated functional pathways. Gene and pathway association tests examined coding disruptive/broad and regulatory narrow/broad sets independently.

All single-variant and gene-level tests were corrected for sex, sequencing platform/array, and the first 5 PCA components for each subject. PLINK/SEQ pathway SMP tests were corrected for sequencing platform/array version. The top pathway association results were also separately controlled with alcohol dependence as a covariate due to an enrichment of genes that are known to be involved in alcohol metabolism in several of the top pathways.

All gene, pathway, and the overall functional variant enrichment tests were performed using rare variants within 2 minor allele frequency (MAF) thresholds: ≤0.05 MAF (MAF05) and ≤0.01 MAF (MAF01). Variant thresholds were defined by the allele frequency within our complete data set in addition to European 1,000 genomes [52], Non-Finnish European subjects from the Exome Aggregation Consortium (ExAC 0.3; http://exac.broadinstitute.org), and the National Heart Lung and Blood Institute (NHLBI) GO Exome Sequencing Project (http://evs.gs.washington.edu/EVS/) data, if available. To control for multiple testing, all test results were compared against conservative Bonferroni thresholds: p < 5 × 10^-8 (∼1,000,000 tests) for single variant, conservative p < 1.0 × 10^-6 (49,452 tests) and a more liberal p < 2.5 × 10^-6 (∼20,000 genes assessed) for gene level, and p < 1.4 × 10^-5 (3,654 tests) for pathway association tests. In addition, PLINK/SEQ SMP pathway association analyses were further assessed via permutation by swapping attempter/nonattempter labels 500 times and re-running the analyses to develop an empirical corrected p value for each pathway result. Levels of suggestive significance for all analyses were defined as any p value within approximately 1 magnitude of the established threshold for significance for that test, as defined above. Additional statistical analysis details may be found within the online supplementary materials and methods.

This study utilized high-quality exome data for 1,018 BP subjects over a targeted region of approximately 26, 36, or 54 Mb depending on the array version used. A mean depth of coverage of 72.6× across captured targeted sites was achieved within this data set. In addition, an average of 93.9% of the target reads per subject reached at least 10× depth, the quality threshold we imposed for inclusion in our analyses (see online suppl. Table S5 for additional sample set metric details). A total of 494,475 variants were detected with 59.9% (295,959 variants) being found within coding sites and 40.1% (198,516 variants) residing in regulatory regions. Additionally, 20% of all detected variants (99,303 variants) represent novel variation not found within existing databases. Finally, the expected Ti/Tv values for exome arrays are 3.0-3.5 and 2.07-2.10 for whole-genome sequencing [41]. Our data represent expected values for exome sequencing (Ti/Tv range of 2.96-3.08 for the SeqCap EZ v1 and v2) and exome plus some selected noncoding targets in the SeqCap EZ v2+ array (Ti/Tv = 2.83; online suppl. Table S5).

No individual variant in our data was associated with attempted suicide at study-wide significance (see online suppl. Table S6 and Fig. 1). The top individual variant result was identified as rs2215955, a common nonsynonymous variant within the amphiphysin (AMPH) gene with odds ratio (OR) = 0.61 for the minor allele and nominal p = 2.8 × 10^-5. This result fell well short of the Bonferroni threshold of 5.0 × 10^-8 required for a significant result. A QQ-plot was then generated to assess the overall pattern of our individual variant tests at varying MAFs. This plot demonstrated our results followed the expected null distribution (online suppl. Fig. S2). As a result, we focused on assessing sets of variants in gene and pathway analyses. We examined variants that were bioinformatically predicted to be coding disruptive (n = 6,316), regulatory narrow (n = 4,830), coding broad (n = 129,735), or regulatory broad (n = 44,846), separately comparing each set of variants between suicide attempters and nonattempters.

No association test with any individual gene achieved study-wide significance following correction for multiple testing. Many genes, however, achieved nominal significance (p < 0.05) within the data. A total of 1,613 genes reached nominal significance with 909 genes identified in the burden analyses and 1,011 genes identified in the SKAT analyses (overlap of 307 genes reaching nominal significance in both test types). Online supplementary Table S7 outlines those genes that achieved a nominal p value <0.01, referred to as our “top genes.”

Two genes reached a level of suggestive significance: CFAP70 and SLC6A13 (Table 1). Functionally, CFAP70 [53] is suspected to play a role in cilia function which, in turn, is essential for many key developmental and cellular maintenance processes. SLC6A13 is known to encode a γ-aminobutyric acid (GABA) reuptake transporter that plays a role in regulating GABA neurotransmission [54].

We also closely examined all gene-based test signals within genes previously associated with suicidal behavior (see online suppl. Table S1). This collection of genes was assembled, in part, using the list published by Perlis et al. [17] in 2010 along with any more recent genes identified via a PubMed search with the keywords “suicide,” “genetic,” and “association.” Two previously identified genes, BDNF and DISC1, yielded nominal evidence for association (p < 0.05) of OR = 1.8 with p = 4.5 × 10^-3 for BDNF and OR = 0.72 with p = 3.5 × 10^-2 for DISC1 with both signals arising in rare regulatory broad variants.

The lack of any single gene that could survive correction for multiple testing encouraged us to examine groups of genes that operate together within biological pathways. Pathway tests afford greater detection power across biologically related genes that are not significantly associated individually, and have been successful at detecting broad patterns of rare disruptive variation in other complex diseases such as schizophrenia [44]. We performed a primary examination of 33 biological pathways of potential relevance to suicidal behavior risk (online suppl. Table S4). The majority of these primary pathways were selected based on general neuronal development, maintenance, and function. In addition, several pathways were selected that have been specifically hypothesized to be involved with suicidal behavior risk, including the hypothalamic-pituitary-adrenal (HPA) axis, glutamatergic, and serotonergic pathways. None of the primary pathway association tests generated results that survived correction for multiple testing (p < 0.0015 for 33 primary pathways), with the best result being identified within extremely rare disruptive variants within the HPA axis pathway [55] with an OR of 6.5 and a p value of 0.040.

As no primary pathway survived correction for multiple testing, we chose to broaden our pathway investigation to include a total of 3,621 additional pathways obtained from the MSigDB [49]. Specifically, we selected well-described biological process pathways such as those found within KEGG [51] and GO [50] data sets, as well as pathways of genes that are regulated by specific transcription factors or microinhibitory RNAs. These tests generated 1 promising result, the KEGG Limonene and Pinene Degradation pathway (OR = 1.7, p = 1.6 × 10^-5, permutation corrected p = 0.10). Top contributing variants and genes for this pathway are presented in online supplementary Tables S8 and S9, respectively. In addition, several suggestive (p < 0.001) results were identified (Table 2) with 3 of these results, including our top result, being driven particularly by a set of 5 overlapping aldehyde dehydrogenase genes (KEGG Limonene and Pinene Degradation, KEGG Histidine Metabolism, and KEGG Beta Alanine Metabolism; p = 1.6 × 10^-5 to 4.7 × 10^-4; Fig. 2). As the aldehyde dehydrogenase genes serve several important metabolic roles, including functioning as key players in the metabolism of ingested alcohol, we included subject history of alcohol dependence as a covariate in the analyses. Correcting for alcohol dependence produced essentially no change in the magnitude of the signals for these pathways (p = 1.0 × 10^-5 to 6.9 × 10^-4).

This study represents a systematic, large-scale next-generation sequencing effort to examine the phenotype of suicidal behavior. The goal of this study was to provide an unbiased view of the contribution of exome-wide rare functional variation to suicidal behavior risk. This goal was accomplished through the use of a broad range of analysis tools on a variety of variant subsets in BP subjects with and without a history of suicide attempt.

Recent efforts to explore rare functional variation within severe psychiatric disorders, such as schizophrenia [44,56] have had success in identifying associations with disease risk. These findings incorporated a large number of rare variants that were too weak in association signal to be identified individually or even at the level of individual genes, but together implicated common biological pathways.

We chose to utilize a similar approach in order to maximize our analytical power. We performed a hypothesis-driven investigation of 33 pathways, including a number of general neuronal pathways and a selection of pathways particularly suspected to play a role within suicidal behavior risk (online suppl. Table S4). In addition, we performed pathway assessments on an additional 3,621 pathways primarily selected from the MSigDB [49]. Though this analysis did not identify results that could survive conservative Bonferroni correction for multiple testing, several potentially interesting pathways arose with 1 reaching near-significance, the KEGG [51] Limonene and Pinene Degradation pathway.

This pathway result, along with 2 other top pathway results, was driven by an enrichment of rare functional coding variants in suicide attempters within aldehyde dehydrogenase (ALDH) genes (Fig. 2). This is of potential interest to the field due to the well-established contribution of alcohol use to suicidal behavior risk [57]. Indeed, within the sample of attempters and nonattempters investigated in this study, alcohol dependence is significantly more common within attempters (OR = 2.1; p = 4.9 × 10^-8). This prompted us to include the history of alcohol dependence for each subject as a covariate within a secondary assessment of the top pathways. Including the alcohol dependence covariate resulted in essentially no change in the magnitude of the top pathways (Limonene and Pinene Degradation Pathway OR = 1.7, p = 1.0 × 10^-5) suggesting that alcohol dependence per se is not the primary driver of this signal trend.

It must be noted that none of the 33 primary pathways, including several hypothesized to play a role within suicidal behavior risk, assessed within this study showed strong evidence for association within our analyses. These primary pathways included the heavily researched serotonergic gene pathway. The serotonergic pathway has long been studied in suicidal behavior due to the central role this system plays in the treatment of psychiatric disease coupled with early evidence associating differential serotonin metabolite levels with suicidal behavior [58]. A large number of subsequent genetic studies have been carried out that have identified associations including SLC6A4, MAOA, TPH, TPH2, HTR1A, and HTR2A of the serotonergic system with suicidal behavior (see online suppl. Table S1).

There are several possibilities for why all of these previously associated serotonergic genes are not strongly enriched for functional variation within our exome data. The candidate gene studies that identified these associations were primarily focused on common variants that were either synonymous coding changes or fell within noncoding regions of the genome. These attributes are important for 2 reasons. First, the associated noncoding sites were not covered within our exome-sequencing targets. Second, synonymous variants would not be considered within our gene and pathway analyses due to our focus on sites with evidence for functional effects. Ultimately, additional expanded sequencing efforts to increase sample sizes as well as to encompass regulatory regions may be required to identify causal variants behind previously identified signals.

Existing high-throughput analyses in suicidal behavior, including GWAS, expression studies, and DNA methylation have identified a broad array of potentially important genes and regions, many of which are not among the set of candidate genes that has traditionally been the focus of suicidal behavior research efforts (online suppl. Table S1). For example, GWAS techniques have had success identifying nominally associated sites throughout the genome as well as 2 potential genome-wide significant associations (p ∼ 5.0 × 10^-8) with common variants in at least 1 population each within the 2p25 region [19] and the ABI3BP gene [17]. Genome-wide gene expression investigation efforts have similarly implicated pathways of potential importance such as polyamine regulation [6], general nervous system development [10], and glutamatergic and GABAergic synaptic signaling [11]. Genome-wide DNA methylation investigations in suicidal behavior have also identified general associations with many genes involved with cognitive processes [59] as well as the recent implication of SKA2 as a potential suicide candidate gene [60].

In like manner, our study sought to address the knowledge gap regarding the contribution of rare, functional variation throughout the human genome to suicidal behavior risk. A few of the suggestive results from this study show general overlap with some of the previously identified pathways, such as in the case of 1 of our top genes, SLC6A13, implicating GABAergic systems, and also the appearance of central nervous system development within our top pathway findings. Despite these examples, however, the results from this effort demonstrate minimal overlap with the genes and regions implicated by the above-referenced high-throughput efforts. This may be a reflection of fundamental differences in study design between the presented work and these other genome-wide investigations.

Specifically, the effects of rare functional sites are unlikely to have been robustly detected within the GWAS framework due to a focus on common variation sites. In addition, the top regions identified within suicidal behavior GWAS are found within regulatory regions not targeted within the design of our sequencing study. Similarly, methylation and expression analyses are focused on exploring the contribution of perturbed regulation of genes rather than coding variation directly damaging to gene function as was the focus of our analyses. Therefore, the variation examined within our whole-exome sequencing effort is largely unexamined by previous efforts, yielding a novel and valuable data resource that serves as a complement to existing efforts.

In total, 20% (99,303 sites) of the high-quality variants identified within our data set were previously undescribed in existing sequencing databases, including the 1000 Genomes [52], DBSNP version 142 [61], Exome Sequencing Project (http://evs.gs.washington.edu/EVS/), or the ExAC (http://exac.broadinstitute.org/). The majority of the novel sites (70%, 69,214 variants) are bioinformatically predicted to have potential functional affects and can be added to existing variant databases to be further assessed within future complex disease studies and meta-analyses of suicidal behavior. These resources, when combined with growing sample sizes and existing data, may assist in the creation of convergent lines of evidence that could lead to a more comprehensive and empirically supported understanding of the basis for suicidal behavior risk than is appreciated at this time.

Our study has several limitations. First, our sample had limited power to detect variants and gene signals of small to moderate effect size. We estimated 80% power to detect independent rare variants with a relative risk of 2.9 (MAF 0.05) to 7.2 (MAF 0.01) at study-wide significance (see online suppl. Fig. S3). We also estimated 80% power to detect a study-wide significant gene-burden signal with a relative risk of 2.1-3.9 depending on the collapsed functional variant frequency within the gene locus by subjects, noting that power levels will be reduced for burden testing if not all collapsed variants for a given gene contribute to risk (see online suppl. Fig. S4).

Second, the pathways we chose to examine represent a subset of all available pathways for analysis. Additionally, the examined pathways are not all-inclusive representations of the genes that may be important to the given processes, and may be missing key genes that would help explain that pathway's role in suicidal behavior. For example, one of the primary pathways we assessed, the KEGG [51] glutamatergic synapse pathway, did not include several key synaptic structural proteins such as NRXN or LRRTM family genes which have been implicated in suicidal behavior [19]. Furthermore, analyses geared toward any single pathway may fail to capture a signal that crosses several potentially associated pathways.

Third, while whole-exome sequencing offers an economical way to deeply examine the variation most likely to directly impact gene structure and function, this method suffers from certain disadvantages as well. Exome sequencing arrays are designed to target well-described transcripts within the genome. Therefore, exome sequencing data do not comprehensively cover all coding transcripts of the genome, potentially missing important regions. In addition, existing genome-wide studies of complex disease have frequently identified associated loci that are not within the coding regions of the genome, but are enriched for regions near coding genes [37]. While a fraction of the genes examined in this study did have limited coverage of regulatory regions as part of our customized array, our design largely did not examine the vast nonexomic regions of the genome. Further assessment of the whole genome or within regulatory regions near candidate genes may identify further key loci for suicidal behavior risk.

Fourth, our study design did not consider all available sources of captured variation within our gene and pathway-based analyses. For example, there is evidence that synonymous variation, which we did not include in these analyses, can contribute to disease risk [62]. In addition, sites of common variation (MAF >0.05) were excluded from these analyses. It remains to be seen how these sources of variation, combined with rare functional variation, contribute to suicidal behavior risk.

Fifth, we did not include consideration of environmental factors, such as individual abuse history, in our association tests as we did not have this information for all of the subjects. Such factors are known to increase the risk for suicidal behavior and may act synergistically with genetic variation to enhance risk.

Sixth, our study population was composed entirely of subjects with a history of BP. While such homogeneity is useful to limit signals associated with phenotypes other than suicidal behavior, it also limits the generalizability of our findings to those with a BP background. Additional cohorts with differing psychiatric backgrounds will need to be assessed in order to identify the generalizability of these findings to all cases of suicidal behavior. In addition, living subjects were sequenced for this study. It remains to be seen whether these findings can be replicated in those who have died from suicide.

This study represents the first effort to examine rare functional variation across the exome through sequencing in suicidal behavior. By using comprehensive analyses and investigating variation largely unexamined by existing studies, we have provided a view of the potential contribution of rare functional variation within suicidal behavior. Our results offer suggestive evidence of many genes and pathways that may prove to be informative for suicidal behavior risk. This study provides a wealth of novel information that can be considered within future investigations. Larger sample sizes and studies focused on whole-gene/-genome sequencing will further illuminate which signals are particularly important in suicide risk.

This work was supported by NIH grant R01MH079240 (Dr. Willour), using data generated by NIH grants R01MH087979 (Dr. Potash) and R01MH087992 (Dr. McCombie). Additional funding was provided by the University of Iowa Medical Scientist Training Program (MSTP) training grant 5 T32 GM007337 (E.T. Monson), American Foundation for Suicide Prevention (AFSP) PDF-0-067-12 (Dr. Breen), and NIH grant K01MH093809 (Dr. Pirooznia). We also wish to acknowledge the support of the University of Iowa Interdisciplinary Graduate Program in Genetics (E.T. Monson, S.C. Gaynor, Dr. Willour, and Dr. Potash). We gratefully acknowledge the preparation and distribution of samples used in this study by RUCDR Infinite Biologics. We also wish to express our gratitude to the many individuals who participated in the diagnosis and interviewing of subjects included in this study. Finally, we would like to thank the many families who offered their support and time to the study.

All subjects included in this study supplied institutional review board-approved written informed consent.

W.R. McCombie is a founding member of the plant genomics and cancer genetics company, Orion Genomics, and retains shares within this company. W.R. McCombie has also been provided costs associated with travel as well as honoraria for presenting at Illumina and Pacific Biosciences sponsored meetings over the past several years. None of these companies (Orion Genomics, Illumina, nor Pacific Biosciences) have played any part in the decision making, direct support, data generation, analysis, or any other part of this study. Mr. Monson, Dr. Pirooznia, Dr. Parla, Ms. Kramer, Dr. Goes, Dr. Breen, Ms. Gaynor, Ms. de Klerk, Dr. Jancic, Dr. Karchin, Dr. Zandi, Dr. Potash, and Dr. Willour reported no biomedical financial interests or potential conflicts of interest.