Biological Insights From 108 Schizophrenia-Associated Genetic Loci

128 Independent Associated Loci

We obtained genome-wide genotype data from which we constructed 49 ancestry matched, non-overlapping case-control samples (46 of European and three of East Asian ancestry, 34,241 cases and 45,604 controls) and 3 family-based samples of European ancestry (1,235 parent affected-offspring trios) (Supplementary Table 1 and Supplementary text). These comprise the primary PGC GWAS dataset. Genotypes from all studies were processed by the PGC using unified quality control procedures followed by imputation of SNPs and insertion-deletions using the 1000 Genomes Project reference panel.²⁵ In each sample, association testing was conducted using imputed marker dosages and principal components (PCs) to control for population stratification. The results were combined using an inverse-weighted fixed effects model.²⁶ After quality control (imputation INFO score ≥ 0.6, MAF ≥ 0.01, and successfully imputed in ≥ 20 samples), we considered around 9.5 million variants. The results are summarized in Figure 1. To enable acquisition of large samples, some groups ascertained cases via clinician diagnosis rather than a research-based assessment and provided evidence of the validity of this approach (Supplementary text).^11,13Post-hoc analyses revealed the pattern of effect sizes for associated loci was similar across different assessment methods and modes of ascertainment (Extended Data Figure 1), supporting our a priori decision to include samples of this nature.

For the subset of LD (linkage disequilibrium)-independent SNPs with P < 1×10⁻⁶ in the meta-analysis, we next obtained results from deCODE genetics (1,513 cases and 66,236 controls of European ancestry). We define LD-independent SNPs as those with low LD (r² < 0.1) to a more significantly associated SNP within a 500 kb window. Given high LD in the extended MHC region spans ~8 Mb, we conservatively include only a single MHC SNP to represent this locus. The deCODE data were then combined with those from the primary GWAS to give a dataset of 36,989 cases and 113,075 controls. In this final analysis, 128 LD-independent SNPs surpassed genome-wide significance (P ≤ 5×10⁻⁸) (Supplementary Table 2).

As in meta-analyses of other complex traits which identified large numbers of common risk variants, ^27,28 the test statistic distribution from our GWAS deviates from the null (Extended Data Figure 2). This is consistent with the previously documented polygenic contribution to schizophrenia. ^10,11 The deviation in the test statistics from the null (λ_GC = 1.47, λ1000 = 1.01) is only slightly less than expected (λ_GC = 1.56) under a polygenic model given fully informative genotypes, the current sample size, and the lifetime risk and heritability of schizophrenia.²⁹

Multiple additional lines of evidence allow us to conclude the deviation between the observed and null distributions in our primary GWAS indicates a true polygenic contribution to schizophrenia. First, applying a novel method ³⁰ that uses LD information to distinguish between the major potential sources of test statistic inflation, we found our results are consistent with polygenic architecture but not population stratification (Extended Data Figure 3). Second, for 78% of 263 LD-independent SNPs with P <1×10⁻⁶ in the case-control GWAS, the risk alleles from were overrepresented in cases in the independent samples from deCODE (P = 5.6 ×10⁻²¹). This degree of consistency between the primary GWAS and the replication data is highly unlikely to occur by chance (P = 1×10⁻⁹). (Note the tested alleles surpassed the P < 10⁻⁶ threshold in our GWAS before we added either the trios or deCODE data to the meta-analysis. This trend test is therefore independent of the primary case-control GWAS). Third, analysing the 1,235 parent-proband trios, we again found excess transmission of the schizophrenia-associated allele at 69% of the LD-independent SNPs with P <1×10⁻⁶ in the case-control GWAS. This is again unlikely to occur by chance (P = 1×10⁻⁹) and additionally excludes population stratification as fully explaining the associations that surpassed our threshold for seeking replication. Fourth, we used the trios trend data to estimate the expected proportion of true associations at P <1×10⁻⁶ in the discovery GWAS, allowing for the fact that half of the index SNPs are expected to show the same allelic trend in the trios by chance, and that some true associations will show opposite trends given the limited number of trio samples (Supplementary Material). Given the observed trend test results, around 67% (95% CI 64%-73%) or N = 176 of the associations in the scan at P <1×10⁻⁶ are expected to be true, and therefore the number of associations that will ultimately be validated from this set of SNPs will be considerably more than those that now meet genome-wide significance. Taken together, these analyses indicate that the observed deviation of test statistics from the null primarily represents polygenic association signal, and that the considerable excess of associations at the tail of extreme significance largely correspond to true associations.

Independently associated SNPs do not translate to well-bounded chromosomal regions. Nevertheless, it is useful to define physical boundaries for the SNP associations to identify candidate risk genes. We defined an associated locus as the physical region containing all SNPs correlated at r² > 0.6 with each of the 128 index SNPs. Associated loci within 250 kb of each other were merged. This resulted in 108 physically distinct associated loci, 83 of which have not been previously implicated in schizophrenia and therefore harbour potential novel biological insights into disease aetiology (Supplementary Table 3; regional plots in Supplementary Figure 1). The significant regions include all but 5 loci previously reported as genome-wide significant in large samples (Supplementary Table 3).

Characterization of Associated Loci

Of the 108 loci, 75% include protein-coding genes (40% a single gene) and a further 8% are within 20 kb of a gene (Supplementary Table 3). Notable associations relevant to major hypotheses of the aetiology and treatment of schizophrenia include DRD2 (the target of all effective antipsychotic drugs) and multiple genes (e.g. GRM3, GRIN2A, SRR, GRIA1) involved in glutamatergic neurotransmission and synaptic plasticity. In addition, associations at CACNA1C, CACNB2, and CACNA1I, which encode voltage-gated calcium channel subunits, extend previous findings implicating members of this family of proteins in schizophrenia and other psychiatric disorders. ^11,13,31,32 Genes encoding calcium channels, and proteins involved in glutamatergic neurotransmission and synaptic plasticity have been independently implicated in schizophrenia by studies of rare genetic variation, ^33-35 suggesting convergence at a broad functional level between studies of common and rare genetic variation. We highlight in the supplementary text genes of particular interest within associated loci with respect to current hypotheses of schizophrenia aetiology or treatment (although we do not imply these genes are necessarily the causal elements).

For each of the schizophrenia associated loci, we identified a credible causal set of SNPs (for definition see Supplementary text). ³⁶ In only 10 instances (Supplementary Table 4) was the association signal credibly attributable to a known nonsynonymous exonic polymorphism. The apparently limited role of protein coding variants is consistent both with exome sequencing findings ³³ and with the hypothesis that most associated variants detected by GWAS exert their effects through altering gene expression rather than protein structure^37,38 and with the observation that schizophrenia risk loci are enriched for expression quantitative trait loci (eQTL).³⁹

To try to identify eQTLs that could explain associations with schizophrenia, we merged the credible causal set of SNPs defined above with eQTLs from a meta-analysis of human brain cortex eQTL studies (N=550) and an eQTL study of peripheral venous blood (N=3754)⁴⁰ (Supplementary Text). Multiple schizophrenia loci contained at least one eQTL for a gene within 1 Mb of the locus (Supplementary Table 4). However, in only 12 instances was the eQTL plausibly causal (two in brain, and nine in peripheral blood, one in both). This low proportion suggests that if most risk variants are regulatory, available eQTL catalogs do not yet provide power, cellular specificity, or developmental diversity to provide clear mechanistic hypotheses for follow-up experiments.

Brain and Immunity

To further explore the regulatory nature of the schizophrenia associations, we mapped the credible sets (N=108) of causal variants onto sequences with epigenetic markers characteristic of active enhancers in 56 different tissues and cell lines (Supplementary Text). Schizophrenia associations were significantly enriched at enhancers active in brain (Figure 2) but not in tissues unlikely to be relevant to schizophrenia (e.g., bone, cartilage, kidney, and fibroblasts). Brain tissues used to define enhancers consist of heterogeneous populations of cells. Seeking greater specificity, we contrasted genes enriched for expression in neurons and glia using mouse ribotagged lines⁴¹. Genes with strong expression in multiple cortical and striatal neuronal lineages were were enriched for associations, providing support for an important neuronal pathology in schizophrenia (Extended Data Figure 4) but this is not statistically more significant than, or exclusionary of, contributions from other lineages⁴².

Schizophrenia associations were also strongly enriched at enhancers that are active in tissues with important immune functions, particularly B-lymphocyte lineages involved in acquired immunity (CD19 and CD20 lines, Figure 2). These enrichments remain significant even after excluding the extended MHC region and regions containing brain enhancers (enrichment P for CD20 <10⁻⁶), demonstrating that this finding is not an artefact of correlation between enhancer elements in different tissues and not driven by the strong and diffuse association at the extended MHC. Epidemiological studies have long hinted at a role for immune dysregulation in schizophrenia, the present findings provide genetic support for this hypothesis. ⁴³

Aiming to develop additional biological hypotheses beyond those that emerge from inspection of the individual loci, we further undertook a limited mining of the data through gene-set analysis. However, since there is no consensus methodology by which such analyses should be conducted, nor an established optimal significance threshold for including loci, we sought to be conservative, using only two of the many available approaches^44,45 and restricting analyses to genes within genome-wide significant loci. Neither approach identified gene-sets that were significantly enriched for associations after correction for the number of pathways tested (Supplementary Table 5) although nominally significantly enrichments were observed among several predefined candidate pathways (Extended Data Table 1). A fuller exploratory analysis of the data will be presented elsewhere.

Overlap with rare mutations

CNVs associated with schizophrenia overlap with those associated with ASD and intellectual disability⁹ as do genes with deleterious de novo mutations. ³⁴ Here, we find significant overlap between genes in the schizophrenia GWAS associated intervals and those with de novo nonsynonymous mutations in schizophrenia (P=0.0061) (Extended Data Table 2), suggesting that mechanistic studies of rare genetic variation in schizophrenia will be informative for schizophrenia more widely. We also find evidence for overlap between genes in schizophrenia GWAS regions and those with de novo nonsynonymous mutations in intellectual disability (P=0.00024) and ASD (P=0.035), providing further support for the hypothesis that these disorders have partly overlapping pathophysiologies.^9,34

Polygenic Risk Score Profiling

Previous studies have shown that risk profile scores (RPS) constructed from alleles showing modest association with schizophrenia in a discovery GWAS can predict case-control status in independent samples, albeit with low sensitivity and specificity.^10,11,16 This finding was robustly confirmed in the present study. The estimate of Nagelkerke R² (a measure of variance in case-control status explained) depends on the specific target dataset and threshold (P_T) for selecting risk alleles for RPS analysis (Extended Data Figure 5 & 6a). However, using the same target sample as earlier studies and P_T =0.05, R² is now increased from 0.03 (reference¹⁰) to 0.184 (Extended Data Figure 5). Assuming a liability-threshold model, a lifetime risk of 1%, independent SNP effects, and adjusting for case-control ascertainment, RPS now explains about 7% of variation on the liability scale⁴⁶ to schizophrenia across the samples (Extended Data Figure 6b), about half of which (3.4%) is explained by genome-wide significant loci.

We also evaluated the capacity of RPS to predict case-control status using a standard epidemiological approach to a continuous risk factor. We illustrate this in three samples, each with different ascertainment schemes (Figure 3). The Danish sample is population-based (i.e. inpatient and outpatient facilities), the Sweden sample is based on all cases hospitalized for schizophrenia in Sweden, and the MGS study was ascertained specially for genetic studies from clinical sources in the US and Australia. We grouped individuals into RPS deciles and estimated the odds ratios for affected status for each decile with reference to the lowest risk decile. The odds ratios increased with greater number of schizophrenia risk alleles in each sample, maximizing for the tenth decile in all samples: Denmark 7.8 (95% CI 4.4-13.9), Sweden 15.0 (95% CI 12.1-18.7), and MGS 20.3 (95% CI 14.7-28.2). Given the need for measures that index liability to schizophrenia,^47,48 the ability to stratify individuals by RPS offers new opportunities for clinical and epidemiological research. Nevertheless, we stress that the sensitivity and specificity of RPS do not support its use as a predictive test. For example in the Danish epidemiological sample, the area under the receiver operating curve is only 0.62 (Extended Data Figure 6c, Supplementary Table 6).

Finally, seeking evidence for non-additive effects on risk, we tested for statistical interaction between all pairs of 125 autosomal SNPs that reached genome-wide significance. P-values for the interaction terms were distributed according to the null, and no interaction was significant after correction for multiple comparisons. Thus, we find no evidence for epistatic or non-additive effects between the significant loci (Extended Data Figure 7). It is possible that such effects could be present between other loci, or occur in the form of higher order interactions.

Discussion

In the largest molecular genetic study of schizophrenia ever conducted, or indeed of any neuropsychiatric disorder, we demonstrate the power of GWAS to identify large numbers of risk loci. We also show that the use of alternative ascertainment and diagnostic schemes designed to rapidly increase sample size does not inevitably introduce a crippling degree of heterogeneity. That this is true for a phenotype like schizophrenia where there are no biomarkers or supportive diagnostic tests provides grounds to be optimistic about applying GWAS to other disorders where diagnosis is similarly clinician-based, and recruitment traditionally constrained by the need for expensive and time consuming assessments.

We further show that the associations are not randomly distributed across genes of all classes and function; rather they converge upon genes that are expressed in certain tissues and cellular types. The findings include molecules that are the current, or the most promising, targets for therapeutics, and point to systems that align with the predominant aetiological hypotheses of the disorder. This suggests that the many novel findings we report also provide a aetiologically relevant foundation for mechanistic and treatment development studies. We additionally find overlap between genes affected by rare variants in schizophrenia and those within GWAS loci, and broad convergence in the functions of some of the clusters of genes implicated by both sets of genetic variants, particularly genes related to abnormal glutamatergic synaptic and calcium channel function. How variation in these genes impact function to increase risk for schizophrenia cannot be answered by genetics, but the overlap strongly suggests that common and rare variant studies are complementary rather than antagonistic, and that mechanistic studies driven by rare genetic variation will be informative for schizophrenia.

Extended Data

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=β0+β1X1+β2X2+β3∗X1∗X2+β4X4+β5X5X₁ 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.

Supplementary Material