Genetics of educational attainment aid in identifying biological subcategories of schizophrenia

Higher educational attainment (EA) is known to have a protective effect regarding the severity of schizophrenia (SZ). However, recent studies have found a small positive genetic correlation between EA and SZ. Here, we investigate possible causes of this counterintuitive finding using genome-wide association results for EA and SZ ( n = 443,581) and a replication cohort (1,169 controls and 1,067 cases) with high-quality SZ phenotypes. We find strong genetic overlap between EA and SZ that cannot be explained by chance, linkage disequilibrium, or assortative mating. Instead, our results suggest that the current clinical diagnosis of SZ comprises at least two disease subtypes with non-identical symptoms and genetic architectures: One part resembles bipolar disorder (BIP) and high intelligence, while the other part is a cognitive disorder that is independent of BIP.


INTRODUCTION
Schizophrenia (SZ) is the collective term used for a severe, highly heterogeneous and costly psychiatric disorder that is caused by a complex interplay of environmental and genetic factors [1][2][3][4] . The latest genome-wide association study (GWAS) by the Psychiatric Genomics Consortium (PGC) identified 108 genomic loci that are associated with SZ 5 . These 108 loci jointly account for ≈3.4% of the variation on the liability scale to SZ 5 , while all single nucleotide polymorphisms (SNPs) that are currently measured by SNP arrays capture ≈64% (s.e. = 8%) in the variation in liability to the disease 6 . This suggests that most of the genetic variants contributing to the heritability of SZ have very small effects and that they have not been isolated yet. This could be due in part to the fact that the clinical disease classification of SZ spans across many different syndromes (e.g., catatonia, paranoia, grandiosity, difficulty in abstract thinking, thought blocking, social withdrawal, hallucinations) that may not have identical genetic architectures. Therefore, identifying additional genetic variants and understanding through which pathways they influence specific SZ syndromes is an important step in understanding the etiologies of the 'schizophrenias' 7 . However, GWAS analyses of specific SZ syndromes would require very large sample sizes to be statistically well-powered, and the currently available datasets on deeply phenotyped SZ individuals are not large enough yet for this purpose.
Here, we use an alternative approach that combines data for SZ with another cognitive phenotype that can be studied in very large GWAS samples-educational attainment (EA). The relationship between SZ and EA is peculiar: There are contradictory results on the relationship between SZ and EA from phenotypic and genetic data that can be used as an avenue to further our understanding about SZ. Phenotypic data seem to suggest a negative correlation between EA and SZ 8 . For example, SZ patients with lower EA typically show an earlier age of disease onset, higher levels of psychotic symptomatology, and worsened global cognitive function 8 . In fact, EA has been suggested to be a measure of premorbid function and predictor of outcomes in SZ. Moreover, it has been forcefully argued that retarded intellectual development during childhood and bad school performance should be seen as core features of SZ and early indicators of the disease that precede the development of psychotic symptoms 9,10 . Furthermore, credible genetic links between SZ and impaired cognitive performance have been found 11 .
In contrast to these findings, recent studies using GWAS results identified a small, but positive genetic correlation between EA and SZ 12,13 . Here, we explore possible reasons for this contradictory result using the largest, non-overlapping GWAS samples on cognitive traits to date, totaling 443,581 individuals of European descent (the vast majority of observations coming from EA). For follow-up analyses, we use data from an independent replication sample that has exceptionally detailed measures of SZ symptoms, the GRAS (Göttingen Research Association for Schizophrenia) data collection 4,7,14 .
As a first step, we used the proxy-phenotype method (PPM) to illustrate the genetic overlap between EA and SZ. As a side-result, this approach may isolate novel empirically plausible candidate genes for SZ, comparable to similar studies using PPM that have demonstrated this for cognitive performance 15 , Alzheimer's disease, intracranial and hippocampal volume 13 , depression and neuroticism 16 . PPM is a two-stage approach that increases statistical power by using genetic association results from a large, independent sample for a related phenotype to limit the multiple testing burden for the phenotype of interest 15 . Previous evidence suggests a strong genetic overlap between EA and SZ, which implies that EA could be used as a proxyphenotype for SZ because EA can be studied in much larger samples 13,16 . However, compared to the present work, these previous studies used substantially smaller and partially overlapping samples and did not have access to an independent cohort that could be used for replication and follow-up analyses.
There are several possible reasons why EA-associated SNPs may also be associated with SZ. One possibility is that a set of genes that is generally important for all brain-related phenotypes is driving this enrichment. This hypothesis suggests that the set of genetic loci that our proxy-phenotype analysis identifies should be generally enriched for association with all brain-related phenotypes, but not for non-brain-related outcomes. To investigate this possibility, we test genetic loci that are jointly associated with EA and SZ for enrichment across 21 additional traits (Supplementary Note).
Second, enrichment could also be a generic consequence of EA-associated SNPs exhibiting above average linkage disequilibrium (LD) with neighboring SNPs. This would increase the probability that these SNPs "tag" other genetic variants that are associated with SZ, or any other disorder 12 . To investigate this possibility, we propose a measure that tests for enrichment beyond what is expected for each EA related SNP given its LD with its neighbors (Supplementary Note).
A third possible cause of strong enrichment and weak genetic correlation is heterogeneity in SZ-i.e., sub-types of the disease having different biological causes and varying genetic correlations with EA. Heterogeneity in the disease may also be a reason why previous studies did not succeed in predicting specific syndromes of SZ using a "normal" polygenic score (PGS) that was derived from large-scale GWAS on SZ, which implicitly assumed that all SZassociated SNPs influence all syndromes in the same way 4,17 . If heterogeneity in the disease is causing the observed enrichment of EA with SZ, the sign concordance pattern of SNPs with both traits may contain relevant information that is pertinent to specific SZ syndromes. We tested this by constructing PGS in our replication cohort with high-quality SZ phenotypes that take the sign concordance of SNPs for EA and SZ into account (Supplementary Note). As a robustness check, we repeat this analysis excluding patients diagnosed with schizoaffective disorder.
A fourth possible cause of enrichment is that other phenotypes are genetically correlated with both EA and SZ. Previous studies indicated a particularly strong positive genetic correlation between SZ and bipolar disorder (BIP), which may influence the genetic overlap of both diseases with related phenotypes such as EA, childhood intelligence (IQ), and neuroticism 12,13,18 . We use genome-wide inferred statistics (GWIS) that allow controlling for the genetic correlation between SZ and BIP to investigate how "unique" SZ (controlling for BIP) and "unique" BIP (controlling for SZ) are related to EA, childhood IQ, and neuroticism 18 .
A fifth possible cause may be assortative mating, which has been demonstrated both for EA 19 and SZ 20 . We use simulations to explore if independent assortative mating for the two phenotypes may induce a spurious genetic overlap.
This list of potential causes for the genetic overlap between EA and SZ may not be exhaustive and several of these factors may be at work simultaneously. Figure 1 presents an overview of the proxy-phenotype analyses. The first-stage GWAS on EA (Supplementary Note) identified 506 loci that passed our predefined threshold of P EA < 10 −5 (https://osf.io/dnhfk/); 108 of them were genome-wide significant (P EA < 5 × 10 −8 , see Supplementary Table 5.1). Of the 506 EA lead-SNPs, 132 are associated with SZ at nominal significance (P SZ < 0.05), and 21 of these survive Bonferroni correction (P SZ < 0.05 506 = 9.88 × 10 −5 ). LD score regression results suggest that the vast majority of the association signal in both the EA 13 and the SZ 5 GWAS are truly genetic signals, rather than spurious signals originating from uncontrolled population stratification. Figure 2a shows a Manhattan plot for the GWAS on EA highlighting SNPs that were also significantly associated with SZ (red crosses for P SZ < 0.05, green crosses for P SZ = 9.88 × 10 −5 ).

Proxy-phenotype analyses
A Q-Q plot of the 506 EA lead SNPs for SZ is shown in Figure 2b. Although the observed sign concordance of 52% is not significantly different from a random pattern ( = 0.40), we find 3.23 times more SNPs in this set of 506 SNPs that are nominally significant for SZ than expected given the distribution of the P values in the SZ GWAS results (raw enrichment = 6.87 × 10 −10 , Supplementary Note). The observed enrichment of the 21 EA lead SNPs that pass Bonferroni correction for SZ ( < 0.05 506 = 9.88 × 10 −5 ) is even more pronounced (27 times stronger, = 5.44 × 10 −14 ).

Bayesian credibility of the results
The effect sizes of these 21 SNPs on SZ are small, ranging from Odds = 1.02 (rs4500960) to Odds = 1.11 (rs4378243) after winner's curse correction (Table 1). However, Bayesian calculations with reasonable prior beliefs (e.g., 1% or 5%, Supplementary Note) suggest that most of these 21 SNPs are likely or virtually certain to be truly associated with SZ.

Prediction of future genome-wide significant loci for schizophrenia
Of the 21 variants we identified, 12 are in LD with loci previously reported by the PGC 5 and 2 are in the major histocompatibility complex (MHC) region on chromosome 6 and were therefore not separately reported in that study. Three of the variants we isolated (rs7610856, rs143283559, rs28360516) were independently found in a recent meta-analysis of the PGC results 5 with another large-scale sample 21 . We show in the Supplementary Note that using EA as a proxy-phenotype for SZ helped to predict the novel genome-wide significant findings reported in that study, which illustrates the power of the proxy-phenotype approach. Furthermore, two of the 21 variants (rs756912, rs7593947) are in LD with loci recently reported in a study that also compared GWAS findings from EA and SZ using smaller samples and a less conservative statistical approach 22 (Supplementary Note). The remaining 2 SNPs we identified (rs7336518 on chr13 and rs7522116 on chr 1) add to the list of empirically plausible candidate loci for SZ. Figure 3 and Supplementary Table 5. 2 show the LD-aware enrichment of the SNPs that are jointly associated with EA and SZ across 22 traits. We find significant joint LD-aware enrichment of this set of SNPs for SZ, BIP, neuroticism and childhood IQ, and for inflammatory bowel disease and age at menarche. However, we find no LD-aware enrichment for other brain-traits that are phenotypically related to SZ, such as depressive symptoms, subjective well-being, autism, and attention deficit hyperactivity disorder. We also do not find LD-aware enrichment for most traits that are less obviously related to the brain (e.g., BMI, coronary artery disease) and our negative controls (e.g., fasting insulin, birth weight, birth length). Furthermore, one of the novel SNPs we isolated shows significant LDaware enrichment both for SZ and for BIP (rs7522116).

Replication in the GRAS sample
A PGS based on the 132 loci jointly associated with both EA and SZ (SZ_132) adds ∆ 2 = 7.54% − 7.01% = 0.53% predictive accuracy for the SZ case-control status to a PGS (SZ_all) derived from the GWAS on SZ alone ( = 1.7 × 10 −4 , Table 2, Model 3). The SZ_132 score also significantly adds ( = 3.4 × 10 −4 ) to the predictive accuracy of the SZ case-control status when all other scores we constructed are included as control variables. In addition to SZ_132, PGS for SZ (SZ_all) and for BIP (BIP_all ) also predict case-control status, jointly reaching an adjusted ΔR 2 of ≈ 9% ( Table 2, Model 9 and Supplementary Note).

Polygenic prediction of schizophrenia measures in the GRAS sample
We find that the number of years of education is phenotypically correlated with later age at prodrome, later onset of disease, and less severe disease symptoms among SZ patients in the GRAS sample (Supplementary Note, Supplementary Table 8.1 and Supplementary Fig.  1). The EA_all score is associated with years of education ( = 2.6 × 10 −6 ) and premorbid IQ ( = 2.3 × 10 −4 ) among SZ patients (Supplementary Note and Table 3). Consistent with earlier results 4 , we find that none of the SZ measures can be predicted by the normal SZ PGS (SZ_all, Supplementary Table 8.2). Importantly, by utilizing GWAS results from both EA and SZ, we show that it is possible to predict specific features of SZ (Global Assessment of Functioning (GAF), Clinical Global Impression of Severity (CGI-S), and Positive and Negative Syndrome Scale (PANSS)) from genetic data. In a multiple regression analysis 23 that allows a "ceteris paribus" interpretation of the included variables, we find that the EA_all score is associated with less severe disease outcomes only if we condition on the effects of the Concordant and Discordant scores. And conditional on the EA_all score, the Concordant and Discordant scores are associated with more (less) severe positive and negative symptoms as measured by the PANSS scale, respectively ( Table 3). The best predictive accuracy of SZ readouts using these scores is currently observed for GAF (R 2 = 1.38%). Of note, several of the symptoms measured by PANSS are also symptoms of BIP. The degree and composition of symptoms varies with the phase at evaluation (manic or depressive) and the general disease severity. We repeated these analyses excluding patients who were diagnosed with schizoaffective disorder (SD) and found similar results, implying that the genetic heterogeneity in SZ that we identify is not only due to SD (Supplementary Note, Supplementary Table 8

Controlling for the genetic overlap between schizophrenia and bipolar disorder
None of the EA-associated lead SNPs (P EA < 5 × 10 −8 ) are significantly associated with "unique" SZ(min BIP) after Bonferroni correction (Supplementary Table 9.1, Supplementary Note). The sign concordance of the EA lead SNPs with "unique" SZ(min BIP) was 44.5% ( = 0.046). Supplementary Figure 2 shows a Q-Q plot of the EA lead-SNPs for "unique" SZ(min BIP). Although we find 1.6 times more EA-associated SNPs with < 0.05 than expected by chance (raw enrichment = 0.02, Supplementary Note), the enrichment is much weaker than in the main SZ GWAS results that did not control for the genetic overlap between SZ and BIP. The genetic correlations between EA SZ(min BIP), and IQ and SZ(min BIP) are negative and significant (rg = -0.16, P = 3.88×10 -04 and rg = -0.31, P = 6.00×10 -03 respectively), which is in line with the idea of SZ being a cognitive disorder 9 . Furthermore, the genetic correlations of EA and IQ with BIP(min SZ) remain positive and get somewhat stronger (rg = 0.31, P = 2.87×10 -07 and rg = 0.33, P = 3.18×10 -02 respectively) compared with the ordinary BIP GWAS results. However, controlling for the genetic overlap of SZ and BIP does not affect the genetic correlations with neuroticism (Figure 4).

Simulations of assortative mating
Our simulations for assortative mating were based on relatively extreme assumptions that increased our chance of finding spurious enrichment of EA loci for SZ. The results suggest it is unlikely that assortative mating is a major cause for the genetic overlap we observe between EA and SZ (Supplementary Fig. 3).

Biological annotations
Biological annotation of the 132 SNPs that are jointly associated with EA and SZ using DEPICT identified 111 significant reconstituted gene sets (Supplementary  Figure 5a). All significantly enriched tissues are related to the nervous system and sense organs (Figure 5b). Furthermore, "Neural Stem Cells" is the only significantly enriched cell-type (Supplementary Table 10.3). DEPICT prioritized genes that are known to be involved in neurogenesis and synapse formation (Supplementary Table  10.4). Some of the genes, including SEMA6D and CSPG5, have been suggested to play a potential role in SZ 24,25 . For the two novel candidate SNPs reported in this study (rs7522116 and rs7336518), DEPICT points to the FOXO6 (Forkhead Box O6) and the SLITRK1 (SLIT and NTRK Like Family Member 1) genes, respectively. FOXO6 is predominantly expressed in the hippocampus and has been suggested to be involved in memory consolidation, emotion and synaptic function 26,27 . Similarly, SLITRK1 is also highly expressed in the brain 28 , particularly localized to excitatory synapses and promoting their development 29 , and it has previously been suggested to be a candidate gene for neuropsychiatric disorders 30 .

DISCUSSION
We explored the genetic overlap between EA and SZ using the largest currently available GWAS sample on human cognitive traits to date. Using EA as a proxy-phenotype, we identified 21 genetic loci for SZ and showed that this approach helps to predict future GWAS hits for SZ. We isolated two additional candidate genes for SZ, FOXO6 and SLITRK1. Our results show that EA-associated SNPs are much more likely to also be associated with SZ than expected by chance. However, these genetic loci do not influence both traits with a systematic sign pattern that would correspond to a strong positive or negative genetic correlation.
The results of our follow-up analyses are most consistent with two hypotheses that complement each other: First, the genetic overlap between EA and SZ is to some extent induced by pleiotropic effects of many genes that affect not only EA and SZ but also other traits such as BIP and IQ. Second, different syndromes of SZ (e.g., low cognitive performance and psychosis) seem to be driven by different genetic effects. The clinical diagnosis of SZ aggregates over these different syndromes. In particular, our results suggest that the current clinical diagnosis of SZ comprises at least two disease subtypes with nonidentical symptomatology and genetic architectures: One part resembles bipolar disorder (BIP) and high intelligence, while the other part is a cognitive disorder that is independent of BIP. Consistent with this idea, we find that PGS that take the sign concordance of SNPs with EA and SZ into account begin for the first time to predict specific SZ features from genetic data (R 2 between 0.4% and 1.4%), while this was not possible with "ordinary" PGS for SZ.
Other mechanisms that we explored, in particular LD-patterns of the EA-associated SNPs and assortative mating, do not seem to be major drivers of the genetic overlap between EA and SZ. Furthermore, the loci we identified in our PPM analysis do not seem to be associated with all brain-related phenotypes, suggesting some degree of phenotype-specificity of the results. We note that the enrichment for age at menarche of the SNPs that are jointly associated with EA and SZ may be related to the final stage of brain development which coincides with the onset of puberty [31][32][33][34] .
The highly complex genetic architecture of the "schizophrenias" that our results point to implies that most patients will have individual-specific genetic loads for either subtype of the disease, contributing to individual differences in symptoms. The genetic heterogeneity we identified could imply that treatments will vary in their effectiveness across disease subtypes.
Overall, our study corroborates that EA is a useful proxy-phenotype for psychiatric outcomes. Specifically, combining GWAS results from EA and SZ led to the identification of two seemingly distinct subcategories of SZ. Even though each of them may still harbor highly heterogeneous disease subgroups, the new subcategories can pave the way for further biological subgroup analyses. Therefore, a psychiatric nosology that is based on biological causes rather than pure phenotypical classifications may be feasible in the future. Studies that combine well-powered GWAS from several diseases and from phenotypes that represent variation in the normal range such as EA are likely to play an important part in this development. However, deep phenotyping of large patient samples will be inevitable to link GWAS results from complex outcomes such as EA and SZ to specific biological disease subgroups.

ONLINE METHODS
All reported statistical results are based on two-sided tests, unless indicated otherwise. Our proxy-phenotype analyses and our replication strategy followed a pre-registered analysis plan (https://osf.io/dnhfk/). The full description of all materials and methods is provided in the Supplementary Note.

GWAS on educational attainment.
The EA sample excluded all cohorts that participated in the GWAS on SZ described below, yielding a sample size of n = 363,502 individuals of European descent 13 . The GRAS replication sample was not part of the GWAS on EA, either.

GWAS on schizophrenia.
The SZ sample consisted of n = 34,409 cases and n = 45,670 controls, diagnosed with SZ or schizoaffective disorder 5 . We excluded the GRAS data collection from the GWAS on SZ.

Proxy-phenotype look-up.
Analyses were carried out using 8,240,280 autosomal SNPs that passed quality controls in both GWAS and additional filters described in the Supplementary Note. We selected approximately independent lead SNPs from the EA GWAS results using the clumping procedure in PLINK 35 . We looked up the SZ results of all approximately independent EA lead-SNPs that passed the pre-defined significance threshold of P EA < 10 −5 .
We tested if the observed sign concordance between EA and SZ is different from 50% using the binomial probability test 36 . "Raw" enrichment factors and "raw" enrichment p-values of the EA lead-SNPs on SZ were calculated by taking the actual distribution of P values in the SZ GWAS result files into account but ignoring the LD scores 12,37 .

LD-aware enrichment across different traits.
We developed an enrichment test that corrects for the LD score of each SNP (Supplementary Note). We conducted this test for the 132 SNPs that are jointly associated with EA and SZ in our proxy-phenotype analyses (P EA < 10 −5 and P SZ < 0.05). LD scores were obtained from the HapMap 3 European reference panel. We investigated SZ and 21 additional traits for which GWAS results were available in the public domain. Some of the traits were chosen because they are phenotypically related to SZ (e.g., BIP), while others were less obviously related to SZ (e.g., age at menarche) or served as negative controls (e.g., fasting insulin). If one of the 132 candidate SNP was not available in the reference panel or the GWAS results of the other traits, we tried to use a good proxy, yielding 79 to 105 available SNPs per trait.

Phenotypic correlations.
We explored the correlations between the number of years of education with 7 quantitative measures of SZ in the GRAS sample of SZ cases: Age at prodrome, age at disease onset, premorbid IQ, GAF, CGI-S, and PANSS positive and PANSS negative scores.

Replication and Bayesian credibility of results.
Our replication uses a PGS in the GRAS data collection, which is based on the 132 independent EA lead-SNPs that are also nominally associated with SZ ( < 10 −5 and < 0.05). This PGS (called SZ_132) was constructed using the regression coefficient estimates of the SZ GWAS as weights. In addition to this polygenic replication strategy, we further probed the credibility of our results using a heuristic Bayesian calculation.

Polygenic prediction of schizophrenia symptoms in the GRAS sample.
We predicted the number of years of education and 7 quantitative measures of SZ in the GRAS sample of SZ cases. For each phenotype, we separately compared the predictive performance of several PGS: Scores constructed from the full GWAS result on SZ, EA, BIP, and neuroticism (called SZ_all, EA_all, BIP_all, Neuro_all, respectively); scores constructed using only the 132 SNPs that are jointly associated with EA and SZ (called EA_132 and SZ_132, using EA and SZ GWAS coefficients as weights, respectively); and two scores that split the SZ_all score into two parts based on sets of SNPs that either have concordant or discordant effects on EA and SZ (called Concordant and Discordant). Genetic outliers of self-reported non-European descent (n = 13 cases) were excluded from the analysis.

Controlling for the genetic overlap between schizophrenia and bipolar disorder.
We estimated GWIS 18 to obtain SNP regression coefficients that are unique to SZ, which are corrected for the genetic overlap between SZ and BIP. The SZ samples used in the GWIS are not overlapping with the samples used in the EA GWAS and they exclude our replication sample (GRAS). BIP GWAS results were obtained from the PGC 38 . We refer to the set of obtained summary statistics as "unique" SZ(min BIP). We then repeated the look-up of the EAassociated lead SNPs in those summary statistics as described above. Similarly, we obtained GWIS results for "unique" BIP(min SZ) using the same method and data. We computed genetic correlations of these GWIS results with EA, childhood intelligence, and neuroticism using bivariate LD score regression 12 and compared the results to those obtained using ordinary SZ and BIP GWAS results.

Simulations of assortative mating.
We conducted simulations to test if strong assortative mating on EA and SZ can induce a spurious genetic overlap between the two traits.

Biological annotation.
To gain first insights into possible biological pathways that are indicated by the genetic loci identified by our PPM analysis, we applied DEPICT 13,39 using a false discovery rate (FDR) threshold of ≤ 0.05. To identify independent biological groupings, we used the Affinity Propagation method based on the Pearson distance matrix for clustering 40 .

Data availability.
GWAS meta-analysis results for EA and SZ as well as GWIS results for unique" SZ(min BIP) and "unique" BIP(min SZ) can be downloaded from the SSGAC website (https://www.thessgac.org/). For information about the GRAS data collection, contact the principal investigator of the study: Hannelore Ehrenreich (ehrenreich@em.mpg.de).

Code availability.
Computer code used to generate LD-aware enrichment and GWIS results can be downloaded from the SSGAC website (https://www.thessgac.org/).   The reported effects are the standardized beta values of linear probability model (LPM) for schizophrenia status. Models 1 -9 differed only in the inclusion of the variables displayed in this table. All models also included the first 10 genetic principal components as control variables. *denotes significance at P < 0.05. ** denotes significance at P < 0.001. The Δ R² is the difference between the adjusted R² of the model compared to the baseline model that only included control variables but no polygenic scores. Linear regression using the first 10 genetic principal components as control variables. 1 : Age of onset was included as covariate. 2 : Medication was included as covariate. *denotes significance at P < 0.05. ** denotes significance after Bonferroni correction (P < 0.05/40 = 0.00125). The Δ R² is the difference between the adjusted R² of the model compared to the baseline model that only included control variables but no polygenic scores.

Figure 1: Workflow of the proxy-phenotype analyses
Notes: Educational attainment (EA) and schizophrenia (SZ) GWAS results are based on the analyses reported in ref. 5,13 . All cohorts that were part of the SZ GWAS were excluded from the meta-analysis on EA. The GRAS data collection was not included in either the SZ or the EA meta-analysis. Proxy-phenotype analyses were conducted using 8,240,280 autosomal SNPs that passed quality control. Genetic outliers of non-European descent (n = 13 cases) were excluded from the analysis in the GRAS data collection (Supplementary Note). Notes: Panel a: The x axis is the chromosomal position, and the y axis is the significance on the −log10 scale. The black dashed line shows the suggestive significance level of 10 −5 that we specified in our preregistered analysis plan. Red and green crosses identify EA-associated lead-SNPs that are also associated with SZ at nominal or Bonferroni-adjusted significance levels, respectively. Panel b: SNPs with concordant effects on both phenotypes are pink, and SNPs with discordant effects are blue. SNPs outside the grey area (21 SNPs) pass the Bonferroni-corrected significance threshold that corrects for the total number of SNPs we tested (P < 0.05/506 = 9.88×10 -5 ) and are labelled with their rs numbers. Observed and expected P values are on the −log10 scale. For the sign concordance test: P = 0.40. Notes: Color codes illustrate the degree of LD-aware enrichment. Red and blue represent stronger and weaker LD-aware enrichment than expected, respectively. Darker colors illustrate more substantial deviation from expectation. A star (*) indicates that the observed LD-aware enrichment is significant after Bonferroni correction for the number of SNPs that were tested for the specific trait (ranging from 79 to 105). Note that the p-values of this test are not the same as those reported in the GWAS and proxy-phenotype analyses because the hypothesis that was test here is different (Supplementary Note). The lines on the left side represent hierarchical clustering using the euclidean distance and the complete agglomeration method. Note that we restricted our analysis to HapMap3 SNPs because these are available in most publically available GWAS summary statistics. If a candidate SNP from our proxy-phenotype results was not included in HapMap3, we replaced it by the best available proxy SNP that was available ( 2 > 0.8 and a maximum distance of 500 kb to our missing EA lead SNPs, choosing the proxy with the highest 2). See Supplementary Note for more details.

Figure 4: Genetic correlations of GWAS and GWIS results that are central to the relationship between SZ and EA.
Notes: The heatmap displays the genetic correlations across 7 sets of GWAS or GWIS summary statistics. Genetic correlations were estimated with LD score regression 12 . The color scale represents the genetic correlations ranging from -1 (red) to 1 (blue). Asterisk denotes that the genetic correlation is significant at P value < 0.01, without adjustment for multiple comparisons.

GWAS on educational attainment
We obtained GWAS summary statistics on educational attainment (EA) from the Social Science Genetic Association Consortium (SSGAC). The results are based on the analyses reported in Okbay et al. 1 , including the UK Biobank. However, our project required nonoverlapping GWAS samples for EA and schizophrenia (SZ). Therefore, all cohorts that were part of the most recent GWAS on SZ by the Psychiatric Genomics Consortium (PGC) 2 (deCODE, DIL, EGCUT, FTC, FVG, H2000, KORA, MGS, WTCCC58C) were excluded from the meta-analysis on EA, yielding a total sample size of n = 363,502. Our replication sample (the Göttingen Research Association for Schizophrenia -GRAS, see Supplementary Note section 6) was not part of the GWAS on EA.

GWAS on schizophrenia
The PGC shared GWAS summary statistics on SZ with us. The results are based on Ripke et al. 2 , but excluded data from our replication sample (GRAS, see Supplementary Note section 6), yielding a total sample size of n = 34,409 cases and n = 45,670 controls.

Quality control
Data sources and quality control procedures for the GWAS on EA and SZ are described in Okbay et al. 1  Before we proceeded with our proxy-phenotype analyses, we applied the following additional quality control steps: 1. To maximise statistical power, we excluded single nucleotide polymorphisms (SNPs) that were missing in large parts of the two samples. Specifically, we continued with SNPs that were available in at least 19 out of 50 cohorts in the SZ results 2,a and in N > 200,000 in the EA meta-analysis 1 . This step excluded 3,778,914 and 6,369,138 genetic markers for EA and SZ, respectively.
2. We dropped SNPs that were not available in both GWAS results files. This step restricted our analyses to the set of available genetic markers that passed the qualitycontrol filters in both the EA and the SZ GWAS results, leaving us with 8,403,560 autosomal SNPs.
3. We dropped 6 SNPs with non-standard alleles (i.e. not A, C, T, or G) and 2 SNPs with mismatched effective alleles. Furthermore, we dropped 163,272 SNPs in the first and the 99 th percentile of the distribution of differences in minor allele frequency (MAF) in the two results files. This final step eliminated SNPs that were a The actual N per SNP was not provided in the SZ GWAS summary statistics.
. CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint likely to be affected by coding errors, strand flips, or substantial differences in MAF in the EA and SZ samples.
The remaining 8,240,280 autosomal SNPs were used in the proxy-phenotype and prediction analyses described in Supplementary Note sections 4, 5, 7 and 8.

Selection of education-associated candidate SNPs
We conducted our proxy-phenotype analyses following a pre-registered analysis plan (https://osf.io/dnhfk/), using 8,240,280 autosomal SNPs that passed quality control (see section 3).

Setting the P value threshold for the proxy-based SNPs
Ideally, proxy-phenotype analyses should use a pre-specified P value threshold that maximises the expected number of true positive results in the look-up stage. The optimal threshold trades off between two opposing effects. On the one hand, a less stringent threshold yields a larger number of candidates that are forwarded to the second stage. A larger set of candidates is more likely to contain true positives. On the other hand, a larger number of candidates requires that a more stringent experiment-wide significance level needs to be applied in the second stage to adjust for multiple testing, which decreases power to pick out the true positives from among the set of candidates 3 . In principal, it is possible to calculate the optimal P value threshold based on the observed distribution of standardised effects sizes and standard errors in the GWAS results of the proxy, the genetic correlation between the two traits, their SNP-based heritabilities, and the sample size of the target phenotype. Given these parameters, one can infer the expected effect size of a genetic variant on the target from the results on the proxy phenotype, calculate the statistical power for the look-up, and approximate the number of expected true positive associations from the look-up at various P value thresholds 4 .
In the current case, the value of such theoretical calculations is limited because the genetic correlation between EA and SZ does not reflect the actual genetic overlap between the two traits adequately. Specifically, Okbay et al. 1 report low, but significant, bivariate LD score regression estimates for the genetic correlation between EA and SZ of 0.08 (P = 3.2×10 -4 ). In addition, the 74 genome-wide significant EA loci have only 51% sign concordance with the SZ results. This sign concordance pattern cannot be differentiated from what would be expected by chance for two traits that exhibit no genetic overlap. Yet, the same 74 EAassociated loci are strongly enriched for association with SZ (enrichment P value < 0.002), which strongly rejects the hypothesis of no genetic overlap between the two traits. This implies that standard formulas 4 to calculate the expected effect size of a genetic variant on SZ from the observed results on EA will be too conservative because they ignore the specific pattern of split sign concordance but strong enrichment for this pair of traits.
Instead of calculating a noisy theoretical optimum, we follow Rietveld et al. 3 and selected 10 -5 as the default P value threshold prior to carrying out the proxy-phenotype analyses (https://osf.io/dnhfk/).

Look-up, sign concordance and enrichment
We used the GWAS results on EA and SZ (see Supplementary Note sections 1-2) that passed our quality control (Supplementary Note section 3) for the analyses described here.

Look-up
To select approximately independent SNPs from the EA GWAS results, we applied the clumping procedure in PLINK version 1.9 5,6 using r 2 > 0.1 and 1,000,000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel 7 to estimate linkage disequilibrium (LD) among SNPs. This algorithm assigns the SNP with the smallest P value as the lead SNP in its "clump". All SNPs in the vicinity of 1,000,000 kb around the lead SNP that are correlated with it at r 2 > 0.1 are assigned to this clump. The next clump is formed around the SNP with the next smallest P value, consisting of SNPs that have not been already assigned to the first clump. This process is iterated until no SNPs remain with P < 10 -5 , leading to 506 approximately independent EA-associated lead SNPs. 108 of the 506 EAassociated lead SNPs are genome-wide significant (P < 5 × 10 -8 ).
We looked up the SZ GWAS results (Supplementary Note section 2) for these 506 EAassociated lead SNPs. Results for all 506 SNPs are reported in Supplementary Table 5.1. Figure 2a shows a Manhattan plot for the GWAS on EA highlighting SNPs that were also significantly associated with SZ (red crosses for PSZ < 0.05, green crosses for PSZ < 0.05/506 = 9.88 × 10 -5 ). Figure 2b presents a Q-Q plot of the look-up. 132 SNPs are associated with SZ at nominal significance (P < 0.05) and 21 of these survive Bonferroni correction (PSZ < 9.88 × 10 -5 ).
In order to investigate the novelty of the findings, we extracted all the SNPs in LD with these 21 SNPs at r 2  0.1 with a maximum distance of 1,000 kb using the 1000 Genomes phase 1 European reference panel. Of these 21, 12 are in LD with loci previously reported by Ripke et al. and 2 are in the major histocompatibility complex (MHC) region on chromosome 6 and were therefore not separately reported in that study 2 . Three of the 21 loci that were not identified yet by Ripke et al. 2 were independently found in a recent meta-analysis of the PGC results with another large-scale cohort, yielding a total sample size of n = 40,675 SZ cases and n = 64,643 controls 8 , strengthening the credibility of these loci. Furthermore, two are in LD with loci recently reported by Hellard et al. 9 who used a false discovery rate (FDR) method that is less conservative than our approach. The results of that latter study were also based on a much smaller GWAS sample for EA (n = 95,427) and partially overlapping samples between EA and SZ. Thus, our finding lends additional credibility to the suggestive associations reported in Hellard et al 9 .

Bayesian credibility of results
We probed the credibility of our proxy-phenotype association results using a heuristic Bayesian calculation following Rietveld 3 . We focus on the 21 EA-associated lead SNPs that are also associated with SZ after Bonferroni correction. Bayes' Rule implies that the probability that an association is true given that we observe significance is given by Power", as well as the significance test, are 2-sided, is the prior belief that the SNP is truly associated, and is the significance threshold used for testing (in our case, = 0.05 506 = 9.88 × 10 −5 ).
To calculate power for each SNP, we computed the winner's curse corrected odds ratio using the procedure described in Rietveld  An important question is which prior beliefs are reasonable starting points for these Bayesian calculations. For an arbitrarily chosen SNP, the most conservative reasonable prior would assume that each truly associated SNP has the same effect size as the strongest effect size that was actually observed in the data. If one divides the SNP-based heritability of the trait by that effect size in R 2 units, one obtains a lower bound for the number of SNPs that can be assumed to be truly associated. To aid this line of thinking, we converted the winner's curse corrected odds ratios of our 21 SNPs into R 2 using where is Cohen's d, which is calculated as and is a correction factor that adjusts for the MAF of the SNP. This correction factor is calculated as where 1 = × and 2 = × (1 − ), see Borenstein et al. 12 The largest effect size in 2 that we observe in our results is rs4378243 with 0.044%. The SNP-based heritability of SZ is ≈21% 13 . Thus, if all causal SZ SNPs would have an effect of 2 = 0.044%, we would expect that ≈500 truly causal loci exist. The chance of finding any one of them by chance from a set of ≈500,000 independent loci in the human genome is . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint ≈0.1%. b However, in reality most truly associated loci for SZ will surely have smaller effects than that. Thus, a prior belief of ≈0.1% is certainly too conservative. Furthermore, the SNPs we investigate are not arbitrary but selected based on their association with another, genetically related cognitive trait (EA) in a very large, independent sample. Thus, a prior belief of 1% or 5% that these SNPs are also associated with SZ is probably more reasonable. As an upper bound, we assume that 10% of all loci are causal. Thus, the chance to pick any one of them by chance would be 10%. Table 1 displays the winner's curse corrected effect size of the 21 EA-associated lead SNPs that are also associated with SZ after Bonferroni correction. It also shows the posterior probability that these SNPs are truly associated with SZ given our results for prior beliefs ranging from 0.1%, 1%, 5%, to 10%. Thirteen of these SNPs have posterior probabilities of being true positives of >50% for even the most conservative prior. For a more realistic prior belief of 5%, all 21 SNPs are likely or almost certain to be true positives.

Sign concordance
We compared the signs of the beta coefficients of the 506 EA lead SNPs ( < 10 −5 ) with the beta coefficients for SZ. If the signs were aligned, we assigned a "1" to the SNP and "0" otherwise. By chance, sign concordance is expected to be 50%. We tested if the observed sign concordance is different from 50% using the binomial probability test 14 . 263 of the 506 SNPs have the same sign (52%, P = 0.40, 2-sided).

Enrichment
Because EA and SZ are highly polygenic, we tested for enrichment by taking the actual distribution of P values in the GWAS result files into account.

Raw enrichment factor (not corrected for LD score of SNPs)
Due to the polygenic architecture of both traits, it is expected to find some EA-associated SNPs that are also associated with SZ just by chance even if both traits are genetically independent. Under this null hypothesis, the expected number of EA-associated lead SNPs that are also significantly associated with SZ is where , is the total number of independent lead SNPs in the EA GWAS results, and and are the shares of SNPs in , that have P values for EA and SZ below a certain threshold, respectively. b It is typically assumed that GWAS data for European populations contain ≈1,000,000 independent loci. However, the quality-control procedures for GWAS summary statistics in studies like ours decreases the number of independent loci to <1,000,000 41,62 . In fact, clumping the post-QC GWAS results for SZ without a P value threshold, an R 2 LD<0.1, and a LD-window of 1,000,000 kb with the 1000 Genomes phase 1 version 3 European reference panel 42 leads to only 223,065 independent loci. Thus, assuming 500,000 independent loci in these calculations is conservative.

Raw enrichment P value (not corrected for LD score of SNPs)
Following Okbay et al. 4 , we performed a non-parametric test of joint enrichment that probes whether the EA lead SNPs are more strongly associated with SZ than randomly chosen sets of SNPs with MAF within one percentage point of the lead SNP. To perform our test, we randomly drew 10 matched SNPs for each of the 506 EA lead SNPs with < 10 −5 .
We then ranked the 506×10 randomly matched SNPs and the original 506 lead EA SNPs by P value and conducted a Mann-Whitney test 15 of the null hypothesis that the P value distribution of the 506 EA lead SNPs are drawn from the same distribution as the 506×10 randomly matched SNPs. We reject the null hypothesis with = 6.872 × 10 −10 ( = 6.169, 2-sided). As a negative control test, we also calculated the raw enrichment P value of the first randomly drawn, MAF-matched set of SNPs against the remaining 9 sets, yielding P = 0.17.
Repeating this raw enrichment test for the subset of 21 EA-associated SNPs that remained significantly associated with SZ after Bonferroni correction (threshold < 0.05 506 = 9.88 × 10 −5 )yields = 5.44 × 10 −14 ( = 7.521, 2-sided). The negative control test based on the raw enrichment P value of the first randomly drawn, MAF-matched set of SNPs against the remaining 9 sets yields P = 0.34.

LD-aware enrichment test
The "raw" enrichment reported in Supplementary Note section 5.4.1 could in principle be due to the LD-structure of the EA lead SNPs that we tested. Specifically, if these EA lead SNPs have stronger LD with other SNPs in the human genome than expected by chance, this could cause the observed enrichment of this set of SNPs on SZ and other traits because higher LD increases the chance these SNPs would "tag" causal SNPs that they are correlated with. To test whether the observed enrichment is due to LD, we developed an LD-aware enrichment test. Furthermore, we used this LD-aware enrichment test to probe if the observed enrichment of our EA lead SNPs can also be observed for other traits.
We investigated the SZ GWAS results described in Supplementary Note section 2 and 21 additional traits for which GWAS summary statistics were available in the public domain (Supplementary Table 5.2). Some of the traits were chosen because they are phenotypically related to SZ (bipolar disorder (BIP), neuroticism, depressive symptoms, major depressive disorder, autism, and childhood intelligence (IQ)), while others were less obviously related to SZ (e.g. intracranial volume, cigarettes per day) or to the brain (e.g. age at menarche, inflammatory bowel disease). Finally, we included five traits as negative controls (height, birth weight, birth length, fasting (pro)insulin).
We calculated the LD score regression intercept and slope of the traits using LDSC 16 . For SNP i in trait j, the expected chi-square statistic can be calculated as where N is the sample size of the target trait j, h 2 is the heritability of trait j, = ∑

=1
for SNP i is calculated using HapMap3 SNPs from European-ancestry, M is the number of SNPs included in the calculation of the LD score (n = 1,173,569 SNPs), 2 is the squared correlation between SNPs j and k in the HapMap3 reference panel, and 1 + Na is the LD score regression intercept for trait j. We used precomputed LD scores available from the LDSC software 16 .
As recommended by Bulik-Sullivan et al. 16 , we restricted our analysis to HapMap3 SNPs (using the --merge-alleles flag) because these seem to be well-imputed in most studies. Out of 132 SNPs with < 1 × 10 −5 and < 0.05, only 30 SNPs are directly present in HapMap3 SNP list. Therefore, we extracted proxy SNPs with 2 > 0.8 and a maximum distance of 500 kb to our missing EA lead SNPs and chose the one with the highest 2 as a proxy. After this step, we could include 105 (out of 132) SNPs in our analyses. For each of these 105 SNPs, we observed the Z-statistics in the publically available GWAS results of the traits. Z-statistics were converted into Chi 2 statistics by squaring them. The LD score corrected enrichment per SNP for each trait is the ratio of the observed to the expected Chi 2 . The results are shown in Figure 3 (and in Supplementary Table 5.2). To test whether a particular realization is significantly larger than expected (and thus the ratio ℎ 2 / ℎ 2 is significantly greater than one), we test each particular observed Chi 2 against a non-central Chi 2 distribution with = 1 + degrees of freedom and a non-centrality parameter of × ℎ 2 × / . The intuition behind using the LD score regression slope as the non-centrality parameter is as follows: True genetic signal, not caused by stratification, covaries with LD 16 . True genetic signal also results in non-central Chi 2 statistics (i.e. ≠ 0 leads to non-centrality in ( ) 2 ). Therefore, under a polygenic model, SNPs with a high LD score are expected to tag several true effects, raising our expectation of their Chi 2 statistic. Here we account for the LD score specific expectation of effect size. We show the selected SNPs influence our traits of interest over and above what is expected based on LD and, importantly, our results suggest that this is not the case for several other phenotypes, revealing a certain amount of specificity in the enriched SNPs.
Furthermore, since the SNPs considered for enrichment are independent, their Chi 2 and noncentrality parameters are additive. This additivity allows us to formulate an expected distribution of the sum of ℎ 2 , based on the sum of non-centrality and the sum of  Our LD-aware enrichment test has two limitations. First, LD score regression assumes that allele frequency (AF) does not correlate with effect size, an assumption which has been empirically shown to be violated for low-frequency alleles 17 . Second, our test assumes the absence of selection on the trait. Variation in AF and the degree of negative selection could explain excess signal in low LD SNPs 18 . However, our raw enrichment P value (see Supplementary Note section 5.4.2) is robust to this because it takes the AF of the candidate SNPs explicitly into account. Figure 3 summarise the results. We find that the enrichment of EA-associated SNPs for SZ cannot be explained by the LD scores of these SNPs: The set of 105 SNPs is jointly associated with SZ after controlling for their LD scores (P < 4 × 10 -14 ) and 15 of these SNPs are individually associated with SZ after Bonferroni correction. We also observe LD score corrected enrichment of these SNPs with several other phenotypes, most noticeably with BIP (joint P = 1.1 × 10 -16 ). Four out of 93 tested SNPs are significantly associated with BIP after Bonferroni correction, including one of the SNPs that our proxyphenotype analysis isolated as a novel candidate locus for SZ (rs9575628, a proxy for rs7336518, see Supplementary Table 5.3). We also observe weaker, but still significant joint LD score corrected enrichment of this set of SNPs for inflammatory bowel disease, neuroticism, age at menarche, and childhood IQ. Note that several brain-phenotypes do not show significant enrichment, including depressive symptoms, major depressive disorder, ADHD, and Alzheimer's disease. This implies that the set of SNPs we are testing exhibits some phenotype-specificity is not simply involved in all brain-related outcomes. Also note that none of the negative control phenotypes we included shows significant LD score corrected enrichment (height, birth weight, birth length, fasting (pro)insulin).

Prediction of future genome-wide significant loci for schizophrenia
As reported above (Supplementary Note section 5.1), three of the SNPs that our proxyphenotype approach identified after Bonferroni correction have been reported as novel, genome-wide significant loci for SZ in an effort 8 that was ongoing parallel to ours. Overall, 50 novel loci for SZ were reported in that study. This provides us with the opportunity to ask if our proxy-phenotype approach using GWAS results from EA was able to predict "future" GWAS findings for SZ.
We are using a simple proportions test for this purpose, which compares the ratio of novel SNPs included in our list of 132 loci that are jointly associated with EA and SZ ( < 10 −5 and < 0.05) with the ratio observed in all remaining approximately independent loci with < 0.05.
To identify LD partners and to clump our GWAS results, we used a threshold of 2 > 0.1 and a 1,000,000 kb window in the 1000 Genomes phase 1 version 3 European reference panel. Our SZ summary statistics contained 51,721 approximately independent SNPs with < 0.05. We identified 21,430 SNPs in LD with the 50 novel SNPs 8 and 54,425 SNPs in LD with the 128 genome-wide significant loci that were previously reported. 2 We removed SNPs in LD with the previously GWAS hits from our analyses because those SNPs could (by definition) not be identified as novel. The remaining set of 51,528 approximately independent . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint SNPs with < 0.05 in our SZ GWAS results contained one proxy for each of the 50 novel SNPs 8 . 110 SNPs with < 0.05 also exhibited < 10 −5 in the independent EA GWAS sample. Of those 110 SNPs, six were identified as novel SZ loci in the most recent GWAS dataset expansion 8 . Using Fisher's exact test, we rejected the null hypothesis that the proportion of novel SNPs is equal in the two sets (P = 4.1 × 10 -8 , 2-sided). Furthermore, as a robustness check, we performed the analysis again by excluding the SNPs with MAF ≤ 0.1 and found similar results (P = 1.2 × 10 -6 ). Thus, we conclude that conditioning GWAS results on SZ with independent GWAS evidence on EA significantly outperforms pure chance in predicting GWAS results on SZ from even larger samples.

6
The GRAS data collection All parts of GRAS data collection comply with the Helsinki Declaration and were approved by the ethical committee of the Georg-August-University of Göttingen (master committee) as well as by the respective local regulatory/ethical committees of all collaborating centres.

GRAS schizophrenia and schizoaffective patients
The GRAS data collection has been established over the last 10

GRAS healthy controls
Healthy controls, who gave written informed consent, were voluntary blood donors, recruited by the Department of Transfusion Medicine of the George-August-University of Göttingen (Germany) according to national guidelines for blood donation. As such, they widely fulfil health criteria, ensured by a broad predonation screening process containing standardised questionnaires, interviews, haemoglobin, blood pressure, pulse, and body temperature determinations 19 . Participation as healthy controls for the GRAS data collection was anonymous, with information restricted to age, gender, blood donor health state and ethnicity. For the present study 1,169 subjects were included, 62.2% were male (n = 727) and 37.8% female (n = 442). The average age was 37.44 ±13.27 years (range 18-69).

Genotyping
Genotyping of the GRAS patients and control sample was done with a semi-custom Axiom myDesign genotyping array (Affymetrix, Santa Clara, CA, USA), based on a CEU (Caucasian residents of European ancestry from Utah, USA) marker backbone including 518,722 SNPs, and a custom marker set including 102,537 SNPs. The array was designed using the Axiom Design centre, applying diverse selection criteria 21 . Genotyping was done by Affymetrix on a GeneTitan platform. Several quality control steps were used (SNP call rate >97%, Fisher's linear discriminant >3.6, heterozygous cluster strength off set > -0.1, and . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint homozygote ratio off set > -0.9). These steps were completed with use of either genotyping console software (Affymetrix) or R. In a subsequent step, markers in X, Y, and mitochondrial chromosomes and those with Hardy-Weinberg equilibrium P < 1 × 10 -6 (GRAS healthy controls) or P < 1 × 10 -10 (GRAS patients) were removed, leaving 589,921 SNPs available for analyses.

Imputation and estimation of genetic principal components
The full details of genotype imputation and the estimation of the genetic principal components in the GRAS data collection are described elsewhere 2 . Briefly, imputation was done with the prephasing or imputation approach implemented in IMPUTE2 and SHAPEIT (chunk size 3 Mb) 22,23 . A version of the phase 1 integrated variant set release (v3) from the full 1000 Genomes Project dataset (March 2012) that is limited to variants with more than one minor allele copy ("macGT1"; Aug 26, 2012) was used as imputation reference dataset (INFO value > 0.1 and MAF >0.005). Ten principal components of the genetic data in the GRAS sample were obtained using the standard PGC protocol 2 .

Phenotyping procedures
Patients from the GRAS data collection were examined by the GRAS team of travelling investigators after giving written informed consent, own and/or authorised legal representatives. The GRAS team of travelling investigators consisted of 1 trained psychiatrist and neurologist, 3 psychologists and 4 medical doctors/last year medical students. All investigators had continuous training and calibration sessions to ensure the highest possible agreement on diagnoses and other judgments as well as a low interrater variability regarding the instruments applied. A full description of the GRAS data collection standard operating procedures is provided elsewhere 20 .
Deep clinical phenotyping data were available for all the GRAS patients. For the purpose of the present study the following phenotypes for SZ were selected: i) age at prodrome which precedes SZ onset and is characterized by cognitive decline, social withdrawal, and depression; ii) age at disease onset defined as onset of first psychotic episode; iii) premorbid IQ using the MWT-B (Mehrfachwahl-Wortschatz-Intelligenztest-B) which estimates the intellectual functioning of a person prior to known or suspected onset of disease 24 ; iv) positive and negative symptoms using the Positive and Negative Syndrome Scale (PANSS) 25 . Each domain ranges from 7 to 49 with higher scores indicating severe symptoms; v) the Global Assessment of Functioning (GAF) which subjectively rates the social, occupational, and psychological functioning of an individual, e.g., how well one is meeting various problems-in-living, on a continuum ranging from 1 to 100 26 . Poorer functioning is indicated by lower GAF scores; vi) the Clinical Global Impression of Severity (CGI-S) scale which rates illness severity from 1 to 7 with higher scores indicating more severe illness. 27 ; vii) years of education was measured based on the highest degree obtained, converted into US-schooling year equivalents based on the 1997 International Standard Classification of Education (ISCED) of the United Nations 28 . Education was assessed retrospectively at the time of patient recruitment for the GRAS data collection. At this time, 75% of the GRAS patients were older than 29 years.
An important feature of SZ patients that may influence their everyday functioning and performance, and result in a considerable number of side effects, is their antipsychotic medication. Medication was assessed as the dose of present antipsychotic medication of each patient at the moment of the interview, expressed as chlorpromazine equivalents 29 .

Replication in the GRAS data collection
We showed in our pre-registered analysis plan that our replication sample (GRAS) is not large enough to replicate individual SNPs (https://osf.io/dnhfk/). Instead, we decided at the outset to attempt replication of the proxy-phenotype analysis results using a polygenic score (PGS) that consists of the >80 most strongly associated, independent SNPs. The set that best meets this criterion are the 132 independent EA lead-SNPs that are also nominally associated with SZ ( < 0.05), see Supplementary Note section 5. PGS for this set of 132 candidate SNPs were constructed using either the coefficient estimates of the EA or the SZ GWAS meta-analysis, resulting in two different scores (named EA_132 and SZ_132 in Supplementary Tables 7.1 In addition, we also constructed PGS for EA, SZ, BIP, and neuroticism in the GRAS data collection using all available SNPs as control variables for multivariate prediction analyses

Polygenic score calculations
PGS were calculated using PLINK version 1.9 5,6 . We calculated eight different scores, which are described below.

Schizophrenia scores
We received the GWAS summary statistics for SZ from the PGC excluding the data from our replication sample (GRAS). We constructed a PGS using the 132 EA lead-SNPs ( < 10 −5 )that are also nominally associated with SZ ( < 0.05). This score (SZ_132) is used for replication of the proxy-phenotype analyses described in Supplementary Note section 5.
Furthermore, we constructed a PGS using all 8,240,280 SNPs that survived quality control (SZ_all, see Supplementary Note section 3). Next, we applied the clumping procedure using r2 > 0.1 and 1,000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 349,357 SNPs ready for profile scoring.
For secondary analyses on the prediction of SZ symptoms (Supplementary Note section 8.3), we constructed two PGS using the 4,147,926 SNPs and 4,092,354 SNPs that have concordant (+ and +; orand -) or discordant signs (+ and -; orand +) for EA and SZ, respectively. Clumping resulted in 260,441 and 261,062 independent SNPs with concordant or discordant, respectively. We used these approximately independent SNPs for profile scoring and call the resulting PGS Concordant and Discordant.
These four scores were calculated using the -score function in PLINK using the natural log of the odds ratio of the SNPs for SZ.   11 for the derivation. Using these betas, we constructed a PGS using the 132 EA lead-SNPs ( < 1 × 10 −5 )that are also nominally associated with SZ ( < 0.05). The resulting score is called EA_132.
Furthermore, we constructed a PGS using all 8,240,280 SNPs that survived quality control (see Supplementary Note section 3). Next, we applied the clumping procedure using r 2 > 0.1 and 1,000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate linkage disequilibrium among SNPs, eventually leaving a set of 348,429 SNPs ready for profile scoring. The resulting score is called EA_all.

Bipolar disorder score
We obtained GWAS summary statistics on BIP from the PGC 30 . We used the LD-pruned GWAS summary from PGC ("pgc.bip.clump.2012-04.txt") with a set of 108,834 LD-pruned SNPs ready for profile scoring. PGS for the GRAS data collection were calculated by the application of the -score function in PLINK using the natural log of the odds ratio. The resulting score is called BIP_all.

Neuroticism score
We obtained GWAS summary statistics on Neuroticism from the SSGAC. The results are based on the analyses reported in Okbay et al. 4 containing 6,524,432 variants. We applied the clumping procedure using r 2 > 0.1 and 1,000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 232,483 SNPs ready for profile scoring. PGS for the GRAS data collection were calculated by the application of -score function in PLINK using the Neuroticism beta values. The resulting score is called Neuro_all. Note that our replication sample (GRAS) was not included in the GWAS summary statistics of any of these traits.

Predicting case-control status using PGS
To investigate if our proxy-phenotype results help to predict SZ case-control status, we estimated a linear probability model (LPM) in Stata 31,32 . Genetic outliers (n = 13) based on self-reported non-European ancestry were excluded from all prediction analyses. Following the PGC protocol described in Ripke et al. 2 , we included 10 principal components as covariates (but no other variables). The proportion of variance explained (Adjusted R 2 ) was computed by the comparison of a full model (covariates + PGS) to a reduced model (covariates only).
Results of the predictions with different models are summarised in Table 2. As effect sizes, we are reporting standardised regression coefficients. SZ_132 scores significantly predict the case-control status in our replication sample ( = 5.4 × 10 −0 ) ( Table 2 Model 1) and remain significant even if we include SZ_all score as a control variable ( = 1.7 × 10 −0 ) ( Table 2 Model 3). The EA_132, EA_all and Neuro_all scores did not predict case-control status in our replication sample. The BIP_all score significantly predicts case-control status (Table 2 Model 7), which is expected given the genetic correlation between BIP and SZ 33,34 . Interestingly, the SZ_132 score still significantly predicts case-control status when all other scores are included as control variables (P = 3.4 × 10 −0 ) ( Table 2 Model 9), highlighting the importance of these 132 SNPs. Furthermore, as a robustness check we excluded the 132 SNPs from the construction of the SZ_all and EA_all scores. The prediction results did not change except a slight decrease of variance explained with SZ_all scores, which is expected given that 132 SNPs were excluded from the construction of the score (Supplementary Table  7.2).
As an additional robustness check, we also ran a logistic regression with standardised PGS to predict SZ case-control status using the same explanatory variables as in Models 1-9 in Table  2. In all models, the P values of the coefficients were very similar to the ones obtained by LPM. In the equivalent of Model 9, we find that increasing the PGS from its mean by one standard deviation increases the probability of being an SZ case in the GRAS sample by 4.3% for the SZ_132 score ( = 3.1 × 10 −0 ), 15.2% for the SZ_all score ( = 3.3 × 10 −0= ), and 7.1% for the BIP_all score ( = 4.1 × 10 −0= ), respectively.

8
Polygenic prediction of schizophrenia symptoms

Phenotypic correlations
Phenotypic correlations between years of education and different phenotypes of SZ in the GRAS data collection (see Supplementary Note section 6) were calculated using Pearson's correlation. Target SZ phenotypes included were age at prodrome, age at disease onset, premorbid IQ, GAF, CGI-S, and positive and negative symptoms of the PANSS.
Results from the phenotypic correlations between the target phenotypes are reported in Supplementary Table 8.1 and Supplementary Figure 1. In summary, years of education was significantly correlated (all P values <0.0003) with all the quantitative traits of interest for the present study measuring the psychopathology and level of functioning of the patients (-0.212 < r < 0.435). The correlation was positive for age at prodrome and age at disease onset, indicating that the earlier the disease started the lower the level of education achieved. Years of education was also positively correlated with premorbid IQ and GAF and negatively . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint correlated with CGI-S and PANSS positive and negative scores, indicating that higher levels of education were associated with less severe disease outcomes. Our results are in line with previous studies suggesting that less educated SZ patients are at higher risk of a poorer course of the disease 35 . Moreover, we found a strong positive correlation between age at prodrome and age at disease onset (r = 0.918; P < 0.0001). Premorbid IQ was positively correlated with GAF (r = 0.200; P < 0.0001) and negatively correlated with CGI-S, PANSS positive and negative (-0.277 < r < -0.110; all P values < 0.0009). Measures of disease severity (CGI-S, PANSS positive, PANSS negative) were positively correlated with each other (0.478 < r < 0.638; all P values < 0.0001) and negatively correlated with an assessment of global functioning (GAF, -.579 < r < -0.824; all P values <0.0001).

Based on the 132 EA lead-SNPs
To investigate if our proxy phenotype results can predict specific SZ features, we predicted eight quantitative outcomes among the SZ cases in the GRAS data collection: years of education, age at prodrome, age at disease onset, premorbid IQ, GAF, CGI-S, PANSS positive and PANSS negative scores. These phenotypes are described in Supplementary Note section 6.
SZ and BIP, although being classified as two different disorders, share several symptoms. The PANSS is a clinical instrument principally developed for use in SZ to identify the presence and severity of psychopathology symptoms. However, BIP patients also present some symptoms measured by this scale manifesting the phenotypic overlap between both diseases 36,37 . Thus, we used the sum score of the positive symptoms (delusions, conceptual disorganization, hallucinations, hyperactivity, grandiosity, suspiciousness/persecution, hostility) and the sum score of the negative symptoms (blunted affect, emotional withdrawal, poor rapport, passive/apathetic social withdrawal, difficulty in abstract thinking, lack of spontaneity and flow of conversation, stereotyped thinking) included in PANSS to test if the genetic associations we identified using our proxy-phenotype approach predict symptoms that SZ and BIP share.
For the prediction of each phenotype a linear regression model was used including the PGS described in Supplementary Note section 7 (SZ_132, SZ_all, EA_132, EA_all, BIP_all and Neuro_all). Each regression included 10 principal components as covariates. The regressions for the prediction of years of education and premorbid IQ also included the age of onset as a covariate. Furthermore, the regressions for the prediction of GAF, CGI-S and the PANSS scores also controlled for medication because medication significantly affects these outcomes. We calculated the marginal R 2 of each PGS by squaring its standardised beta coefficient. Results of the predictions have been summarised in Supplementary Table 8.2.a. We found that both the EA_all (stand. beta = 0.17, = 8 × 10 −0 ) and the EA_132 score (stand. beta = 0.08, = 9 × 10 −0 ) were associated with years of education in the GRAS sample of SZ patients. However, the EA_132 score did not survive correction for multiple testing (8 phenotypes and 6 PGS = 48 tests; thus the Bonferroni-adjusted P value is 0.05 48 = 1.04 × 10 −0 . However, we note that neither the phenotypes nor the PGS are strictly independent. Therefore, the Bonferroni correction is likely to be too conservative to obtain a family-wide error rate of 0.05). Since age of disease onset, disease severity and progress of the disease can have an effect on the level of education achieved by a SZ patient, the years of education assessed in the GRAS data collection is not as solid as the measure used in the most recent EA GWAS 1 . The potential measurement error of EA in the GRAS data collection may therefore lead to an underestimation of the predicted association. The EA_all score was also associated with premorbid IQ (stand. beta = 0.14, = 1.8 × 10 −0 ). None of the PGS we constructed were associated with any of the other phenotypes tested (age at prodrome, age at disease onset, GAF, CGI-S, PANSS positive and negative scores) at < 0.03.
As a robustness check, we excluded the 132 EA lead-SNPs from the proxy-phenotype analyses from the construction of the EA_all and SZ_all scores to correct for double-counting of these SNPs [EA_all (without 132) and SZ_all (without 132)]. All prediction models have been analysed again with these PGSs, yielding virtually identical results as our main model specification (Supplementary Table 8

Based on the sign concordance between EA and SZ GWAS results
If heterogeneity in the genetic architecture of SZ subtypes is causing the observed enrichment of EA-associated loci with SZ, the sign concordance pattern of SNPs with both traits may contain relevant information that is pertinent to specific SZ symptoms. We tested this by constructing PGS that take the sign concordance of SNPs for both traits into account. Specifically, we took SNPs and SZ GWAS results that were used to construct the SZ_all score and sorted the SNPs into two sets based on their sign concordance with EA. The resulting two sets of SNPs were used to construct the Concordant and Discordant scores (for details, see Supplementary Note section 7.1.1). A linear regression model was used for the prediction of each phenotype including the PGS described in Supplementary Note section 7 (Concordant, Discordant, EA_all, BIP_all and Neuro_all). Due to the high correlation between the SZ_all score with the Concordant (r = 0.858, P < 0.0001) and Discordant (r = 0.873, P < 0.0001) scores in the GRAS sample of patients (Supplementary Table 7.1.b), we excluded the SZ_all from these prediction models to avoid multicollinearity. Each regression included the first 10 genetic principal components as covariates. The regressions for the prediction of years of education and premorbid IQ also included the age of onset as a covariate. Furthermore, the regressions for the prediction of GAF, CGI-S and the PANSS scores also controlled for medication because medication significantly affects these outcomes.
Results of the predictions are summarised in Table 3. As expected, the EA_all score was associated with years of education and premorbid IQ accounting for 4.2% ( = 2.6 × 10 −0 ) and 3% ( = 2.3 × 10 −0 ) of the variance, respectively. Interestingly, with this new model we could also predict SZ symptoms such as GAF and PANSS (Table 3). For example, all five PGS jointly achieve an ∆ 2 = 1.2% for PANSS negative and ∆ 2 = 1.4% for GAF. Conditional on the effects of the Concordant and Discordant scores, the EA_all score is now associated with less severe disease outcomes, consistent with the observed phenotypic correlations. And conditional on the EA_all score, the Concordant score is now associated with more severe positive and negative symptoms as measured by the PANSS scale, worse global functioning measured by GAF, and higher illness severity measured by the CGI-S. Although a Bonferroni correction for multiple testing is too conservative for our models given that PGS and SZ phenotypes are not independent, 5 of the PGS survive Bonferroni correction ( = 0.05/(5 × 8) = 1.25 × 10 −0 , see Table 3).
Results from linear regression models excluding the EA_all scores (Supplementary Table 8.3) suggest that the association between EA_all and disease outcomes is conditional on the effects of Concordant and Discordant scores. In the same way, the associations between Concordant and Discordant and disease outcomes are conditional on the effects of EA_all scores.
Although the correlation between EA_all and the Concordant (r = 0.256) and Discordant (r = -0.377) scores was only moderate (Supplementary Table 7.1.b), we checked for possible multicollinearity in all our prediction models. The variance inflation factor was less than 3 for all the models, suggesting that multicollinearity is not a concern for the prediction results reported in Table 3 38 .
Since the GRAS data collection includes SZ and schizoaffective disorder (SD) patients, we repeated analyses described above excluding patients that were diagnosed with SD (n = 198). We found that the 95% confidence intervals of the estimated regression coefficients of both model specifications overlapped in all cases, implying that the genetic heterogeneity in SZ that we identify is not only due to SD (Supplementary Table 8.4.a and b).

9
Controlling for genetic overlap between schizophrenia and bipolar disorder

GWIS schizophreniabipolar disorder
One possible reason for the observed genetic overlap between EA and SZ is that both phenotypes could be jointly genetically correlated with other outcomes. In fact, Nieuwboer et al. 39 suggest that the genetic correlation between EA and SZ is probably induced by the genetic correlation between SZ and BIP as well as the genetic correlation between BIP and EA. If that is indeed the case, the EA-associated lead SNPs should not show enrichment for association in "unique" SZ GWAS results that are "purged" of the genetic correlation between SZ and BIP.
To test this hypothesis, we estimated genome-wide inferred statistics (GWIS) 39 to obtain SNP regression coefficients that are unique to SZ, corrected for BIP. We refer to this set of summary statistics as "unique" SZ(min BIP). We then repeated the look-up of the EA-associated lead SNPs in those summary statistics and note that the EA-and "unique" SZ(min BIP) results have been derived from independent samples, similar to our main look-up described in Supplementary Note sections 1-3.
A GWIS infers genome-wide summary statistics for a (non-linear) function of phenotypes for which GWAS summary statistics are available 39 . Here, in particular, we wish to infer for each SNP the effect on SZ, conditioned upon its effect on BIP. One possible approximation involves a GWIS of the following linear regression function: where the parameter is estimated from the genetic covariance between SZ and BIP and the genetic variance in BIP as = . The residual ( ) is actually our trait of intrest, for which we use the term SZ(min BIP). Using GWIS we infer the genome wide summary statistics for SZ(min BIP) given the most recent PGC GWAS results for SZ (omitting the GRAS data collection) 2 and BIP 40 . The effect size with respect to SZ(min BIP) for a single SNP is computed as: The standard error for each SNP effect is approximated using the delta method and accounts for the possible effect of sample overlap between the SZ and BIP GWAS.
As data input, we used the GWAS results on schizophrenia (excluding the GRAS data collection) described in Supplementary Note section 2. GWAS results for BIP 40 were obtained from the website of the PGC (https://www.med.unc.edu/pgc/files/resultfiles/pgc.cross.bip.zip).

Look-up in GWIS results schizophreniabipolar disorder
The "unique" SZ(min BIP) results obtained from the GWIS were processed and merged with the EA GWAS results using the same procedures described in Supplementary Note section 3, leading to 1,153,214 SNPs that passed the quality control thresholds. This number is substantially lower than in the main look-up reported in Supplementary Note section 3 because the BIP GWAS was based on HapMap 2 imputation 41 not on 1000 Genomes imputation 42 like the SZ 2 and EA GWAS 1 .
We repeated the clumping and look-up of the EA-associated lead SNPs in the cleaned and merged "unique" SZ(min BIP) results following the steps described in Supplementary Table 9.1 reports the full results.

Sign concordance
We compared the signs of the beta coefficients of the 346 EA lead SNPs ( < 10 −5 ) with the beta coefficients of the "unique" SZ(min BIP) results. If the signs were aligned, we assigned a "1" to the SNP and "0" otherwise. By chance, sign concordance is expected to be 50%. We tested if the observed sign concordance is different from 50% using the binomial probability test 14 . 154 of the 346 SNPs had the same sign (44.5%, P = 0.046, 2-sided). This result is consistent with a negative genetic correlation between the most strongly EA-associated SNPs and "unique" SZ(min BIP) and contrasts with the positive genetic correlation between EA and SZ reported in Obkay et al 1 .

Raw enrichment factor (not corrected for LD score of SNPs)
We calculated the raw enrichment factor of the EA-associated SNPs in the "unique" SZ(min BIP) results using the approach described in Supplementary  Thus, the observed enrichment of the EA-associated SNPs in the "unique" SZ(min BIP) results is weaker than in our main look-up reported in Supplementary Note section 5.4.1, but it still deviates from the expectation under the null hypothesis.

Raw enrichment P value (not corrected for LD score of SNPs)
We repeated the procedures described in Supplementary Note section 5.4.2 and found a raw enrichment P value of 0.02 ( = 2.317, 2-sided). Thus, although the enrichment of the EAassociated top SNPs is unlikely to be drawn from the same distribution as all "unique" SZ(min BIP) results with the same MAF distribution, the enrichment is weaker than in the main SZ GWAS results that did not control for the genetic overlap between SZ and BIP.

GWIS bipolar disorder -schizophrenia
Using the GWIS method and the data sources described above (Supplementary Note section 9.1), we also "purged" the genetic association results for BIP of their overlap with SZ, obtaining "unique" BIP(min SZ) results.

Genetic correlations of GWAS and GWIS results
To test if the genetic overlap of SZ with EA is partially due to their genetic correlation with other traits, we computed genetic correlations of SZ and EA with three other phenotypes of particular relevance-BIP, childhood IQ, and neuroticism. Given that SZ is sometimes referred to as a cognitive disorder 43,44 , it is somewhat puzzling that previous studies did not find a significant (negative) genetic correlation between SZ and childhood IQ 33 . Furthermore, the personality trait neuroticism has been demonstrated to correlate across various psychiatric disorders and is positively associated ( ≈ 0.4) with a general psychopathology factor (p) 45 . In addition, moderate and strong negative genetic correlations of neuroticism have been reported for EA 1 and depressive symptoms 4 , respectively, raising the possibility that neuroticism may contribute to the genetic overlap between EA and SZ. Finally, Nieuwboer et al. 39 suggest that the genetic correlation between EA and SZ may be induced by the genetic correlation between SZ and BIP as well as the genetic correlation between BIP and EA. Extending this logic, it could also be that the lack of a clear negative genetic correlation between SZ and childhood IQ is induced by the genetic correlation between SZ and BIP.
Our analyses are facilitated by the fact that large-sample GWAS summary statistics are available in the public domain for all five traits (see data sources in Supplementary Table  9.2). In addition to analysing GWAS summary statistics, we also included the GWIS results described above in Supplementary Note sections 9.1 and 9.3 (GWIS SZ(min BIP) and GWIS BIP(min SZ)). Figure 4 display the results. When using the GWAS summary statistics, we obtain results very closely resembling those in earlier studies. In particular, we find a strong positive genetic correlation between SZ and BIP (rg = 0.72, P = 8.57 × 10 -60 ) 33,39 , a positive genetic correlations of EA with childhood IQ (rg = 0.74, P = 2.22 × 10 -30 ) and BIP (rg = 0.27, P = 1.75 × 10 -11 ), as well as a negative genetic correlation between EA and neuroticism (rg = -0.25, P = 7.08 × 10 -17 ) 1 . Also in line with previous reports is the insignificant genetic correlation between SZ and childhood IQ (rg = -0.03, P = 0.61) 33 as well as the positive genetic correlation between SZ and EA (rg = 0.09, P = 1.04 × 10 -04 ) 1 , both of which are counter-intuitive results. Furthermore, we find some positive genetic correlation of neuroticism with SZ (rg = 0.19, P = 4.5 × 10 -07 ) and BIP (rg = 0.10, P = 0.06). However, these results are too weak to justify a GWIS approach that would purge the SZ and BIP results of their genetic correlation with neuroticism.

Supplementary Table 9.2 and
Interestingly, the genetic correlations change substantially when we purge the SZ results of their genetic overlap with BIP (GWIS SZ(min BIP)) and vice versa (GWIS BIP(min SZ)). The . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint genetic correlations between EA and IQ with SZ(min BIP) are now negative and significant (rg = -0.16, P = 3.88×10 -04 and rg = -0.31, P = 6.00 × 10 -03 respectively), which is more in line with the idea of SZ being a cognitive disorder. 43 Furthermore, the genetic correlations of EA and IQ with BIP(min SZ) remain positive and get somewhat stronger (rg = 0.31, P = 2.87 × 10 -07 and rg = 0.33, P = 3.18 × 10 -02 respectively). The genetic correlations of SZ(min BIP) and BIP(min SZ) with neuroticism, however, remain quite stable.
These results suggest two things: First, the genetic overlap between SZ and EA reported in Supplementary Note section 5 and the small, positive genetic correlation between the two traits is indeed to some extent caused by their genetic overlap with both BIP and childhood IQ (and not by their overlap with neuroticism). Second, once we purge the genetic association with SZ of their overlap with BIP, we see that the remaining part of "unique" SZ(min BIP) has negative genetic correlations with childhood IQ and EA. Thus, SZ diagnoses that have been used in large-scale GWAS analyses until now seem to comprise of at least two subtypes of the disease that have different genetic components: One part resembles a cognitive disorder which does not overlap with BIP, and one part does overlap with BIP but is not characterised by cognitive deficiencies.

Simulating assortative mating
Previous studies found substantial assortative mating for EA, with spousal correlations in the range between 0.45 and 0.66, which led to a genetic resemblance among spouses for EAassociated genetic markers 46 . There is also evidence for substantial assortative mating for psychiatric disorders (in particular, SZ) with spousal correlations around 0.4 47 . One possibility is that strong, simultaneous assortative mating on EA and SZ may cause an enrichment of EA-associated loci for SZ. If this happens in the absence of phenotypic and genetic correlations between the two traits, one may call such an enrichment "spurious" because the genetic variants for EA would have no actual influence on SZ, even if they would be found to be robustly associated with SZ. We ran simulations to test how likely it is that our results may be driven by strong, assortative mating.

Description
We simulated an initial generation of n = 25,000. For each individual, two sets of 5,000 genetic markers were drawn from a binomial distribution with a 50% chance to be either 0 or 1. The sum of the two copies of each marker is the genotype of the individual. We assumed that two non-overlapping subsets of markers (500 each) were causal for either EA or SZ. We further assumed that both EA and SZ are 100% heritable, binary outcomes. The frequency of high educational attainment was set to 0.3, and the frequency of SZ to 0.2. The phenotype of each individual was determined by examining the value of a polygenic score based on the known causal markers for each trait. We determined the relevant cut-off point of the score such that it matched the assumed frequencies of both traits in the population.
The initial generation went through a matching algorithm where each person was matched to a spouse and spousal correlations for both traits was assumed to be 0.6. Next, each couple had exactly two offspring and each offspring's genotype was drawn from the genotypes of the parents. The offspring generation was then also matched to a spouse and the process was . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint repeated for 50 generations. After each generation, we tested the causal EA markers for association with SZ and computed the raw enrichment test described in Supplementary Note section 5.2.
Our assumptions were not necessarily meant to represent the real world. Instead, our aim was to simulate a set-up with relatively extreme assumptions about heritability, assortative mating, and the prevalence of SZ, all of which would increase our chance of finding spurious enrichment. However, we note that our model was also deliberately simple: Implicitly, we assumed that there is no selective pressure on any genetic marker and that the effect sizes of the markers are constant over time. In reality, this may not be the case and both the phenotype and its genetic architecture may evolve.

Power
The power of our enrichment test is a function of the power to detect associations of individual genetic markers with SZ. To calculate power, we assumed an odds ratio of 1.117. This odds ratio was calculated using the fact that all causal markers had equal effect sizes and individual markers were drawn from a Bernoulli distribution, such that the polygenic score can be seen as a draw from a binomial distribution with parameters = 1,000 and = 0.5. As described above, this score was used to determine the SZ type. Specifically, a simulated individual had SZ if the score was larger than 513, which gives an overall chance of 20 percent to get SZ. The odds ratio was calculated using the following formula: Where 1 is the chance of getting SZ given that you have at least one of the causal variants increasing the chance of SZ, , and 0 is the chance of getting SZ given that you have the opposite variant. 1 and 0 can calculated from the binomial distribution. 1 is the chance that a draw from a binomial distribution, with parameters = 999 and = 0.5, is larger than 512. 0 is the chance that a draw from the same distribution is larger than 513.
We had 93.3% power to detect the causal markers associated with SZ at nominal significance and 18.1% after Bonferroni correction ( = 0.05 5000 = 10 −5 ). Furthermore, we had 80% power to detect (spurious) effects of ≥ 1.093 at = 0.05 and of ≥ 1.182 after Bonferroni correction. The raw enrichment test had even more power because it took the entire distribution of tested P values into account, not only the effects of one single SNP. Thus, we were well-powered to detect spurious enrichment in our simulation, even if the effects of individual genetic markers are relatively small.

Results
If assortative mating would cause the genetic overlap between EA and SZ, one would expect the absolute value of the Z-statistic from the enrichment test to increase over time. At some point, the difference between the two subsets would become statistically significant. The Zstatistic of the test for each generation is plotted in Supplementary Figure 3. There is no persistent trend in the Z-statistics over time and the mean of the Z-statistics is not significantly different from 0 (P = 0.776). Furthermore, we cannot reject the hypothesis that the Z-statistics are simply drawn from a standard normal distribution (Shapiro-Wilk test, = 0.075, one-sided). Nevertheless, the Z-statistic drops below -1.96 in two simulated . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint generations, indicating that the causal EA markers have lower P values than the non-causal markers. Yet, this result does not persist over time. In conclusion, it may be possible to find (spurious) enrichment in some generations by chance, but it is unlikely that assortative mating is a major cause for the strong genetic overlap that we observed between EA and SZ. 11 Biological annotation 11.1 Prioritisation of genes, pathways, and tissues/cell types with DEPICT To gain first insights into possible biological pathways that are indicated by our genetic associations, we applied Data-driven Expression Prioritized Integration for Complex Traits (DEPICT). DEPICT is a novel data-driven integrative method that uses reconstituted gene sets based on massive numbers of experiments measuring gene expression to (1) prioritise genes and gene sets and (2) identify tissues and cell types were prioritised genes are highly expressed. The method has been described in detail elsewhere 1,4 . The input for our analyses (DEPICT version 1 release 194) were the 132 EA lead-SNPs that are also nominally associated with SZ.

Significant reconstituted gene sets
DEPICT identified 111 significant reconstituted gene sets at an FDR below 5% (Supplementary Table 10.1). To identify independent biological groupings, we computed the pairwise Pearson correlations of all significant gene sets using the "network_plot.py" script provided with DEPICT. Next, we used the Affinity Propagation method on the Pearson distance matrix to cluster the findings 48 . The Affinity Propagation method automatically chooses an exemplar for each cluster (Supplementary Table 10.2). npBAF complex (7-set cluster) is named after the GO cellular component (GO:0071564) defined as "A SWI/SNF-type complex that is found in neural stem or progenitor cells". The prefix np stands for neural progenitor. The npBAF complex is essential for the selfrenewal/proliferative capacity of multipotent neural stem cells.
Transcription cofactor activity (12-set cluster) is named after the GO molecular function (GO:0003712) defined as "Interacting selectively and non-covalently with a regulatory transcription factor and also with the basal transcription machinery in order to modulate transcription. Cofactors generally do not bind the template nucleic acid, but rather mediate protein-protein interactions between regulatory transcription factors and the basal transcription machinery".
REACTOME_TRANSMISSION_ACROSS_CHEMICAL_SYNAPSES (15-set cluster) is named after the Reactome pathway centred on genes involved in transmission across chemical synapses. Chemical synapses are specialised junctions that are used for communication between neurones, neurones and muscle or gland cells. The pre-synaptic neurone communicates via the release of neurotransmitter which binds the receptors on the postsynaptic cell.
Dendrite (21-set cluster) is named after a large and heterogeneous GO cellular component (GO: 0030425) defined as "A neuron projection that has a short, tapering, often branched, . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint morphology, receives and integrates signals from other neurons or from sensory stimuli, and conducts a nerve impulse towards the axon or the cell body. In most neurones, the impulse is conveyed from dendrites to axon via the cell body, but in some types of unipolar neurone, the impulse does not travel via the cell body". Abnormal cerebral cortex morphology (5-set cluster) is named after the Mammalian Phenotype category (MP:0000788) defined as "any structural anomaly of a thin layer of grey matter on the surface of the cerebral hemisphere that folds into gyri".
Dendritic spine morphogenesis (3-set cluster) is named after the GO biological process (GO:0060997) defined as "The process in which the anatomical structures of a dendritic spine are generated and organised. A dendritic spine is a protrusion from a dendrite and a specialised subcellular compartment involved in synaptic transmission". REACTOME_AXON_GUIDANCE (4-set cluster) is named after the Reactome pathway centred on genes involved in axon guidance. Axon guidance or axon pathfinding is the process by which neurones send out axons to reach the correct targets.
GRIN2A PPI subnetwork (9-set cluster) is named after the gene GRIN2A (glutamate ionotropic receptor NMDA-type subunit 2A). This gene encodes a member of the glutamategated ion channel protein family. The encoded protein is an N-methyl-D-aspartate (NMDA) receptor subunit. The most significantly enriched member set is "protein serine/threonine phosphatase complex", which is named after GO cellular component (GO:0008287) defined as "A complex, normally consisting of a catalytic and a regulatory subunit, which catalyzes the removal of a phosphate group from a serine or threonine residue of a protein". SNW1 PPI subnetwork (3-set cluster) is named after the gene SNW1 (SNW Domain Containing 1) encodes a coactivator that enhances transcription from some Pol II promoters.
DLGAP3 PPI subnetwork (3-set cluster) is named after the gene DLGAP3 (Discs Large Homolog Associated Protein 3) that plays a role in the molecular organisation of synapses and neuronal cell signalling.
Neurone recognition (3-set cluster) is named after the GO biological process (GO:0008038) defined as "The process in which a neuronal cell in a multicellular organism interprets its surroundings".
WRB PPI subnetwork (3-set cluster) is named after the gene WRB (Tryptophan Rich Basic Protein) also known as CHD5 (Congenital Heart Disease 5) that has a potential role in the pathogenesis of Down syndrome congenital heart disease.
Site of polarised growth (3-set cluster) is named after the GO cellular component (GO:0030427) defined as "Any part of a cell where non-isotropic growth takes place".
Central nervous system neurone axonogenesis (3-set cluster) is named after the GO biological process (GO:0021955) defined as "Generation of a long process from a neurone whose cell body resides in the central nervous system. The process carries efferent (outgoing) action potentials from the cell body towards target cells".
Regulation of neurone projection development (6-set cluster) is named after the GO biological process (GO:0010975) defined as "Any process that modulates the rate, frequency or extent of neurone projection development. Neurone projection development is the process whose specific outcome is the progression of a neurone projection over time, from its formation to the mature structure. A neurone projection is any process extending from a neural cell, such as axons or dendrites (collectively called neurites)".
. CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint WDR37 PPI subnetwork (7-set cluster) is named after the gene WDR37 (WD Repeat Domain 37) that encodes a member of the WD repeat protein family and may facilitate the formation of heterotrimeric or multiprotein complexes.
REACTOME_PI3K_EVENTS_IN_ERBB4_SIGNALING (2-set cluster) is named after the Reactome pathway centred on genes involved in PI3K events in ERBB4 signalling. ERBB4 is a member of the Tyr protein kinase family and the epidermal growth factor receptor subfamily. It is for the normal development of the embryonic central nervous system, especially for normal neural crest cell migration and normal axon guidance.
Regulation of skeletal muscle fibre development (1-set cluster) is named after the GO biological process (GO:0048742) defined as "Any process that modulates the frequency, rate or extent of skeletal muscle fibre development".
HMGB2 PPI subnetwork (1-set cluster) contains only one single gene set. It is named after the gene HMGB2 (High Mobility Group Box 2) that encodes a member of the non-histone chromosomal high mobility group protein family.

Significant tissue/cell types
DEPICT determines the enrichment of expression in particular tissues and cell types by testing whether the genes overlapping the GWAS loci are highly expressed in any of 209 Medical Subject Heading (MeSH) annotations. Interestingly, all the significantly enriched tissues (FDR < 0.05) are related to the nervous system except retina, which is annotated to sense organs (Fig. 5.b). Furthermore, we observed only 1 significantly enriched cell-type, namely "Neural Stem Cells". All the significantly enriched tissues and cell-type along with the gene names are listed in Supplementary Table 10.3.

Significant gene prioritisation
Any particular locus centred on a SNP may contain multiple genes. One straightforward approach is to nominate a gene that is closest to the SNP. But this approach does not consider if the expression of the gene is likely to be altered or regulated by the causal site in the locus. Therefore, we used DEPICT to map genes to associated loci, which prioritise important genes that share similar annotations in bioinformatic databases. For our 132 lead SNPs, DEPICT significantly prioritized (FDR < 0.05) 56 genes (Supplementary Table 10.4). For the two novel SNPs reported in this study (rs7336518 and rs7522116), DEPICT points to the FOXO6 (Forkhead Box O6) and the SLITRK1 (SLIT and NTRK Like Family Member 1) genes. FOXO6 is predominantly expressed in the hippocampus and has been suggested to be involved in memory consolidation, emotion and synaptic function 53,54 . Similarly, SLITRK1 is also highly expressed in the brain 55 , particularly in the frontal lobe, and has previously been suggested as a candidate gene for neuropsychiatric disorders 56 . In particular, SLITRK1 is also associated with Tourette syndrome, which is characterised by persistent involuntary vocal and motor tics and often occurs together with Obsessive-Compulsive disorder and ADHD 57-59 .

GWAS catalog lookup
In order to investigate the novelty of the 21 SNP associations that were found significant for SZ after Bonferroni correction, reported in Table 1, we performed a lookup in the GWAS catalogue with the SNPs and all their "LD partners" (i.e. all SNPs with an r 2 > 0.5 within a 250kb window). The LD partners were extracted with PLINK 5 using a version of the 1000G . CC-BY-NC-ND 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted March 8, 2017. . https://doi.org/10.1101/114405 doi: bioRxiv preprint reference panel specifically harmonised to combine 1000G phase 1 and phase 3 imputed data 60 , and the reference panel has been described previously 4 . The result of the GWAS catalogue lookup is reported in Supplementary Table 10.5. We searched the GWAS catalog 61 (revision 2016-08-25, downloaded on 2016-08-29) c to see if any of the associated SNPs or their LD partners have been reported to be associated with a phenotype previously. 13 out of the 21 SNPs (or their LD partners) have reported associations with SZ in the GWAS catalogueeight out of the 21 SNPs have no reported associations with SZ. Combining the associations reported in the GWAS catalog and those reported in Ripke et al. 2 , Pardinal et al. 8 , and Hellard et al. 9 , we find two SNPs that have not previously been found to be associated with SZ at genome-wide significance (P < 5×10 -8 ) in any of these sources -rs7522116 and rs7336518. c URL: https://www.ebi.ac.uk/gwas/api/search/downloads/full  Notes: SNPs with concordant effects on both phenotypes are pink, and SNPs with discordant effects are blue. SNPs outside the grey area would have passed the Bonferroni-corrected significance threshold that corrects for the total number of SNPs we tested (P < 0.05/346 = 1.445×10 -4 ). Observed and expected P values are on the −log10 scale. For the sign concordance test: P = 0.046, 2-sided.