Uncovering the Hidden Risk Architecture of the Schizophrenias: Confirmation in Three Independent Genome-Wide Association Studies

Complex diseases, such as schizophrenia, may be influenced by hundreds or thousands of genetic variants that interact with one another in complex ways, and consequently display a multifaceted genetic architecture (1). The genetic architecture of heritable diseases refers to the number, frequency, and effect sizes of genetic risk alleles and the way they are organized into genotypic networks (2). In complex disorders, the same genotypic networks may lead to different clinical outcomes (a concept known as multifinality, which is called pleiotropy in genetics), and different genotypic networks may lead to the same clinical outcome (equifinality, which is also described as heterogeneity) (1, 3). In general, geneticists must expect the likelihood that many genes affect each trait and each gene affects many traits (4). Consequently, research on complex heritable disorders like schizophrenia is likely to yield weak and inconsistent results unless the complexity of their genetic and phenotypic architecture is taken into account (5).

For example, twin and family studies of schizophrenia consistently indicate that the variability in risk of disease is highly heritable (81%) (6, 7), but only 25% of the variability has been explained by specific genetic variants identified in genome-wide association studies (GWAS) (8). This is not surprising for complex disorders like schizophrenia because current GWAS methods have been unable to characterize the gene-gene interactions (Figure 1A) that influence the developing clinical profiles (Figure 1B) in complex ways (10). The frequent failure to account for most of the heritability of complex disorders has been called the “missing” (11) or “hidden” (12) heritability problem.

In past studies of schizophrenia, the missing heritability problem has been approached by analyzing the explained variance in large individual samples or by using meta-analysis to combine data sets (9, 13, 14). Efforts have also been made to consider the impact of variation related to ethnicity, sex, chromosomes, functional observations, or allele frequency (8). Nevertheless, most of the heritability of schizophrenia remains unexplained (8).

We have chosen to measure and characterize the complexity of both the genotypic and the phenotypic architecture of schizophrenia (Figure 1C). Past studies have generally ignored variation in clinical features, categorizing people as either having or not having schizophrenia, and they have looked only at the average effects of genetic variants, ignoring their organization into interactive genotypic networks. We postulate that schizophrenia heritability is not missing but is distributed into different networks of interacting genes that influence different people (15–17). Unlike previous studies that neglected clinical heterogeneity among subjects with schizophrenia (14, 18, 19), we characterized the clinical phenotype in detail. We also allowed for possible developmental complexity, including equifinality (or heterogeneity) and multifinality (or pleiotropy).

We investigated the architecture of schizophrenia in the Molecular Genetics of Schizophrenia (MGS) study, in which all subjects had consistent and detailed genotypic and phenotypic assessments (9). We then replicated the results in two other independent samples in which comparable genotypic and phenotypic features were available: the Clinical Antipsychotic Trial of Intervention Effectiveness (CATIE) and the Portuguese Island studies from the Psychiatric Genomics Consortium (PGC) (19–23).

Method

We first identified sets of interacting single-nucleotide polymorphisms (SNPs) that cluster within subgroups of individuals (SNP sets) regardless of clinical status in the MGS Consortium study, employing our generalized factorization method (24–27) combined with non-negative matrix factorization to identify candidates for functional clusters (17) (see Figures S1 and S2 in the data supplement that accompanies the online edition of this article). This approach performs an unsupervised co-clustering of subjects together with distinguishing genotypic/phenotypic features based on the empirical data alone. We combined the Genetic Association Information Network (GAIN) and non-GAIN samples of the MGS study, which constitute one GWAS (9). The 4,196 cases and 3,827 controls in the MGS study were combined to identify SNP sets. We had data of good quality on 696,788 SNPs on these cases and controls, and from these we preselected 2,891 SNPs that had at least a loose association (p values <1.0×10⁻²) with a global phenotype of schizophrenia (see the data supplement). SNP sets were labeled by a pair of numbers based on the order in which they were chosen by the algorithm (see the data supplement). Each SNP set was composed of a particular group of subjects described by a particular set of homozygotic and/or heterozygotic alleles; subjects and/or SNPs may be present in more than one set (17, 24, 25). The SNP sets identified by our generalized factorization method are optimal clusters of SNPs in particular subjects that encode AND/OR interactions between SNPs and subjects (Figure 2A–F, Table 1; see also Figure S3 and the Method section in the data supplement). These SNP sets and their relations with one another characterize the genetic architecture of schizophrenia-associated SNPs in all subjects, including cases and controls (Figure 1A).

Second, we examined the risk of schizophrenia for each SNP set and identified those with high risk. The statistical significance of the association of SNP sets with schizophrenia was calculated using the SNP-Set Kernel Association Test (SKAT) program, which properly accounts for multiple comparisons (15–17).

Third, we checked for significant overlap among SNP sets in terms of subjects and/or SNPs using hypergeometric statistics (24, 25, 28) (see Figures S1 and S2 in the online data supplement). This allowed us to characterize the relations among SNP sets and to identify SNP sets that were connected to each other by having certain SNPs or subjects in common, thereby composing genotypic networks. Disjoint networks shared neither SNPs nor subjects, as expected if schizophrenia is a heterogeneous group of diseases.

Fourth, we identified sets of distinct clinical features that cluster in particular cases with schizophrenia (i.e., phenotypic sets or clinical syndromes) without regard for their genetic background (29), again using non-negative matrix factorization (17). Ninety-three clinical features of schizophrenia from interviews based on the Diagnostic Interview for Genetic Studies (30), as well as the Best Estimate Diagnosis Code Sheet submitted by GAIN/non-GAIN to dbGaP, were initially considered with the MGS sample (see references 31, 32; see also Appendix I in the online data supplement). The Diagnostic Interview for Genetic Studies was utilized for the Portuguese Island samples. Corresponding features were extracted in CATIE from the Positive and Negative Syndrome Scale, the Quality of Life Questionnaire, and the Structured Clinical Interview for DSM-IV (23). These phenotypic sets and their relations with one another characterize the phenotypic architecture of schizophrenia (Figure 1B).

Fifth, we tested whether SNP sets were associated with distinct phenotypic sets in the MGS sample, and we tested the replicability of these relations in the two other independent studies. Replication was evaluated in terms of replication of the SNP sets and their corresponding risk, as well as the relationships between SNP sets and phenotypic sets. In the samples that used the Diagnostic Interview for Genetic Studies (the MGS and Portuguese Island samples), the specific phenotypic features can be compared. Since the CATIE study did not use the Diagnostic Interview for Genetic Studies, we estimated the corresponding symptoms from available phenotypic data (based on the Positive and Negative Syndrome Scale, the Quality of Life Questionnaire, and the Structured Clinical Interview for DSM-IV). Genotypic and phenotypic data were available for 738 cases in CATIE and 346 cases in the Portuguese Island study (see the online data supplement). The significance of cohesive relations among SNP sets and clinical syndromes was tested using hypergeometric statistics (17, 24, 25, 28). The relations between the genotypic and phenotypic clusters characterize the genotypic-phenotypic architecture (Figure 1C).

Methodological details and references for the “divide and conquer” algorithm that we developed and used are available in the online data supplement (24–27). Our web server application PGMRA (17), for identifying genotype-phenotype relations in GWAS, is online at http://phop.ugr.es/fenogeno. Statistical analysis was performed by SKAT (15, 16), also accessible through PGMRA.

Results

Identifying Many SNP Sets as Candidates for Schizophrenia Risk

We first investigated the genotypic architecture of schizophrenia in the MGS study to identify SNP sets without knowledge of the subject’s clinical status (i.e., case or control) (9). Our exhaustive search uncovered 723 nonidentical and possibly overlapping SNP sets in the MGS samples. The SNP sets varied in terms of numbers of both subjects and SNPs. For example, one group contains 70 subjects and 24 SNPs, as expected because few subjects can share a large number of SNPs. Conversely, another group contains 258 subjects and three SNPs, as expected because a large number of subjects are likely to share only a few SNPs. Initially, we retained a large number of SNP sets merely to identify the genotypic clusters in all subjects whether they had schizophrenia or not.

SNP Sets Vary Greatly in Risk for Schizophrenia

Second, we computed the risk for schizophrenia in carriers of each SNP set (33) (Figure 2A–F; see also Figure S3 in the online data supplement). The risk of schizophrenia was normally distributed, as expected when capturing the full range of variability. Ninety-eight of the 723 SNP sets had a risk of schizophrenia greater than 66% and accounted for 90% of schizophrenia cases in the MGS study. Forty-two SNP sets had a risk of schizophrenia ≥70% (Table 1; see also Figure S4 in the data supplement). For example, SNP set 19_2 had a risk of 100%, meaning that all carriers were schizophrenia cases. The ability of SNP sets to predict schizophrenia risk is illustrated in Figure 2G. SKAT showed that the association of schizophrenia with particular SNP sets was stronger than with the average effects of their constituent SNPs (Table 1). For example, the SNP set 81_13 has a p value of 1.46×10⁻¹⁰, whereas the best and average SNPs within this set have p values of 2.15×10⁻⁷ and 5.44×10⁻³, respectively. SKAT and PLINK (34) methods estimated similar p values for the individual SNPs (R²=0.99; p values for F statistics, <3.8×10⁻⁴⁶), showing that SKAT does not inflate results.

The global variance in liability to schizophrenia explained by the average effects of all SNPs simultaneously (8, 35) in our sample was 24%. While individual SNPs were mostly low penetrant, many high-risk SNP sets were highly penetrant (e.g., 100% to 70%; see Table 1) and much more informative in predicting schizophrenia risk.

Relations Among SNP Sets to One Another and to Gene Products

We hypothesized that schizophrenia may be an etiologically heterogeneous group of illnesses in which some genotypic networks are disjoint, that is, share neither SNPs nor subjects. To test this, we first checked for overlap in constituent SNPs and/or subjects among all the SNP sets at high risk for schizophrenia (see Figure S5 in the online data supplement). We found that 17 genotypic networks were disjoint, sharing neither SNPs nor subjects (Figure 3A), suggesting that these are distinct antecedents of schizophrenia. These networks vary in size and complexity: one highly connected network associates 11 SNP sets, whereas eight networks are composed of only a single isolated SNP set.

We also determined that some SNP sets share SNPs but not subjects (e.g., 59_48 and 87_76; Figure 3A), as expected because they involve the same SNPs but with different allele values (both alleles of a SNP can act as risk alleles in different genetic contexts). In contrast, we found that the 58_29 and 41_12 SNP sets do not share SNPs, but independently specify almost the same individuals (Figure 3A), as expected when, for example, distinct subsets of genotypic features influence a common developmental pathway. Finally, some SNP sets overlap in both SNPs and subjects, suggesting that one is a subset within the other (e.g., 88_64 and 81_13; see Figure S3A,C in the online data supplement). Therefore, the genotypic networks display distinct topologies differing in the way constituent SNPs and subjects are related.

When evaluating whether different genotypic networks operate through distinct mechanisms, we found that high-risk SNP sets mapped to various classes of genes (e.g., protein coding, ncRNA genes, and pseudogenes) related to known functions and causing different effects on their products (Figure 3A; see also Tables S1–S3 and Figure S6 in the online data supplement). We identified distinct pathways as exemplified in Table 2. Notably, all of these pathways are interconnected by the overlapping gene products that include genes previously associated with schizophrenia by GWAS, as well as genes known to be abnormally expressed in the brains of schizophrenia patients (see Table S4, Figure S7, and the Pathways section in the data supplement). The emerging picture is suggestive of a possible pathophysiology in which abnormal brain development interacts with environmental events triggering abnormal or exaggerated immune and oxidative processes that increase risk of schizophrenia.

Complex Genotypic-Phenotypic Relationships in Schizophrenia

Next we examined whether the complex genetic architecture of schizophrenia leads to phenotypic heterogeneity. Using data from the Diagnostic Interview for Genetic Studies (30), as well as from the Best Estimate Diagnosis Code Sheet submitted by GAIN/non-GAIN to dbGaP (see Appendix I, Figures S1 and S2, and the Method section in the online data supplement), we originally identified 342 nonidentical and possibly overlapping phenotypic sets of distinct clinical features that cluster in particular cases with schizophrenia (i.e., phenotypic sets or clinical syndromes) without regard for their genetic background. Different SNP sets were significantly associated with particular clinical syndromes (hypergeometric statistics, p values from 2×10⁻¹³ to 1×10⁻³). However, the genotypic-phenotypic relations were complex (i.e., many-to-many [29]): the same genotypic network could be associated with multiple clinical outcomes (i.e., multifinality or pleiotropy) and different genotypic networks could lead to the same clinical outcome (i.e., equifinality or heterogeneity; Table 3; see also Table S5 in the data supplement). The genotypic-phenotypic relations were highly significant by a permutation test (empirical p value <4.7×10⁻³; Table 3; see also Table S5).

Specifically, we identified a phenotypic set indicating a general process of severe deterioration (i.e., continuous positive symptoms with marked and progressive impairment) that was associated with many SNP sets (e.g., SNP sets 75_67 and 56_30, with p values <2.3×10⁻¹³ and 2.55×10⁻⁵, respectively; Table 3, Figure 3A). Other SNP sets were associated with a general process of moderate deterioration (moderate or fluctuating impairment despite a continuous mixture of symptoms), as in SNP sets 14_6, and 42_37 (p values <5×10⁻⁴; Table 3, Figure 3A).

We identified specific clinical syndromes that were unambiguously associated with particular genotypic networks. For example, specific phenotypic sets differentiate among SNP sets even within the same network, which illustrate similar but not identical forms of multifinality in schizophrenia (e.g., 76_74 and 58_29; Table 3, Figure 3A, blue lines). Particular phenotype sets can also distinguish SNP sets connected only by shared subjects (Figure 3A, red lines). For example, SNP set 76_74 shares subjects with 56_30 and with 81_13; however, the latter SNP sets are associated with a specific phenotypic set not present in 76_74 (Table 3).

Positive and Negative Symptoms Differentiate Classes of Schizophrenia

Genotypic and phenotypic relationships could be grouped into eight classes of schizophrenia, as shown in Figure 3B and Table 3 (31, 32, 36). First, we identified SNP sets involving subjects with predominantly positive symptoms (e.g., 41_12 and 88_64) and few residual symptoms. Second, we identified SNP sets represented by predominantly negative and disorganized symptoms (e.g., 10_4 and 61_39), decreased psychosocial function, and continuous residual symptoms. As discussed in the online data supplement (see the Replicability of the Phenotypic Features section), bizarre delusions and symptoms of cognitive and behavioral disorganization, such as thought insertion and disorganized speech among others, were accepted as fuzzy indicators of either positive or negative classes of schizophrenia but were considered to be more common in negative and disorganized classes (e.g., in Table 3, thought echo and commenting hallucinations in “negative schizophrenia” with phenotypic set 46_29 associated with SNP set 14_6).

Third, several SNP sets harbor mixed positive and negative symptoms (e.g., 59_48 and 54_51). These three classes were enriched by considering the general severe and moderate patterns, which were frequent in several networks (Figure 3B), as described above. Because the latter patterns appear in combination with a set of only positive symptoms (e.g., 81_13), both positive and negative symptoms (e.g., 75_67), and only negative symptoms (e.g., 19_2), we were able to classify schizophrenia into eight classes (Figure 3B). A principal-components analysis of the phenotypic features in the Diagnostic Interview for Genetic Studies confirmed this classification (see Table S6 and the Method section in the online data supplement).

Replication of Results in Two Independent Samples

We tested the replicability of our findings in the MGS study by carrying out the same analyses of the genotypic and phenotypic architecture of schizophrenia in the CATIE (19, 22, 23) and Portuguese Island (19, 21) samples. A total of 1,303 SNPs were shared between the selected SNPs in the MGS (see the Data Cleaning section in the online data supplement) and CATIE samples, and 1,234 SNPs between the MGS and Portuguese Island samples. Imputed variants were not considered, to avoid possible biases.

We found that 31 and 30 of the 42 SNP sets selected in the MGS sample were also identified in the CATIE and Portuguese Island samples, respectively (see Tables S7 and S8 in the online data supplement). Together, both samples reproduced at least 81% of the SNP sets at risk (see Table S9 in the data supplement). In addition, most of the SNP sets replicated in the two PGC samples achieved risk values as high as those of the MGS sample (>70%) (see Table S8): 70% of those identified exhibit >70% risk, and 90% show >60% risk. Some SNP sets exhibited slightly higher risk values than those in the MGS sample.

The genotypic-phenotypic relations in CATIE and the Portuguese Island studies closely matched those observed in the MGS study (hypergeometric statistics, p values 1×10⁻⁷ to 1×10⁻²; see Tables S7 and S8 and the Replicability section in the data supplement). The eight schizophrenia classes exhibited high reproducibility. For example, except for one relation (“–” in the MGS study and “+ and –” in CATIE; see Table S9 in the data supplement), all relations exhibited similar positive and negative symptoms in the MGS study and CATIE. Three relations showed less specific symptoms in CATIE than in the MGS study, as expected because CATIE did not use the Diagnostic Interview for Genetic Studies (see Table S10 and the Replicability section in the data supplement).

We found few differences when comparing the MGS and Portuguese Island studies (see Table S9 in the data supplement), except differences in severity that preserved the sign of the symptoms. Three relations with negative symptoms in the MGS study exhibited negative and positive symptoms in the Portuguese Island sample (see Table S9). Only two SNP sets in the Portuguese Island sample had no significant cross-match with the phenotypic features expected from the MGS study.

Discussion

Our findings indicate that schizophrenia comprises several distinct clinical syndromes associated with many disjoint genotypic networks. Consequently, much of the heritability of schizophrenia has not been detected by approaches that classify subjects only according to whether or not they have schizophrenia. Our purely data-driven analysis shows that the elusive heritability of schizophrenia is not missing, but is encoded in a complex distribution of genotypic-phenotypic relationships.

We found that 42 interactive SNP sets had greater than 70% risk of schizophrenia. The interactive SNP sets explained the risk more fully than the average effects of all SNPs simultaneously and were more strongly related to their particular syndromes of schizophrenia than are their individual SNPs (Table 1). Consequently, identifying the organization of SNPs into interactive SNP sets enabled us to increase the power to detect associations: 98 SNP sets with greater than 66% risk accounted for 90% of cases. The constituent genes in these networks belong to signaling pathways highly associated with schizophrenia (see Figure S7 in the online data supplement). Our findings have broad implications, so we will consider their strengths and limitations carefully.

Strengths and Limitations

Two particular features of our methods merit consideration in terms of their strengths and limitations. First, we concurrently used detailed assessments of both the genotype and the phenotype to identify their associations, thereby combining genomic and phenomic information (29). Other approaches decrease the number of variables before analysis (“data reduction”), even if the biological importance of these variables is not known a priori. The evidence we have that schizophrenia is a heterogeneous group of disorders suggests that reducing clinical information about schizophrenia to a single categorical diagnosis is inadequate.

Despite the detailed phenotypic information we had available about subjects, there are still limitations to data obtained even from reliable structured interviews like the Diagnostic Interview for Genetic Studies. Interview data are based on self-reports that are interpreted and coded by interviewers. Subjects may not be willing or able to report their symptoms accurately. We had obtained information from treatment records and family history reports, but we chose not to use such additional information, except for the resulting best-estimate final DSM ratings of diagnosis, because its extent and quality varied in unmeasured ways between cases. The greatest limitation in the phenotypic assessments in available GWAS databases has been the overreliance on subjective symptoms with an absence of objective measurements, such as cognitive tests, brain electrophysiology, and neuroimaging (37). Subjective symptoms are fuzzy indicators of the underlying pathophysiology. Objective measures could complement the assessment of symptoms and could be applied to both cases and controls, thereby providing a more comprehensive and valid characterization of the phenotype of all subjects. The biggest challenge in GWAS is access to studies with rich phenotypic data about both subjective and objective measures obtained systematically from all subjects.

Our finding of robust replicability based on detailed symptom profiles from interview data alone do have important implications for the size of samples and the scope of phenotypic assessments in genomic studies of complex disorders. We obtained robust replication of results in moderate size samples, such as 738 cases in CATIE and 346 cases in the Portuguese Island study, which shows that it is incorrect to assume that extremely large samples are needed to obtain robust and replicable findings. Difficulty in replication in previous work can be better explained by the neglect of the complexity of genetic and phenotypic architecture rather than by moderate sample size. We identified more information by combining rich data about the complex architecture of genotypes and interview-based symptoms in such moderate size samples than has been obtained in analysis of much larger compilations of multiple samples that relied on additive gene effects on categorical diagnosis (8).

Nearly all previous genetic association studies have relied on patient interviews for clinical description, as detailed objective testing has been impractical in large samples because of cost and the difficult logistics of securing cooperation with time-consuming test batteries. Now that we have shown that replicable results can be obtained in moderate size samples, it is feasible to complement interview data with more objective and thorough assessments. Fundamental research into the causes and characteristics of the schizophrenias is likely to require phenotypic assessment beyond the clinical features needed for clinical diagnosis and treatment, which has given primacy to signs and symptoms assessable by interview alone (37).

Second, we have strived to extract the maximum information available in a single GWAS without making restrictive a priori assumptions. In other words, we proceeded in a data-driven, model-free manner. As a consequence, whatever information emerged from the data mining process (such as the different classes of schizophrenia) is inherent to the data and was not artificially imposed by either an a priori model or previous knowledge of the data (such as the “case or control” status of the subjects). Nevertheless, our initial pool of 2,891 SNPs, preselected for at least loose association with schizophrenia in the MGS study (9), might be missing additional risk SNPs that would eventually show up in an even more exhaustive genomic analysis.

Our findings about the heterogeneity and complexity of schizophrenia (31, 32, 36) require a careful reconsideration of the concept of “replicability.” In order to be meaningful in complex disorders like schizophrenia, efforts to replicate findings must take into account the distributed heritability and developmental complexity of the disease.

Replication: A Lock and Key Combination of Genomics and Phenomics

Replication is always critical, but it is not usually sought within a single large study. Here, internal replicability was addressed by resampling techniques (94% support; see the online data supplement), where the same SNP sets are systematically identified despite the random alteration of the parameters of the method (17) and/or the sample (38). In addition, our biggest challenge was to identify studies with rich phenotypic data for independent external replication. In most GWAS, phenotypic data have been of “secondary” interest, using a variety of structured or even unstructured interviews (14, 18, 19) (see the Replicability section in the data supplement). So why not attempt to replicate the genotypic architecture alone? The same answer applies for any method for validation of associations: genetic variants associated with individuals may be, and in all likelihood often are, completely unrelated to the disease. The only way to make sense of these associations is to cross-match genomics with high-resolution phenomics (29). One can think of it as a “lock and key” combination (or, more precisely, many such combinations), where both pieces of information are needed to be able to interpret the results with confidence. Note that our approach complements meta-analysis (39) and/or pathway analyses (40), focusing the search on the combined genotypic-phenotypic architecture.

Despite the described constraints, we successfully identified more than 81% of the genotypic-phenotypic relationships previously found in the MGS data set in two independent samples. These samples were the only ones where both genotypic and detailed phenotypic features were available and provided by the researchers. Remarkably, the identification was performed with half of the SNPs used in the MGS study, because of the different platforms and our conservative preference to avoid external imputations. The success of our replication efforts strongly supports the validity and power resulting from combining genomic and phenomic information in association studies.

Overall, we believe our approach is a pioneering effort to specify complex but manageable patterns of gene-gene interaction underlying the polygenic risk of schizophrenia. In addition, our results hold promise for the emergence of a new era in clinical psychiatry in which person-centered treatment of complex disorders can be guided by reliable assessments of well-validated clinical syndromes and their specific causes.