Introduction
A major goal in mental health research is to identify endophenotypes that are neurobehavioral measures of psychiatric disorder liability. While complex in their application, endophenotypes are conceptually simple: in the causal chain between genes and clinical outcome, endophenotypes are a laboratory-based measure of a component of that chain that are causally closer to gene action than the eventual clinical outcome. Endophenotypes are usually defined as laboratory-based measures that are associated with a clinical outcome, heritable, manifest in an affected individual even when the disorder is not active, and co-segregate with the associated clinical outcome in families (Gottesman & Gould, Reference Gottesman and Gould2003). We have also argued previously that to qualify as an endophenotype, versus a putative endophenotype, the laboratory-based measure should also show robust and reliable associations with genetic variants through genetic association studies (Iacono et al. Reference Iacono, Malone and Vrieze2016). Endophenotypes are also thought to be more genetically homogeneous than their associated clinical outcomes (Cannon & Keller, Reference Cannon and Keller2006), and so effects of individual genetic variants on endophenotypes are suspected to be larger and possibly easier to detect than the effect of those variants on the clinical outcome. Recent work suggests, however, that the effects of individual genetic variants on endophenotypes may not be appreciably larger than for complex and distal psychiatric or medical outcomes (Iacono et al. Reference Iacono, Malone and Vrieze2016). Using a community sample of over 4900 individuals from the Minnesota Twin Family Study (MTFS) (Iacono et al. Reference Iacono, Malone, Vaidyanathan and Vrieze2014a ), we have previously evaluated the genetic basis of 17 putative psychophysiological endophenotypes, including antisaccade eye movements, P300, electrodermal activity (EDA), resting electroencephalogram (EEG), and acoustic eyeblink startle, all of which were found to have moderate to high heritability and which, in past research, have demonstrated associations with various clinical outcomes, including schizophrenia. We conducted genome-wide association analyses (Malone et al. Reference Malone, Burwell, Vaidyanathan, Miller, McGue and Iacono2014a , Reference Malone, Vaidyanathan, Basu, Miller, McGue and Iacono b ; Vaidyanathan et al. Reference Vaidyanathan, Isen, Malone, Miller, McGue and Iacono2014a , Reference Vaidyanathan, Malone, Donnelly, Hammer, Miller, McGue and Iacono b , Reference Vaidyanathan, Malone, Miller, McGue and Iacono c ), analyses of array-based non-synonymous exonic variants (Vrieze et al. Reference Vrieze, Malone, Pankratz, Vaidyanathan, Miller, Kang, McGue, Abecasis and Iacono2014b ) and 11× whole genome sequencing (Vrieze et al. Reference Vrieze, Malone, Vaidyanathan, Kwong, Kang, Zhan, Flickinger, Irons, Jun, Locke, Pistis, Porcu, Levy, Myers, Oetting, McGue, Abecasis and Iacono2014c ). Ultimately, we were neither able to convincingly corroborate any prior finding for any gene or variant previously associated with any of the 17 psychophysiological endophenotypes nor were we able to conclusively identify novel genetic associations. Our findings suggested that endophenotypes will not easily deliver novel genetic discoveries of relevance to clinical psychopathological outcomes.
The present study is a crucial direct extension of our prior work with these 17 endophenotypes. Here, the question is not whether psychophysiological endophenotypes can aid in initial discovery of disorder-relevant genes, but whether such endophenotypes can inform us about the role and mechanisms of action of genes already known to be associated with a major neurodevelopmental mental disorder. To answer this question, we took the results from a recent schizophrenia genetic association meta-analysis by the Psychiatric Genomics Consortium (PGC), which discovered 128 independent genetic variants in 108 genomic loci associated with schizophrenia (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014), and tested them for association with the same 17 endophenotypes described above. Many of these endophenotypes have demonstrated substantial construct validity as endophenotypes for schizophrenia (Iacono et al. Reference Iacono, Malone and Vrieze2016); others less so, but nevertheless have been associated with schizophrenia, found to be heritable in twin studies, and show deviation in schizophrenia patient first-degree relatives (Fukushima et al. Reference Fukushima, Fukushima, Chiba, Tanaka, Yamashita and Kato1988; Braff et al. Reference Braff, Grillon and Geyer1992; Iacono et al. Reference Iacono, Ficken and Beiser1999; Hong et al. Reference Hong, Summerfelt, Mitchell, O'Donnell and Thaker2012; Kam et al. Reference Kam, Bolbecker, O'Donnell, Hetrick and Brenner2013; Chen et al. Reference Chen, Lee, Yang, Landau, Chang, Chen, Lu, David and Bramon2014; Narayanan et al. Reference Narayanan, O'Neil, Berwise, Stevens, Calhoun, Clementz, Tamminga, Sweeney, Keshavan and Pearlson2014; Qiu et al. Reference Qiu, Tang, Chan, Sun and He2014). The only two endophenotypes that have not been linked to schizophrenia are our measures of affectively modulated startle (Iacono et al. Reference Iacono, Malone and Vrieze2016). All 17 of our endophenotypes are measures of basic neurological and physiological function (e.g. resting EEG, EDA, startle) that, broadly speaking, may be expected to be disrupted in the presence of high genetic liability for schizophrenia.
The MTFS is a powerful sample in which to test endophenotype utility, especially for schizophrenia and other relatively rare disorders. Endophenotypes are hypothesized to be present in unaffected family members of probands, and before the advent of cheap genotyping and known risk loci, familial relatedness was an excellent proxy for genetic risk for disorder. With the help of the PGC results we can rank-order individuals in a community sample by their predicted genetic risk according to measured genetic variants. Thus, even though schizophrenia rates in the present sample are nearly zero, we still have an estimate of each individual's genetic liability to schizophrenia and, therefore, to endophenotypes truly related to schizophrenia. In the present article we use this sample combined with array-based genotypes, whole genome sequencing, and a variety of statistical approaches to test for relationships between the 108 PGC schizophrenia-associated loci and 17 putative endophenotypes to address these questions:
-
1. Are any of the 128 PGC schizophrenia-associated single nucleotide polymorphisms (SNPs) also associated with any of our 17 endophenotypes, especially those endophenotypes that show robust phenotypic associations with schizophrenia?
-
2. Is an aggregate index of schizophrenia genetic liability (obtained by summing risk alleles, weighted by strength of association from the PGC findings) associated with any of these endophenotypes?
-
3. Are there rare variants within the 108 PGC-identified loci that have reliable associations with any of these endophenotypes?
Method
Participants
All participants were taken from the Minnesota Twin Family Study (MTFS), a community-representative study of twins and their parents. Participants completed psychiatric assessments, self-report questionnaires and laboratory-based measures (Braff et al. Reference Braff, Grillon and Geyer1992; Iacono et al. Reference Iacono, Ficken and Beiser1999, Reference Iacono, Malone, Vaidyanathan and Vrieze2014a , Reference Iacono, Malone and Vrieze2016). The total number of individuals available with genotypes and at least one psychophysiological measure was 4905. All individuals were genotyped on the Illumina 660W-Quad, with genotyping protocols and quality control described extensively previously (Miller et al. Reference Miller, Basu, Cunningham, Eskin, Malone, Oetting, Schork, Sul, Iacono and McGue2012; Vrieze et al. Reference Vrieze, Feng, Miller, Hicks, Pankratz, Abecasis, Iacono and McGue2014a ). After quality control, the Illumina 660W-Quad contained a total of 527 829 variants.
Of those individuals with array genotypes, 1431 were selected for 11× whole genome sequencing. Among the 1431 sequenced participants, 406 were members of a monozygotic twin pair and, for each monozygotic twin pair, we sequenced only one member of the pair and copied the sequence to the other twin resulting in a total of 1837 individuals for whom whole-genome sequence variant calls were available. Monozygotic twin genomes were copied under the assumption that their genomes are identical or, at the very least, so similar that any genetic differences between the two twins would be indistinguishable from error at 11× depth sequencing. Table 1 contains descriptive information on the sample. Online Supplementary Table S1 shows the 17 × 17 endophenotype correlation matrix. All of the endophenotypes shows moderate to strong twin-based heritability except the two startle difference scores.
Descriptive statistics for the 17 endophenotypes a
CI, Confidence interval; GCTA, genetic complex trait analysis; SNP, single nucleotide polymorphism; s.e., standard error.
a The heritability estimates were from a series of previous papers (Malone et al. Reference Malone, Burwell, Vaidyanathan, Miller, McGue and Iacono2014a , Reference Malone, Vaidyanathan, Basu, Miller, McGue and Iacono b ; Vaidyanathan et al. Reference Vaidyanathan, Isen, Malone, Miller, McGue and Iacono2014a , Reference Vaidyanathan, Malone, Donnelly, Hammer, Miller, McGue and Iacono b , Reference Vaidyanathan, Malone, Miller, McGue and Iacono c ). The relatedness threshold used in that prior work to estimate the SNP heritability was 0.025.
Measures
We evaluated 17 psychophysiological measures, all of which, except affectively modulated startle, have received support as endophenotypes for schizophrenia. All measures have been described in detail previously (Iacono et al. Reference Iacono, Malone, Vaidyanathan and Vrieze2014a ; Vrieze et al. Reference Vrieze, Malone, Pankratz, Vaidyanathan, Miller, Kang, McGue, Abecasis and Iacono2014b ). All measures were corrected for sex, age, parent/child and technical artifacts. Measures were inverse-normalized prior to genetic association to mitigate outlier effects.
Antisaccade
Participants viewed a spot of light in the center of a screen. A light then appeared in the left or right visual field and participants were instructed to look away from the stimulus. Antisaccade here is measured as the proportion of trials in which the participant initially looked toward the spot of light, rather than away from it. Antisaccade eye movement has long been a strong candidate endophenotype for schizophrenia, replicated in dozens of studies over recent decades (Radant et al. Reference Radant, Dobie, Calkins, Olincy, Braff, Cadenhead, Freedman, Green, Greenwood, Gur, Gur, Light, Meichle, Millard, Mintz, Nuechterlein, Schork, Seidman, Siever, Silverman, Stone, Swerdlow, Tsuang, Turetsky and Tsuang2010; Leonard et al. Reference Leonard, Robinson, Kaiser, Hahn, McClenon, Harvey, Luck and Gold2013). Antisaccade eye tracking error is highly heritable (Vaidyanathan et al. Reference Vaidyanathan, Malone, Donnelly, Hammer, Miller, McGue and Iacono2014b ), present in first-degree relatives of individuals with schizophrenia, and co-segregates within families (Calkins et al. Reference Calkins, Curtis, Iacono and Grove2004; Radant et al. Reference Radant, Dobie, Calkins, Olincy, Braff, Cadenhead, Freedman, Green, Greenwood, Gur, Gur, Light, Meichle, Millard, Mintz, Nuechterlein, Schork, Seidman, Siever, Silverman, Stone, Swerdlow, Tsuang, Turetsky and Tsuang2010).
P300 event-related potential (ERP)
Participants completed a rotated heads visual oddball task (Begleiter et al. Reference Begleiter, Porjesz, Bihari and Kissin1984). Amplitude of the P300 wave was obtained at a midline parietal electrode (Pz) from the average ERP waveform derived from all target trials. P3 genetic factor score was generated from a factor analysis of P300 amplitude from multiple scalp electrode sites in four-member families (Malone et al. Reference Malone, Vaidyanathan, Basu, Miller, McGue and Iacono2014b ). P300 has been studied extensively in schizophrenia with visual and auditory tasks in case–control studies producing over 100 journal articles in the last 10 years (Bramon et al. Reference Bramon, Rabe-Hesketh, Sham, Murray and Frangou2004). Individuals with schizophrenia tend to have reduced P300 amplitude during the decision-making process as compared with controls (Jeon & Polich, Reference Jeon and Polich2003).
Resting EEG
Participants were asked to relax with eyes closed for 5 min while listening to soft white noise. EEG was recorded from the vertex of the head (Cz), referenced to the two ears, and from two bipolar derivations: O1-P7 and O2-P8. Power in the alpha, beta, theta and delta frequency bands were obtained from a fast Fourier transform of EEG recorded at the Cz electrode. Alpha power and peak frequency were also obtained for the two bipolar leads and averaged together. Schizophrenia is robustly associated with increased low-frequency power in particular (Kam et al. Reference Kam, Bolbecker, O'Donnell, Hetrick and Brenner2013; Narayanan et al. Reference Narayanan, O'Neil, Berwise, Stevens, Calhoun, Clementz, Tamminga, Sweeney, Keshavan and Pearlson2014), a finding that is also observed in their relatives (Narayanan et al. Reference Narayanan, O'Neil, Berwise, Stevens, Calhoun, Clementz, Tamminga, Sweeney, Keshavan and Pearlson2014).
EDA
EDA was measured from the fingertips as participants viewed scenes from a closed-captioned movie while loud tones were intermittently delivered. Four measures of EDA were derived. Skin conductance level was a measure of the participants’ tonic resting level measured at the end of the session. Skin conductance response frequency is a count of the number of tones to which the participant responded. Skin conductance response amplitude is the mean magnitude for trials that participants had an observable response. EDA factor score is a global measure using a factor score derived from a model fit to the three skin conductance measures. Electrodermal hyporesponsiveness is characteristic of schizophrenia (Bernstein et al. Reference Bernstein, Frith, Gruzelier, Patterson, Straube, Venables and Zahn1982), and deviations in skin conductance activity are present in the first-degree relatives of schizophrenia patients (Iacono et al. Reference Iacono, Ficken and Beiser1999).
Acoustic startle
Participants viewed a series of standard images from the International Affective Picture System (Lang et al. Reference Lang, Bradley and Cuthbert2008). Images are categorized as affectively pleasant or aversive or as affectively neutral. Eye blink reactions to loud noise probes during viewing of the images were recorded from the orbicularis oculi muscle, and three measures of reactivity derived. Overall startle is the overall magnitude of the integrated electromyographic response averaged over all the trials; overall startle is moderately heritable in this MTFS sample (Vaidyanathan et al. Reference Vaidyanathan, Isen, Malone, Miller, McGue and Iacono2014a ) and the acoustic startle response has been reported as deviant in schizophrenia patients (Braff et al. Reference Braff, Grillon and Geyer1992). Aversive difference startle is the z-score difference in mean startle magnitude between aversive and neutral images whereas pleasant difference startle is the difference between the pleasant and neutral images.
Selected PGC variants
All variants were selected based on their inclusion in the PGC's 2014 meta-analysis of schizophrenia. The variant list was downloaded in July 2014 (https://www.med.unc.edu/pgc/downloads).
In their meta-analysis, the PGC conducted imputation using 1000 Genomes phase 1, the standard reference panel available at that time. In the present study, we used 1000 Genomes phase 3 as our imputation reference panel, as phase 3 is an improvement over phase 1. Imputation was conducted using Minimac3 (Das et al. Reference Das, Forer, Schönherr, Sidore, Locke, Kwong, Vrieze, Chew, Levy, McGue, Schlessinger, Stambolian, Loh, Iacono, Swaroop, Scott, Cucca, Kronenberg, Boehnke, Abecasis and Fuchsberger2016) on an imputation webserver (https://imputationserver.sph.umich.edu). In total, 45 994 018 variants were imputed, including 124 of the 128 PGC genome-wide significant variants. Imputation quality is shown in online Supplementary Table S2. Of these 124 variants, 87% had a Minimac R 2 (RSQ) higher than 0.9. The four remaining variants failed because they did not exist in 1000 Genomes phase 3. According to the 1000 Genomes Project Consortium, a total of 2.3 million variants were present in phase 1 but absent from phase 3 (http://www.1000genomes.org/faq/why-are-phase-1-sites-missing-phase-3-dataset/). Differences between phase 1 and 3 are due to platform differences, a few samples from phase 1 not included in phase 3, and improved variant calling pipelines.
Of the 128 variants, 124 were therefore successfully imputed. However, one of the 124 variants (an indel on chromosome 6) was monomorphic among those individuals in the MTFS who participated in the psychophysiology laboratory. To replace the monomorphic and remaining four variants, we used surrogate variants from 1000 Genomes phase 3 that were in high linkage disequilibrium with the missing variants but were still significantly associated with schizophrenia in the PGC results. We used rs281768 (association with schizophrenia in PGC: p = 3.47 × 10−14) as a surrogate for an indel on chromosome 2, position 200 825 237 (p = 1.78 × 10−14); rs3798869 (p = 1.15 × 10−9) for an indel on chromosome 6, position 84 280 274 (p = 8.57 × −10); and an indel on chromosome 22, position 39 987 024 (p = 8.91 × 10−11) for an indel on chromosome 22, position 39 987 017 (p = 2.20 × 10−11). Since two of the unimputed variants (indel at chromosome 10, position 104 957 618 and rs7907645; p = 1.04 × 10−13 and p = 5.82 × 10−11, respectively) were within the same locus and we do not have access to conditional genetic association results from the PGC, we replaced those two variants with a single variant (rs12414028; PGC p = 1.96 × 10−12). In sum, we analysed 123 variants directly and four surrogate variants out of the 128 total target PGC schizophrenia-associated variants.
Whole genome sequencing
Sequencing details are included in the online Supplementary material. We used the integrated GotCloud program for alignment and variant calling (Jun et al. Reference Jun, Wing, Abecasis and Kang2015). After alignment, read clipping and duplicate removal, average sequencing depth genome-wide was 11×. Quality-control information is shown in online Supplementary Table S3. Variant calling was done with GotCloud defaults, using Beagle for linkage disequilibrium refinement of called genotypes (online Supplementary Table S4). To estimate power of 11× sequencing for rare variant discovery we compared the sequence variant calls to Illumina HumanExome array variant calls. The sequenced calls identified 11 484 of the 15 791 singletons on an Illumina HumanExome array genotyped on the same sample (73%), 10 308 of the 10 990 doubletons (94%) and 6488 of 6573 tripletons (99%). Of 151 860 monomorphic sites on the exome array, sequencing called 2951 of these sites polymorphic (about 2%). We next estimated genotype concordance between the sequence variant calls and the Illumina HumanExome and 660W-Quad genotypes. Autosomal genotypes showed 99.91% concordance (online Supplementary Fig. S1). In sum, the sequence data provided accurate genotype calls, even for rare variants including singletons and doubletons.
Single variant association tests
For each of the 127 successfully imputed schizophrenia-associated variants we tested association with each of the 17 endophenotypes. We also conducted single variant tests of all whole-genome-sequenced variants within the 108 PGC-identified loci. In both cases, we used a linear mixed model, EMMAX, that accounts for both family structure and population stratification using an empirical kinship matrix (Kang et al. Reference Kang, Sul, Service, Zaitlen, Kong, Freimer, Sabatti and Eskin2010).
Gene-based aggregate tests
We conducted gene-based burden tests of rare non-synonymous variants identified in the whole genome sequences. Non-synonymous variants were annotated using EPACTS (Kang et al. Reference Kang, Sul, Service, Zaitlen, Kong, Freimer, Sabatti and Eskin2010) against GENCODE v11. Two types of gene-based tests were conducted, a variable threshold burden test (VT), and the sequence kernel association test (SKAT; Price et al. Reference Price, Kryukov, de Bakker, Purcell, Staples, Wei and Sunyaev2010; Wu et al. Reference Wu, Lee, Cai, Li, Boehnke and Lin2011). In the VT approach, each participant is assigned a burden score for each gene. The burden score is simply the sum of the number of non-synonymous minor alleles under a certain minor allele frequency (MAF). The MAF is allowed to vary until the sum score shows the strongest association with the trait, and the resulting p value is corrected for the number of MAF thresholds considered. The SKAT approach is more complex, but the basic intuition is that it conducts a test of the variance of variant effects within a gene, and whether that variance is larger than expected by chance. That is, do the variant effects differ from zero more than expected by chance? In SKAT we tested both MAF = 0.01 and MAF = 0.05. The purpose of using both VT and SKAT is that they each achieve greater statistical power in complementary scenarios. VT is most powerful when all minor alleles within a gene have the same direction of effect and has low power when minor alleles have opposite directions of effect. SKAT has the same power regardless of whether the minor alleles within a gene have opposite or the same directions of effect.
Polygenic risk scores
Finally, we computed polygenic risk scores and tested for association between the score and each endophenotype. A polygenic risk score is a weighted sum of risk alleles across many variants, with the weight of each variant equal to the effect size for that variant in the PGC schizophrenia results. More formally, for individual i, their polygenic risk score for variants j through n is:
$${\rm PRS}_i = \sum\limits_{\,j = 1}^n \ln (\hbox{PGC OR}_{ij}) \times(\hbox{no. of risk alleles}_{ij} ).$$
The odds ratio (OR) and effect allele coding were taken from the PGC results described earlier. First, we considered a polygenic risk score from the imputed PGC genome-wide-significant variants, corresponding to a p value threshold of 5 × 10−8. At this threshold, we included all the variants (and their surrogates) identified by the PGC. We then computed polygenic risk scores with increasing numbers of variants by relaxing the inclusion p value threshold at each order of magnitude below 5 × 10−8. For these scores, we only considered variants where the reference and alternate alleles matched between our imputed genotypes and the PGC-reported genotypes. Of the variants, 93.12% matched, and the majority of variants that did not match were indels or multi-allelic sites. We chose to exclude indels due to their allelic complexity and the multiple ways that indel alleles have been coded in sequence VCF files (Tan et al. Reference Tan, Abecasis and Kang2015). It would be difficult if not impossible to ensure that we accurately coded effect alleles for indels and multi-allelic sites for all the variants. We also removed low-frequency SNPs (MAF < 10%) and low-quality imputed variants (Minimac RSQ < 0.9). For each p value threshold, we then rank-ordered all candidate SNPs by reported p value for association with schizophrenia in the PGC. We considered the most significant SNP and removed all variants in linkage disequilibrium r 2 > 0.1 within 1 megabase of this variant from our rank-ordered list. We then moved on to the next most significant variant and again removed variants in linkage disequilibrium, repeating this process until no more variants remained in the rank-ordered list. This approach is similar to that used in the PGC's own analysis (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). To derive statistically valid p values for the correlations between the polygenic risk scores and each endophenotype in the context of our twin design, we estimated in OpenMx (Boker et al. Reference Boker, Neale, Maes, Wilde, Spiegel, Brick, Spies, Estabrook, Kenny, Bates, Mehta and Fox2011) correlations and their 95% confidence intervals (Neale & Miller, Reference Neale and Miller1997) using a standard twin family model (Martin & Eaves, Reference Martin and Eaves1977) with a shared environmental component estimated in the twins and no shared environment in the parent generation.
Results
Single variant associations of 127 PGC variants
Accounting for multiple testing only within a single endophenotype (0.05/127 = 3.9 × 10−4), we found two SNPs associated with an endophenotype. SNP rs17194490 (reference = G, alternate = T) was associated both with antisaccade and P300. The alternate allele decreased antisaccade errors (β = −0.12, p = 2.3 × 10−4, alternate allele frequency = 0.163) and increased amplitude on the P300 genetic factor (β = 0.12, p = 2.7 × 10−4). The alternate allele in SNP rs7267348 (reference allele = T, alternate = C) decreased beta EEG power (β = −0.11, p = 3.8 × 10−4, alternate allele frequency = 0.25). Neither variant survived multiple testing for more than a single endophenotype. Summary statistics for all 127 variant associations with all 17 endophenotypes are shown in online Supplementary Table S5.
Whole genome sequence fine-mapping of the 108 PGC loci
In the whole genome sequence variant call set we discovered 203 905 variants within the 108 loci identified by and defined by the PGC (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). We again employed a liberal multiple-testing correction (0.05/203 905 = 2.5 × 10−7) that corrected only for the number of SNPs within these loci and ignored multiple phenotypes. Even at this liberal threshold, we found no significant association between any variant in these loci and any of the 17 endophenotypes.
Gene-based associations
We identified 403 genes that are either within or overlap the 108 loci identified by the PGC. Once more at a liberal threshold of 1.2 × 10−4, correcting only for the number of genes and ignoring multiple endophenotypes, the VT identified no significant association between any gene and any of the endophenotypes. SKAT with a MAF threshold of 0.01 also failed to identify associations, although changing the MAF threshold to 0.05, we report one gene (strawberry notch homolog 1; SBNO1) on chromosome 12 that was significantly associated with pleasant difference startle (p = 4.19 × 10−5), although this did not survive correction for multiple testing of all endophenotypes.
Polygenic risk score
Fig. 1 illustrates results from the polygenic risk score analysis. The figure shows, for each p value threshold the number of variants included in the risk score, the p value threshold, and the correlation estimate between the endophenotype and the polygenic risk score. The only polygenic risk score association approaching significance was that between the score based on the top 127 variants and overall startle (p ~ 0.005), which does not survive a correction for multiple testing of 17 endophenotypes (0.05/17 = 0.003), let alone for multiple testing of 17 endophenotypes and nine p value thresholds.
Correlations between schizophrenia polygenic risk scores (PRS) and each of the endophenotypes. To improve readability, the figure is organized into four panels, where each panel shares the x and y axes. No correlation between polygenic risk and any endophenotype was significant after accounting for multiple testing. EEG, Electroencephalogram.
Discussion
One argument for the utility of endophenotypes has been that they may aid in gene discovery for clinical outcomes. Previous genome-wide association work by our group suggests that this is not true for common psychophysiological endophenotypes in large community samples (Iacono et al. Reference Iacono, Vaidyanathan, Vrieze and Malone2014b ; Vrieze et al. Reference Vrieze, Malone, Vaidyanathan, Kwong, Kang, Zhan, Flickinger, Irons, Jun, Locke, Pistis, Porcu, Levy, Myers, Oetting, McGue, Abecasis and Iacono2014c ). Another compelling argument for the utility of endophenotypes is that, for a given gene known to be associated with a clinical disorder, an endophenotype may clarify, pinpoint, or suggest causal neurobehavioral mechanisms by which the gene exerts influence on the clinical outcome. Even if endophenotypes can elucidate neurobehavioral etiology in this way, it will require much larger samples than are routinely used for this purpose. The present study of about 4500 individuals, the largest of its kind to date, produced null results across a wide variety of analyses and genotyping technologies. None of the known schizophrenia variants showed a robust association with any of the endophenotypes. One suggestive single-variant result was for an intronic variant in the gene contactin 4 (CNTN4) associated with both antisaccade and P300, although this association did not survive multiple testing correction for all phenotypes. CNTN4 is a member of the contactin family of immunoglobulins that functions in neuronal network formation and plasticity. Deletion or mutation in CNTN4 has been found to be associated with autism spectrum disorders (Fernandez et al. Reference Fernandez, Morgan, Davis, Klin, Morris, Farhi, Lifton and State2004), and, of course, schizophrenia (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). The other variant was in an intron of the gene for prostaglandin-I2 synthase (PTGIS) which encodes for proteins in the cytochrome P450 enzyme superfamily.
The present study is also the largest of its kind to date to use whole genome sequencing to fine-map the 108 loci identified by the PGC genome-wide association study (GWAS) meta-analysis of schizophrenia. Our one significant gene-based association was with SBNO1, a gene that is highly conserved across species and has a function in regulating development (Egan et al. Reference Egan, St-Pierre and Leow1998), but no obvious biological link to neurobehavioral abnormalities relevant to affective modulated startle. In addition, this particular startle measure is perhaps the least strong endophenotype studied here (heritability = 0.01), casting further doubt on the utility of this finding. The lack of additional findings suggests that there are no rare or low-frequency variants of large effect on any of these 17 endophenotypes within the 108 PGC loci. This finding is perhaps not surprising given difficulties in finding rare variant associations with schizophrenia from other, larger exome sequencing studies (Purcell et al. Reference Purcell, Moran, Fromer, Ruderfer, Solovieff, Roussos, O'Dushlaine, Chambert, Bergen, Kähler, Duncan, Stahl, Genovese, Fernández, Collins, Komiyama, Choudhary, Magnusson, Banks, Shakir, Garimella, Fennell, DePristo, Grant, Haggarty, Gabriel, Scolnick, Lander, Hultman, Sullivan, McCarroll and Sklar2014; Singh et al. Reference Singh, Kurki, Curtis, Purcell, Crooks, McRae, Suvisaari, Chheda, Blackwood, Breen, Pietiläinen, Gerety, Ayub, Blyth, Cole, Collier, Coomber, Craddock, Daly, Danesh, DiForti, Foster, Freimer, Geschwind, Johnstone, Joss, Kirov, Körkkö, Kuismin, Holmans, Hultman, Iyegbe, Lönnqvist, Männikkö, McCarroll, McGuffin, McIntosh, McQuillin, Moilanen, Moore, Murray, Newbury-Ecob, Ouwehand, Paunio, Prigmore, Rees, Roberts, Sambrook, Sklar, Clair, Veijola, Walters, Williams, Sullivan, Hurles, O'Donovan, Palotie, Owen and Barrett2016).
The traditional definition of an endophenotype requires that cases show differences from controls but also that the genetic relatives of cases show differences from controls, in essence a requirement that the endophenotype have a genetic basis. One strength of the present study is that we were able to use a large number of known schizophrenia risk variants to rank individuals according to their genetic risk for schizophrenia in multiple ways. (1) For each of the PGC variants, we tested whether the schizophrenia risk allele was also associated with any of our endophenotypes. (2) We used a polygenic risk score to obtain a weighted sum of the number of risk alleles that each participant harbored. Ultimately, this procedure rank-ordered each individual in our sample according to their genetic risk for schizophrenia, and we correlated this rank ordering with each of the endophenotypes. While these results do not rule out the possibility that these endophenotypes share a genetic basis with schizophrenia, they do indicate that much larger sample sizes than those employed here will be required to detect such an association.
The utility of our 17 endophenotypes for characterizing genetic effects on neurobehavioral mechanisms in clinical samples remains to be seen, and is currently being pursued through work such as that being done by the Consortium on the Genetics of Schizophrenia (COGS). COGS is also evaluating a greater number of putative endophenotypes for study. COGS uses clinical probands to recruit families, thus increasing the variance of endophenotype expression in their sample and ensuring high levels of genetic risk in the unaffected individuals. They have reported on multiple heritable endophenotypes that differentiate cases and controls but it is currently unclear the extent to which they have discovered individual genetic loci robustly associated with any of their schizophrenia-related endophenotypes after one employs standard multiple testing and ancestry corrections (Greenwood et al. Reference Greenwood, Lazzeroni, Murray, Cadenhead, Calkins, Dobie, Green, Gur, Gur, Hardiman, Kelsoe, Leonard, Light, Nuechterlein, Olincy, Radant, Schork, Seidman, Siever, Silverman, Stone, Swerdlow, Tsuang, Tsuang, Turetsky, Freedman and Braff2011, Reference Greenwood, Lazzeroni, Calkins, Freedman, Green, Gur, Gur, Light, Nuechterlein, Olincy, Radant, Seidman, Siever, Silverman, Stone, Sugar, Swerdlow, Tsuang, Tsuang, Turetsky and Braff2016). Case–control samples like COGS and the PGC, because they employ an extreme sampling design, are more powerful than community samples for investigations of genetic association with disorders like schizophrenia. However, when the goal is to test known schizophrenia genetic loci for association with known schizophrenia-relevant endophenotypes, a case–control study of schizophrenia necessarily confounds diagnostic status and endophenotype level, resulting in the possibility of spurious gene–endophenotype correlations (Iacono et al. Reference Iacono, Malone and Vrieze2016). That is, if the gene is associated with schizophrenia, and cases are known to have higher levels of the endophenotype than controls, then the gene will also be associated with the endophenotype regardless of its causal role in schizophrenia. The same will be true for any phenotype (e.g. education level, marital status, etc.) that is systematically different in individuals with schizophrenia.
By current standards in complex trait genetics, our study was small. Our hypothesis was that the size of genetic effects for endophenotypes would overcome the relatively small sample size. Clearly, genetic effects on these endophenotypes are not large enough to be detected in this sample. Given the cost of collecting laboratory-based psychophysiological endophenotypes like those used in the present research it appears that endophenotypes are less efficient for gene discovery than simply assessing and genotyping more cases and controls. However, it remains to be seen the extent to which endophenotypes can illuminate the mechanistic pathways between genes and clinical outcomes. While antisaccade and P300 are convincing putative endophenotypes in schizophrenia and EEG and electrodermal measures have also received significant support, multiple promising putative endophenotypes for schizophrenia were not available in the present research (Allen et al. Reference Allen, Griss, Folley, Hawkins and Pearlson2009). Some of these endophenotypes may be associated with genetic loci of substantially larger effect sizes than any of our 17, although there is little evidence that psychophysiological or imaging-based endophenotypes are associated with larger effects than other complex traits and diseases (Iacono et al. Reference Iacono, Malone and Vrieze2016). Despite the fact that our study could not definitively test all candidate endophenotypes in schizophrenia, we suggest that our findings indicate that mechanism studies of the present kind will not be effective in community samples of the size investigated here, because none of the variants known to be associated with clinical outcome has any appreciable association with any of the endophenotypes. Samples of larger sizes with the endophenotypes measured here may yet uncover that these endophenotypes are indeed associated with schizophrenia risk loci, but the question remains whether such findings warrant the expensive in-person laboratory-based procedures required to measure these particular endophenotypes.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291716003184
Acknowledgements
This research was supported through National Institutes of Health grants AA023974, HG008983, DA037904, DA040177, DA041120, DA036216, DA005147, AA009367 and DA024417.
Declaration of Interest
None.
