Phenotypic Landscape of Schizophrenia-Associated Genes Defines Candidates and Their Shared Functions

Abstract

Genomic studies have identified hundreds of candidate genes near loci associated with risk for schizophrenia. To define candidates and their functions, we mutated zebrafish orthologs of 132 human schizophrenia-associated genes. We created a phenotype atlas consisting of whole-brain activity maps, brain structural differences, and profiles of behavioral abnormalities. Phenotypes were diverse but specific, including altered forebrain development and decreased prepulse inhibition. Exploration of these datasets identified promising candidates in more than 10 gene-rich regions, including the magnesium transporter cnnm2 and the translational repressor gigyf2, and revealed shared anatomical sites of activity differences, including the pallium, hypothalamus, and tectum. Single-cell RNA sequencing uncovered an essential role for the understudied transcription factor znf536 in the development of forebrain neurons implicated in social behavior and stress. This phenotypic landscape of schizophrenia-associated genes prioritizes more than 30 candidates for further study and provides hypotheses to bridge the divide between genetic association and biological mechanism.

Keywords: GWAS; behavior; forebrain; neurodevelopment; neuropsychiatric disorder; prepulse inhibition; schizophrenia; single-cell RNA-sequencing; whole-brain activity; zebrafish.

Figure 1.. Generation and analysis of 132 mutants for schizophrenia-associated genes

(A) Mutants were generated via Cas9 mutagenesis for genes found within and neighboring genomic loci linked to schizophrenia through genome-wide association. Manhattan plot image was adapted from (Schizophrenia Working Group of the Psychiatric Genomics, 2014). The resulting mutants were assessed for changes to brain morphology, brain activity, and behavior. Multiple zebrafish orthologs existed for 33 of the 132 genes (165 individual zebrafish genes), and both copies were mutated and assessed together. (B) Ricopili plot (https://data.broadinstitute.org/mpg/ricopili/) for multi-gene locus #3 of 108 from which five candidates (purple boxes) were selected for mutagenesis. (C) Ricopili plot for a gene (znf536) considered unambiguous because there are no other genes within 0.5 megabase on each side of the associated region, and the gene within the association region is brain-expressed. (D) Mutants made from 79 of 108 associated genomic loci. The locus rank reflects the statistical strength of the genetic association, with 1 being the most significantly associated. A region of 0.5–2 MB around each locus was analyzed, and genes outside of the region of linkage disequilibrium were selected for 19 of the 79 loci. Unambiguously associated genes are implicated strongly by previous literature, such as genes involved in glutamatergic neurotransmission, or are the only genes within or neighboring their locus (Table S1). (E) Mutations generated from Cas9 cleavage. Some mutant alleles consisted of several lesions if multiple gRNAs cleaved the genome independently, and all are included here. A range of mutations was recovered, tending to be either small (<15 bases) when a single gRNA cleaved or large (>100 bases) when a deletion spanned target sites of multiple gRNAs. (F) Protein sequence predicted to remain in mutants, based on sequence alignment identity. This analysis included both orthologs if the gene was duplicated (33 genes), for a total of 165 individual genes (163 in graph, as it does not include mir137 and one gene with unclear wild-type protein sequence length). The four mutants with >75% of the protein remaining did not have frameshifting mutations but did have phenotypes (Table S2), indicating that the protein function was disrupted. (G) Phenotypes in all 132 mutants based on analysis of brain activity signals and 71 behavioral assessments. See also Figure S1 and Figure S2 for the cutoffs for classifying which mutants have phenotypes. (H) Phenotypes in mutants for the 29 unambiguously associated genes (Table S1). (I) Phenotype dimensions affected in mutants for 132 genes from 79 schizophrenia-associated loci. Quantification of brain activity, brain structure, and behavioral differences for mutants designated as having a phenotype (Figure S1, Figure S2, Figure 1G, Figure 3B, Table S2) was scaled for comparison between the three measures, with the weakest phenotype designated as 0 and strongest as 1. Measures below the cutoff for phenotype designation (Figure S1, Figure S2) are displayed in black.

Figure 2.. Behavioral phenotypes in zebrafish mutants.

(A) Behavioral assay protocol. Frequency of movement, features of movement, and location preferences were calculated for each of 14 windows of the data (see STAR Methods for description of each time window). Stimulus responses were quantified from high-speed (285 frames-per-second) 1-second long movies. The exact stimulus paradigms are described in the STAR Methods. Behavioral testing was repeated for 73 of the 132 mutants, 46 of those designated as having a phenotype (Figure S1, Figure S2). (B). Example of altered frequency of movement for clcn3 −/− mutants compared to +/− sibling controls for entire protocol (time windows 1–14), including the approximate location of each time window. P-value = 0.005; N = 32 +/−, 33 −/−. (C) Example of altered features of movement for time window 14. P-value = 0.001; N = 36 +/−, 32 −/−. (D) Example of altered location preference for entire protocol (time windows 1 to 14). P-value = 0.001; N = 17 +/−, 18 −/−. (E) Example of a stimulus-driven high-speed response to a one second dark flash, and quantification of dark flash phenotypes across various modalities of the response. The response graph is an average of larvae in the mutant and control groups for events where a response was observed. Dark flash responses were analyzed for four hour-long blocks of dark flashes by assessing the first ten and last ten flashes each block, totaling eight separate statistical analyses. These eight p-values were combined using Fisher’s method to generate the heatmap of the five mutant phenotypes. Merged p-value for dark flash section (first ten flashes of block 4) for luzp2 = 0.02 (p-values are not calculated for stimulus response traces), with latency being the most significant contributing metric with p-value = 0.0007; N = 32 +/−, 35 −/−. (F) Summary of all 71 behavioral assays for all tested mutants. If a mutant was tested more than once, the lowest combined p-value is shown here. If multiple comparisons were made for a given gene, such as in the case of duplicated genes, the lowest p-value is shown. Significant results are reported for time windows (TW) 1–14 and the significance for the entire protocol (position 15 in heatmap) for each of the slower speed measurements. For high-speed stimulus responses, results presented in the heatmap are as follows: PPI (day early prepulse, e.g. responses in first 20 minutes of stimulus block, about 10–15 events; day all prepulse, e.g. all responses; night early prepulse; night all prepulse), Weak Tap (day early; day all; night early; night all), Strong Tap (day early; day all; night early; night all), Tap Habituation (day block 1; day block 2; day block 3; night block), Light flash (day; night), Dark Flash (block 1 start, e.g. first 10 flashes of 60 minute block; block 1 end, e.g. last 10 flashes of 60 minute block; block 2 start; block 2 end; block 3 start; block 3 end; block 4 start; block 4 end). Specific examples of raw metrics for every mutant are available on genepile.com/scz_gwas108. (G) The 46 mutants shown here had a total behavioral phenotype score greater than the cutoff (Figure S2) and were tested more than once. If the merged p-value (Figure S2) for the assay was significant (p < 0.05) between two independent repeats, the phenotype was counted and the lowest p-value of the two is displayed. Otherwise neither p-value was included. The significant phenotype had to be in the same time window or in the same stimulus assay to be counted (e.g., a significant phenotype observed on night 1 in one repeat and night 2 in another was not counted). Specific examples of raw metrics of the repeatable phenotypes for these 46 mutants are available on genepile.com/scz_gwas108. See also Figure S2. (H) The elfn1 mutants displayed increased responses to light stimulation. Merged p-value for day light flashes = 0.009 with 8/15 significant metrics, although the displayed frequency metric is itself non-significant; N = 31 +/−;+/−, 19 −/−;-/−. (I) The gpm6a mutants displayed a preference for the well center. P-value = 0.002; N = 19 +/−;+/−, 11 −/−;−/− (J) The akt3 mutants displayed reduced prepulse inhibition. The mutant response profile demonstrates the prepulse behavioral paradigm and response, where a weak tap (prepulse) is followed 300 milliseconds later by a strong tap. Merged p-value for the day early prepulse section = 0.018 with 10/36 significant metrics; N = 23 +/+ and +/−, 10 −/−. Plots of mutant compared to control groups in all panels represent mean ± s.e.m.. See also Figure S3. See also Figure S4.

Figure 3.. Whole-brain morphology phenotypes in zebrafish mutants.

(A) Location of several major regions in the zebrafish brain (Randlett et al., 2015). (B) Comparison between brain activity and morphology data for all mutants. The rora mutant (orange triangle) represents the smallest (5.95) structural change designated as a phenotype, shown in panel C. Small signals that were not symmetrical were considered to be noise, compared to small but symmetrical changes as observed in the rora mutant. (C) Examples of structural differences in mutants calculated using deformation-based morphometry and displayed as sum-of-slices projections (Z- and X-axes). Brain images represent the significant differences in signal between two groups (Randlett et al., 2015), most often homozygous mutants versus heterozygous and/or wild-type siblings (Table S2). The number of animals used in all imaging experiments is available in Table S2 and on stackjoint.com/zbrain (website naming conventions for datasets described in STAR Methods). Brain maps were averaged from two independent clutches of larvae if the experiment was repeated and irreproducible signals were eliminated (STAR Methods). Raw imaging data examples (maximum projections) are shown for the foxg1 mutants, demonstrating that the forebrain in homozygous (HOM) mutants is underdeveloped. A reduction in forebrain size of heterozygous (HET) foxg1 mutants when compared to wild-type (WT) siblings can be quantified with deformation-based morphometry (sum-of-slices projection), although it is not readily apparent in the raw projection. (D) Fluorescent RNA in situ images for four mutants with altered brain morphology.

Figure 4.. Whole-brain activity phenotypes in zebrafish mutants.

The number of animals used in all imaging experiments is available in Table S2 (average N is 20 for mutant group, 28 for control group) and on stackjoint.com/zbrain (website naming conventions in STAR Methods). Brain maps were averaged from two independent clutches of larvae if the experiment was repeated and irreproducible signals were eliminated (STAR Methods). The white numbers in the upper right corner of images connect pERK activity maps, both of single genes (panels A and B, and Figure S5) and gene averages (panel E), to the heatmaps and to each other. See also Figure S4. See also Figure S5. (A) Brain activity phenotypes for genes that have been strongly implicated in schizophrenia by previous studies (Table S1). Sum-of-slices projections of significant differences between mutant and control groups of zebrafish larvae (Randlett et al., 2015). The cacna1c mutant phenotype shown is for heterozygous larvae because the homozygous mutant is embryonic lethal (Stainier et al., 1996). (B) Brain activity phenotypes for four genes that have been minimally studied and have unknown functions. (C) Percent overlap between mutant brain activity phenotypes was calculated between each image by comparing each brain activity signal to signal in the same location in all other mutant images. These overlaps were then sorted with hierarchical clustering using average linkage. Genes grouped together by this clustering are labeled. The direction of the change in brain activity was disregarded to maximize identification of affected brain regions, and because the direction of the genetic perturbation in human patients is not clear for most genes. The numbers of image stacks that were compared to calculate significant differences in brain activity for each mutant are available in Table S2. (D) Contribution of each of the four major brain divisions to the overall brain activity phenotype. The signal in each region was divided by the whole brain signal. Prior to dividing the signal in each region by the total signal, the regions were scaled relative to each other based on their respective sizes (rhombencephalon = *1, diencephalon = *1.76, mesencephalon = *1.42, telencephalon = *4.36). The original whole brain signal was separately scaled across all mutants with a phenotype, indicating the relevance of the signal in each of the regions. Measures below the cutoff for phenotype designation (Figure S1, Figure S2) are displayed in black. (E). Averaged signal for mutants with similar brain activity maps. Examples of individual maps that are included in these averages are labeled in panels A and B, as well as in Figure S5. The blue arrow highlights retinal arborization field AF7. (F) Two mutants with brain activity signals that overlap with both the forebrain (group 2) and the tectum (group 4).

Figure 5.. Nominating candidates in multi-gene loci by phenotype.

The number of animals used in all imaging experiments is available in Table S2 and on stackjoint.com/zbrain. (A) Brain activity data (sum-of-slices projection) for lrrn3 and immp2l mutants. (B) Movement frequency for csmd3 mutants. P-value = 0.0001; N = 30 +/−, 29 −/−. (C) Mutants with tectum (green arrow) and retinal arborization field AF7 (blue arrow, Mesencephalon___Retinal_Arborization_Field_7_AF7 on stackjoint.com/zbrain) phenotypes (representative slices). Mutants with retinal arborization field AF7 signal also display signal changes in the same direction in a small subregion of hypothalamus (orange arrow, Diencephalon___Hypothalamus_Gad1b_Cluster_3_Sparse located within cyan Diencephalon___Intermediate_Hypothalamus; see regions on stackjoint.com/zbrain). See also kmt2e, znf804a, cacnb2b, and ambra1 (both areas decreased), and snap91, akt3b, and satb1 (both areas increased) on stackjoint.com/zbrain. (D) Prepulse inhibition phenotypes for five mutants. These mutant phenotypes are specific to the strong prepulse tap (Figure 2J) and do not represent a general increase in tap sensitivity to strong taps (right heatmap). Response features were calculated only on strong tap responses where the weak tap did not elicit movement. See also Figure S3. (E) Habituation phenotype of astn1 mutants. Response frequency to tap events occurring every two seconds is shown in left graph, and the magnitude of responses occurring during the habituation paradigm in the right graph. P-value for frequency metric = 0.0018; merged p-value for the day tap habituation 2 section = 0.015 with 20/47 significant metrics; N = 36 +/−, 18 −/−.

Figure 6.. Neurobiological roles of top candidates.

(A) Example images of larvae after consuming fluorescently labeled paramecia. (B) Quantification of feeding behavior in tcf4 mutants by measurement of paramecia consumed. (C) T-distributed Stochastic Neighbor Embedding (t-SNE) (Hinton and Maaten, 2008) visualization of wild-type single-cell clusters obtained by clustering of 6 dpf forebrain cells. Clusters with substantial differences in znf536 mutants are highlighted in orange, purple, and blue. Cluster counts in mutant and wild type are expressed as percent of the total cell number for each sample. See also Figure S6. (D) Dotplot (confusion matrix) showing the proportion of cells in the znf536 mutant forebrain that were classified to wild-type cluster labels. Each mutant forebrain type was assigned to a wild-type cluster label if > 13% of the trees in the Random Forest (RF) model (Breiman, 2001; Pandey et al., 2018) contributed to majority vote.