Gene expression elucidates functional impact of polygenic risk for schizophrenia

Abstract

Over 100 genetic loci harbor schizophrenia-associated variants, yet how these variants confer liability is uncertain. The CommonMind Consortium sequenced RNA from dorsolateral prefrontal cortex of people with schizophrenia (N = 258) and control subjects (N = 279), creating a resource of gene expression and its genetic regulation. Using this resource, ∼20% of schizophrenia loci have variants that could contribute to altered gene expression and liability. In five loci, only a single gene was involved: FURIN, TSNARE1, CNTN4, CLCN3 or SNAP91. Altering expression of FURIN, TSNARE1 or CNTN4 changed neurodevelopment in zebrafish; knockdown of FURIN in human neural progenitor cells yielded abnormal migration. Of 693 genes showing significant case-versus-control differential expression, their fold changes were ≤ 1.33, and an independent cohort yielded similar results. Gene co-expression implicates a network relevant for schizophrenia. Our findings show that schizophrenia is polygenic and highlight the utility of this resource for mechanistic interpretations of genetic liability for brain diseases.

Figure 1. Enrichment of cis-eQTLs in regulatory and other genomic elements

(a) Enrichments of cis-eQTLs compared to all eQTLs in sequence-defined elements according to the Ensembl annotations implemented in the ANNOVAR (version 2014-07-14) software. The bars illustrate the proportion of SNPs that belong to each category for significant cis-eQTLs (at FDR 5%) compared to all cis-SNPs (within 1 Mb from expressed genes). These categories are illustrated: exonic (fold change (FC) = 2.14); intronic (FC = 1.3); upstream (1 kb region upstream of transcription start site (TSS); FC = 1.48); downstream (1 kb region downstream of transcription end site (TES); FC = 1.52); UTR3 (3’ untranslated region; FC = 2.10); UTR5 (5’ untranslated region; FC = 2.35); splicing (within 2 bp of a splicing junction; FC = 2.51); ncRNA (transcripts without coding annotation in the gene definition, within either the exonic or intronic region; FC = 1.62 or 0.91, respectively); intergenic (FC = 0.69). (^) and (*) indicate significant (I_adjusted < 0.05) depletion or enrichment of cis-eQTLs compared to all cis-SNPs, respectively. (b) Distribution of cis-eQTL location relative to the gene. (c) Enrichment of “max-cis-eQTLs” (most associated eSNP per gene) within enhancer sequences across 98 human tissues and cell lines. Bars represent the Z score for the overlap of max-cis-eQTLs compared to 1,000 sets of random SNPs matched with respect to allele frequency, gene density, distance from the TSS, and linkage disequilibrium density. Brain (red) shows significantly higher enrichment for eQTLs compared to non-brain tissues and cell lines (P = 4.5 × 10⁻⁶) and the strongest enrichment is observed in DLPFC enhancers.

Figure 2. Overlap of GWAS for schizophrenia with eQTL in the DLPFC

(a) eQTL association profiles across two representative SCZ GWAS loci on chromosomes 15 and 4, respectively. SNP-level associations are plotted for the SCZ GWAS (gray), and cis-eQTL association profiles for genes with Sherlock P_corrected < 0.5 (or RTC > 0.9) are plotted in colors, with colors and Sherlock P values noted on top of the graphic (P = 4.07 × 10⁻⁷ and P = 4.07 × 10⁻⁷ for FURIN and CLCN3, respectively). For additional genes in the region with significant eQTL, the single eSNP with minimal eQTL P value (“max-eQTL”) is marked by a black point (corresponding genes names are located above the chromosome marker bar). Locations of regional protein-coding genes and non-coding RNAs without significant eQTL are annotated in gray. Vertical dotted lines mark recombination hotspot boundaries; horizontal dotted lines denote the significance thresholds for eQTL and GWAS, and the ceiling imposed for visualization purposes. Association betas (effect sizes) are plotted for SNP alleles associated with increased SCZ risk, in colors corresponding to genes as above. The red points illustrate the betas for the SCZ risk alleles on expression of the corresponding gene (FURIN and CLCN3, respectively), where values above the 0 line mark up-regulation (CLCN3) and below the line down-regulation (FURIN). (b) The association of expression of FURIN (N = 467, β = −0.071, P = 4.5 × 10⁻¹³) and CLCN3 (N = 467, β = 0.037, P = 1.6 × 10⁻⁹) with SCZ risk allele at the GWAS index SNP in the respective loci from (a), with shape corresponding to diagnosis.

Figure 3. Neuroanatomical phenotypes upon suppression or overexpression of genes at SCZ risk loci

(a) Head size phenotype after suppression of furin_a (3ng MO) or overexpression of TSNARE1, CNTN4, SNAP91 or CLCN3 (200ng). Representative head size images per treatment condition are shown, quantified area is depicted by the dashed white lines in the control image. (b) Quantification of head size phenotype in each treatment condition as compared to control embryos for furin MO (N_control = 76, N_{furin MO} = 66, P = 5.32 × 10⁻²⁰), TSNARE1 (N_control = 78, N_TSNARE1 = 64, P = 4.69 × 10⁻⁵), CNTN4 (N_control = 66, N_CNTN4 = 75, P = 0.018), SNAP91 (N_control = 114, N_SNAP91 = 106, p = 0.57), CLCN3 (N_control = 114, N_CLCN3 = 100, P = 0.40). (c) Representative images of PH3 staining assessing proliferation phenotypes. Dashed blue lines depict the area included in the quantification of cell counts. (d) Quantification of PH3-labeled cells with respect to each treatment condition for furin MO (N_control = 19, N_{furin MO} = 20, P = 7.56 × 10⁻¹⁷), TSNARE1 (N_control = 40, N_TSNARE1 = 40, P = 0.018), CNTN4 (N_control = 39, N_CNTN4 = 38, P = 0.0032), SNAP91 (N_control = 40, N_SNAP91 = 40, P = 0.25), CLCN3 (N_control = 40, N_CLCN3 = 40, P = 0.07). (e) Representative images of TUNEL staining per condition marking cells undergoing apoptosis. Area quantified is depicted within the dashed blue lines. (f) Cell counts of apoptotic cells in each treatment condition as compared to controls for furin MO (N_control = 33, N_{furin MO} = 39, P = 1.10 × 10⁻¹⁰), TSNARE1 (N_control = 33, N_TSNARE1 = 38, P = 9.44 × 10⁻⁶), CNTN4 (N_control = 33, N_CNTN4 = 35, P = 1.98 × 10⁻⁸). Error bars are s.e., * P < 0.05, ** P < 0.005, *** P < 0.0005; MO - morpholino. Scale bar = 100 µm. In all cases, t- tests were used to generate P values.

Figure 4. Decreasing FURIN expression in human NPCs perturbs neural migration

(a) FURIN expression reduction achieved by lentiviral (LV)-FURIN shRNA-PURO, relative to LV-non-hairpin-PURO control (C1: N = 6; P = 4.5 × 10⁻⁴; C2: N = 6, P = 6.2 × 10⁻⁹; C3: N = 5, P = 4.2 × 10⁻⁶). (b) Representative images of the hiPSC NPC neurosphere outgrowth assay after 48 hours of migration, following transduction with LV-FURIN shRNA-PURO and LV-non-hairpin-PURO control. The average distance between the radius of the inner neurosphere (dense aggregate of nuclei) and outer circumference of cells (white dashed line) was calculated. DAPI-stained nuclei (blue), scale bar 100 µm. (c) Across hiPSC NPCs generated from three controls (C1: N_vehicle = 42, N_shRNA-FURIN = 44, 1.16-fold decrease, P < 0.0017; C2: N_vehicle = 49, N_shRNA-FURIN = 53, 1.23-fold decrease, P < 3 × 10⁻⁶; C3: N_vehicle = 56, N_shRNA-FURIN = 63, 1.22-fold decrease, P < 2 × 10⁻⁶), average radial neurosphere migration following transduction with LV-FURIN shRNA-PURO (red bars) or LV-non-hairpin-PURO (gray bars). Error bars are s.e., *P < 0.05, **P < 0.01, ***P < 0.001. In all cases, a t-test was used to generate P values.

Figure 5. Differential expression between schizophrenia cases and controls in the DLPFC

(a) For the N = 693 genes differentially expressed at FDR ≤ 5%, bivariate clustering of individuals (columns) and genes (rows) depicts the case-control differences, as marked by the red-blue horizontal colorbar at top (‘Diagnosis’). An individual’s expression (converted to a z-score per gene) is red for above-average values, and green for below-average values; thus, the top cluster of the plot consists of genes up-regulated in cases versus controls (green in top left; red in top middle), and the bottom cluster of down-regulated genes (red in bottom left; green in bottom middle). In addition to the horizontal colorbar marking case-control status for each sample, additional colorbars denote brain bank (‘Institution’), gender, reported ancestry (‘Ethnicity’), age of death, and RNA quality (‘RIN’), where the latter two use a continuous-values color scale (with low, medium, and high as colored), relative to the range denoted on the figure. (b) Distribution of fold-change of differential expression for 693 differentially expressed genes. Case:control fold-changes for up-regulated genes are plotted in red (N = 332, positive values), and control:case fold-changes for down-regulated genes in green (N = 361, negative values). (c) Binned density scatter plot comparing the t-statistics for case versus control differential expression between the independent HBCC replication cohort assayed on microarrays and the CommonMind RNA-seq data; correlation between the statistics is 0.28 (P < 10⁻¹⁶). (d) For the 10 significantly differentially expressed genes with the largest fold changes (5 up- and 5 down-regulated), the 25 cases and 25 controls of normalized and adjusted gene expression in cases (red) versus controls (blue).

Figure 6. Co-expression network analysis in control DLPFC samples

(a) Control-derived modules were ranked by enrichment [estimated based on Fisher’s exact test (FET)] with differentially expressed genes; number of genes in each module is given in parentheses. Among the 4 modules with strongest overlap (marked in blue), only the M2c module genes are strongly enriched for multiples lines of prior genetic evidence: differential expression (FET: OR = 2.3, Bonferroni adjusted P = 1.9 × 10⁻¹²), SCZ GWAS loci (tested by INRICH: FE [fold-enrichment] = 1.36, P = 0.04), rare CNV (tested by INRICH: FE = 1.52, P = 0.051), and rare nonsynonymous variants (tested by PLINK/Seq and SMP: FE = 1.18, P = 2 × 10⁻⁴). The enrichment of each module with SCZ genetics, cell type-specific markers, neuronal proteome sets (proteins that are localized to the postsynaptic density of neurons), and fragile X mental retardation protein (FMRP) targets is depicted at right. As a control, note the lack of enrichment of M2c with common variants for Alzheimer’s disease (AD) and rheumatoid arthritis (RA). (b) Topological overlap matrix of the differentially connected M2c module in controls (upper right triangle) and SCZ cases (lower left triangle) in the CMC (left) and HBCC (right) cohorts. (c) Circle plot showing connection strengths for the top 50 hub genes of the M2c module, where node size corresponds to intramodular connectivity and nodes are ordered clockwise based on connectivity. Pie chart: SCZ susceptibility genes based on GWAS PGC2-SCZ (green), CNV (orange) or de novo (cyan) studies; Genes that belong in the NMDA (black) or mGluR5 (yellow) signalling pathway; Genes that are differentially expressed in schizophrenia vs. controls at FDR ≤ 5% (red).

Figure 7. Power to detect differential expression

Analysis of power to detect differential expression of a gene for case versus control subjects, where differential expression is expressed as expected log-fold change, the sample size is the total number of cases and controls to achieve significance (50:50 cases:controls), and the significance level for 80% power is 5 × 10⁻⁶. (a) For each gene in the differential expression analysis, we found the cis-eQTL with the smallest P value (see text for additional restrictions). Expected differential expression to achieve 80% power was computed for 10,094 gene-by- cis-eQTL associated pairs. (b) Increased resolution of (a) by limiting the range of differential expression. (c) Standardized log-fold change (80% power) obtained by dividing estimated log-fold change by its estimated standard deviation.