Subjects

The entire study design and procedures were performed in accordance with the Declaration of Helsinki. This study was approved by the Ethics Committee for Genetic Studies of Kobe University and RIKEN.

Individuals who died by suicide

Autopsies on suicide victims were conducted at the Division of Legal Medicine in the Department of Community Medicine and Social Health Science at the Kobe University Graduate School of Medicine. The verdict of “completed suicide” was made through discussion with the Medical Examiner’s Office of the Hyogo Prefecture and the Division of Legal Medicine in the Kobe University Graduate School of Medicine [17]. In order to gather background information on completed suicides, psychological autopsy through their medical records and bereaved family interviews were conducted by professional staff from the Medical Examiner’s Office of the Hyogo Prefecture and the Division of Legal Medicine in the Kobe University, where available.

Non-suicide controls

As non-suicide controls, we used genome-wide genotype data from subjects in the Biobank Japan project who had been genotyped as case subjects for non-psychiatric disorders and from healthy volunteers of the Osaka-Midosuji Rotary Club and the Pharma SNP consortium. The controls were not psychiatrically evaluated [18,19,20].

Genotyping, QC, and imputation

We genotyped 434 individuals who died by suicide between June 1996 and July 2012 and 405 individuals who died by suicide between August 2012 and February 2017 using Illumina HumanOmniExpress and HumanOmniExpressExome BeadChips for the first and second set, respectively. We obtained control data genotyped with the same arrays (N = 7993 and 7136 for the first and second set, respectively) (Table S1).

We performed quality control (QC) using PLINK 1.9 [21]. First, we excluded SNPs with a call rate <0.98 and minor allele frequency (MAF) <0.01, and those with P < 1.0 × 10−6 for Hardy–Weinberg equilibrium (HWE) in controls. Related individuals were excluded (PI_HAT ≥0.175). We then performed principal component analysis (PCA), and excluded samples outside the Japanese main islands cluster [22, 23]. The results of PCA are shown in Fig. S1. The final datasets included 386 suicides and 7458 controls as the first set, and 360 suicides and 6591 controls as the second set. After estimating haplotypes using SHAPEIT2 (v2.r778) [24], we performed genotype imputation by Minimac3 (1.0.13) [25] using ALL samples in the 1000 Genomes Project phase 3v5 [26] as a reference.

To investigate the X chromosome, we called genotypes using GenomeStudio. First, we generated genoplots using only female samples. After that, we added male samples to genoplots and called genotypes. Genotypes called as heterozygotes were treated as missing in male samples. For the genotyping QC, we excluded SNPs with MAF <0.01 and SNP call rate <0.98 in either male or female samples. We also excluded SNPs with P value for HWE <1.0 × 10−6 in female samples. Haplotype phasing and imputation were performed separately for males and females. Allelic dosages were imputed from 0 to 2 in male samples under the assumption of full dosage compensation. The pseudo-autosomal region was excluded from the reference before imputation.

After imputation of autosomal and X chromosomes, we included only SNPs with high imputation quality (r2 ≥ 0.7) and MAF ≥0.01 (Fig. S2).

Association analysis and meta-analysis

For all GWASs identifying genetic variants affecting suicide risk (case–control study; 746 suicides and 14,049 controls) and age at suicide (only suicidal cohorts; 719 suicides), single variant association tests were performed for common variants (MAF ≥0.01) using logistic regression and linear regression based on Wald test implemented in the Rvtests software [27] with a correction for the top 10 PCs as covariates. Meta-analysis was performed with the METAL software [28] using a fixed-effects model with inverse-variance weighted approach. The significance level was set at P < 2.5 × 10−8 due to correction for multiple comparison for two GWASs (case–control GWAS and GWAS for age at suicide). P for heterogeneity between two analyses (first set and second set) was calculated by Cochran’s Q test. Regional association plots were generated using LocusZoom [29]. A QQ plot for each GWAS is shown in Fig. S3.

Evaluation of previously implicated SNPs in prior GWAS of suicidal behavior

Among variants that have previously been associated with suicidal behavior (suicide ideation, suicide attempt, and suicide completion) in the published literature [7, 9, 30,31,32,33,34,35,36], we identified 26 variants with P < 1.0 × 10−6 in previous GWASs and MAF >0.01 in JPT (Japanese in Tokyo, Japan) population (the 1000 Genomes Project phase 3). We looked up the association of each of these SNPs with the results from our case–control GWAS. Only the top SNP in the same region in each reference was selected. For a candidate SNP that was not directly genotyped or imputed in our GWAS, we identified a proxy SNP (r2 > 0.8 in East Asian samples of the 1000 Genomes Project phase 3) whenever possible. For SNPs that reached P < 0.05, we determined whether the direction of the association was consistent between the prior and current studies.

Estimation of the proportion of the variance in completed suicide explained by the genotyped SNPs (SNP-based heritability)

To assess SNP-based heritability (h2 SNP ), we used genome-wide complex trait analysis (GCTA) [37] to generate genetic relatedness matrices (GRMs) among 506,645 SNPs that passed QC in both genotyped datasets, and then performed genomic restricted maximum-likelihood (GREML) analysis. This assumed prevalence rates of 0.1 and 0.5% for completed suicide, considering the reported incidence of completed suicide in Japan (Ministry of Health, Labor, and Welfare of Japan) [3] and the estimates from previous papers [38,39,40]. We strictly controlled the cryptic relatedness of the analyzed samples using the –grm-cutoff option (threshold of 0.05) implemented by GCTA, including 385 cases and 7409 controls for the first set, and 357 cases and 6560 controls for the second set, respectively. We then estimated SNP-based heritability using reml function implemented by GCTA with top 10 PCs as covariates.

PRS analysis

Polygenic risk score (PRS) analyses were performed using PRSice v1.23 [41]. The P threshold (P t ) for selecting “risk” SNPs was sequentially set at 0.1, 0.2, 0.3, 0.4, and 0.5 without SNPs in the major histocompatibility complex region. We then performed linkage disequilibrium (LD) clumping (used by the default setting of the software) to select the eligible SNPs for PRS. To calculate the PRS, we analyzed two discovery/target sets, using the first set (386 suicides and 7458 controls) as the discovery set and the second set (360 suicides and 6591 controls) as the target set, and vice versa. We included the top 10 PCs derived from each genotyped dataset as covariates, respectively. The variance explained for the PRS was estimated based on Nagelkerke’s R2 from a logistic regression model. For the analysis of polygenic effects on age at suicide (first set, 366 suicides; second set, 353 suicides), we also applied the same procedure as the above to the case–control GWAS datasets.

Pathway enrichment analysis

Using the results of a meta-analysis of the case–control GWAS datasets, we ran PASCAL [42] for gene-based enrichment analysis using 1077 gene sets, including KEGG [43], REACTOME [44], and BIOCARTA (http://cgap.nci.nih.gov/Pathways/BioCarta_Pathways). The significance level was set at χ2 P < 4.6 × 10−5 after Bonferroni correction for 1077 tests.