Biological insights from 108 schizophrenia-associated genetic loci

Abstract

Schizophrenia is a highly heritable disorder. Genetic risk is conferred by a large number of alleles, including common alleles of small effect that might be detected by genome-wide association studies. Here we report a multi-stage schizophrenia genome-wide association study of up to 36,989 cases and 113,075 controls. We identify 128 independent associations spanning 108 conservatively defined loci that meet genome-wide significance, 83 of which have not been previously reported. Associations were enriched among genes expressed in brain, providing biological plausibility for the findings. Many findings have the potential to provide entirely new insights into aetiology, but associations at DRD2 and several genes involved in glutamatergic neurotransmission highlight molecules of known and potential therapeutic relevance to schizophrenia, and are consistent with leading pathophysiological hypotheses. Independent of genes expressed in brain, associations were enriched among genes expressed in tissues that have important roles in immunity, providing support for the speculated link between the immune system and schizophrenia.

Extended Data Figure 1. Homogeneity of Effects Across Studies

Plot of the first two principal components (PC) from principal components analysis (PCA) of the logistic regression β coefficients for autosomal genome-wide significant associations. The input data were the β coefficients from 52 samples for 112 independent SNP associations (excluding 3 chrX SNPs and 13 SNPs with missing values in Asian samples). PCAs were weighted by the number of cases. Each circle shows the location of a study on PC1 and PC2. Circle size and colour are proportional to the number of cases in each sample (larger and redder circles correspond to more cases). Most samples cluster. Outliers had either small numbers of cases (“small”) or were genotyped on older arrays. Abbreviations. a500 (Affymetrix 500K); a5 (Affymetrix 5.0). Studies that did not use conventional research interviews are in the central cluster (CLOZUK, Sweden, and Denmark-Aarhus studies).

Extended Data Figure 2. Quantile-quantile plot

Quantile-quantile plot of GWAS meta-analysis. Expected −log10(P) -values are those expected under the null hypothesis. For clarity, we avoided expansion of the Y-axis by setting association P-values < 10⁻¹² to 10⁻¹². The shaded area surrounded by a red line indicates the 95% confidence interval under the null. Lambda is the observed median χ² test statistic divided by the median expected χ² test statistic under the null.

Extended Data Figure 3. LD Score regression consistent with polygenic inheritance

The relationship between marker χ² association statistics and linkage disequilibrium (LD) as measured by the LD Score. LD Score is the sum of the r² values between a variant and all other known variants within a 1 cM window, and quantifies the amount of genetic variation tagged by that variant. Variants were grouped into 50 equal-sized bins based on LD Score rank. LD Score Bin and Mean χ² denotes mean LD score and test statistic for markers each bin. We simulated (Supplementary Text) test statistics under two scenarios: (a) no true association, inflation due to population stratification (b) polygenic inheritance (λ_GC =1.32), in which we assigned independent and identically distributed per-normalized-genotype effects to a randomly selected subset of variants. Panel (c) present results from the PGC schizophrenia GWAS (λ_GC =1.48). The real data are strikingly similar to the simulated data summarized in (b) but not (a). The intercept estimates the inflation in the mean χ² that results from confounding biases, such as cryptic relatedness or population stratification. Thus, the intercept of 1.066 for the schizophrenia GWAS suggests that ~90% of the inflation in the mean χ² results from polygenic signal. The results of the simulations are also consistent with theoretical expectation (see Supplementary Text).

Extended Data Figure 4. Enrichment of Associations in Tissues and Cells

Genes whose transcriptional start is nearest to the most associated SNP at each schizophrenia-associated locus were tested for enriched expression in a) purified brain cell subsets obtained from mouse ribotagged lines. The red dotted line indicates P=0.05.

Extended Data Figure 5. MGS Risk Profile Score Analysis

Polygenic risk profile score (RPS) analyses using the MGS sample as target, and deriving risk alleles from three published schizophrenia datasets (X-axis): ISC (2615 cases and 3338 controls) , PGC1 (excluding MGS, 9320 cases and 10,228 controls) , and the current meta analysis (excluding MGS) with 32,838 cases and 44,357 controls. Samples sizes differ slightly from original publication due to different analytical procedures. This shows the increasing RPS prediction with increasing training dataset size reflecting improved precision of estimates of the SNP effect sizes. The proportion of variance explained (Y-axis; Nagelkerke’s R²) was computed by comparison of a full model (covariates + RPS) score to a reduced model (covariates only). Ten different P-value thresholds (P_T) for selecting risk alleles are denoted by the colour of each bar (legend above plot). For significance testing, see the bottom legend which denotes the P-value for the test that R² is different from zero. All numerical data and methods used to generate these plots are available in Supplementary Tables 6, 7, and Supplementary Text.

Extended Data Figure 6. Risk Profile Score Analysis

We defined 40 target subgroups of the primary GWAS dataset and performed 40 leave-one-out GWAS analyses (see Supplementary Material) from which we derived risk alleles for RPS analysis (X-axis) for each target subgroup. a) The proportion of variance explained (Y-axis; Nagelkerke’s R²) was computed for each target by comparison of a full model (covariates + RPS) score to a reduced model (covariates only). For clarity, 3 different P-value thresholds (P_T) are presented denoted by the colour of each bar (legend above plot) as for Extended Data Figure 5 but for clarity we restrict to fewer P-value thresholds (P_T of 5×10⁻⁸, 1×10⁻⁴, and 0.05) and removed the significance values. (b) The proportion of variance on the liability scale from risk scores calculated at the P_T 0.05 with 95% CI bar assuming baseline population disease risk of 1%. (C) Area under the receiver operating curve (AUC). All numerical data and methods used to generate these plots are available in Supplementary Tables 6, 7, and Supplementary Text.

Extended Data Figure 7. Epistasis

Quantile-quantile plot for all pair-wise (N=7750) combinations of the 125 independent autosomal genome-wide significant SNPs tested for non-additive effects on risk using case-control datasets of European ancestry (32,405 cases and 42,221 controls). We included as covariates the principal components from the main analysis as well as a study indicator. The interaction model is described by: $$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3\ast {X_1}^{\ast }{X}_2+{\beta }_4{X}_4+{\beta }_5{X}_5$$ X₁ and X₂ are genotypes at the two loci, X₁*X₂ is the interaction between the two genotypes modeled in a multiplicative fashion, X₄ is the vector of principal components, X₅ is the vector of study indicator variables. Each β is the regression coefficient in the generalized linear model using logistic regression. The overall distribution of P-values did not deviate from the null and the smallest P-value (4.28×10⁻⁴) did not surpass the Bonferroni correction threshold (p=0.05/7750= 6.45×10⁻⁶). The line x=y indicates the expected null distribution with the grey area bounded by red lines indicating the expected 95% confidence interval for the null.

Figure 1. Manhattan plot

Manhattan plot of the discovery genome-wide association meta-analysis of 49 case control samples (34,241 cases and 45,604 controls) and 3 family based association studies (1,235 parent affected-offspring trios). The x-axis is chromosomal position and the y-axis is the significance of association (−log₁₀(P)).The red line shows the genome-wide significance level (5×10⁻⁸). SNPs in green are in LD with the index SNPs (diamonds) which represent independent genome-wide significant associations.

Figure 2. Enrichment in enhancers

Cell and tissue type specific enhancers were identified using ChIP-seq datasets (H3K27ac signal) from 56 cell line and tissue samples (Y-axis). We defined cell and tissue type enhancers as the top 10% of enhancers with the highest ratio of reads in that cell or tissue type divided by the total number of reads. Enrichment of credible causal associated SNPs from the schizophrenia GWAS was compared with frequency matched sets of 1000 Genomes SNPs (Supplementary Text). The X-axis is the −log₁₀(P) for enrichment. P-values are uncorrected for the number of tissues/cells tested. A −log₁₀(P) of roughly 3 can be considered significant after Bonferroni correction. Descriptions of cell and tissue types at the Roadmap Epigenome website (http://www.roadmapepigenomics.org).

Figure 3. Odds ratio by risk score profile

Odds ratio for schizophrenia by risk score profile (RPS) decile in the Sweden (Sw1-6), Denmark (Aarhus), and Molecular Genetics of Schizophrenia studies (Supplementary text). Risk alleles and weights were derived from ‘leave one out’ analyses in which those samples were excluded from the GWAS meta-analysis (Supplementary text). The threshold for selecting risk alleles was P_T < 0.05. The RPS were converted to deciles (1=lowest, 10=highest RPS), and nine dummy variables created to contrast deciles 2-10 to decile 1 as the reference. Odds ratios and 95% confidence intervals (bars) were estimated using logistic regression with PCs to control for population stratification.