Genomic convergence analysis of schizophrenia: mRNA sequencing reveals altered synaptic vesicular transport in post-mortem cerebellum

Abstract

Schizophrenia (SCZ) is a common, disabling mental illness with high heritability but complex, poorly understood genetic etiology. As the first phase of a genomic convergence analysis of SCZ, we generated 16.7 billion nucleotides of short read, shotgun sequences of cDNA from post-mortem cerebellar cortices of 14 patients and six, matched controls. A rigorous analysis pipeline was developed for analysis of digital gene expression studies. Sequences aligned to approximately 33,200 transcripts in each sample, with average coverage of 450 reads per gene. Following adjustments for confounding clinical, sample and experimental sources of variation, 215 genes differed significantly in expression between cases and controls. Golgi apparatus, vesicular transport, membrane association, Zinc binding and regulation of transcription were over-represented among differentially expressed genes. Twenty three genes with altered expression and involvement in presynaptic vesicular transport, Golgi function and GABAergic neurotransmission define a unifying molecular hypothesis for dysfunction in cerebellar cortex in SCZ.

Figure 1. Characteristics of cDNA sequencing-by-synthesis with Illumina/Solexa instrument.

(A) Run-to-run comparisons of the number of reads aligned per reference transcript for 2 SBS runs of cerebellar cortex #41 (21.6 and 17.0 million reads, respectively). Runs used the same library but were performed on different sequencing instruments on different days. Coefficient of variation was 3.4%. (B) Histogram showing coverage along an “averaged” reference transcript for 1.2 Gb of cerebellar cortex #41 cDNA sequences. Coverage was calculated at 1% intervals along each transcript to which reads aligned. At each interval coverage was averaged across all transcripts and plotted. “Short transcripts” are all transcripts of ≤500 bp to which reads were aligned. “Long transcripts” are all transcripts ≥10 kb to which reads were aligned. Numbers in parentheses are the number of transcripts represented by each category. (C) Histogram showing greater than exponential decline in frequency with increasing SBS coverage on reference transcripts for cerebellar cortex #41 (1.2 Gb). Of 33,938 transcripts, 111 had more than 300× coverage. Maximum coverage was 2323×. Best fit power trendline is shown. (D) Histogram showing the number of reference transcripts covered as a function of the amount of sequence generated. Different levels of minimum average coverage were examined: at least one read aligned (black circles), 1× average coverage (red triangles), 2× average coverage (green “+”), 4× average coverage (dark blue “×”), and 8× average coverage (light blue squares). Data are from SBS of cerebellar cortex #41 (1.2 Gb).

Figure 2. Comparison of gene expression levels measured by two methods.

Comparison of gene expression levels detected by average shotgun, mRNA SBS coverage per transcript (normalized for transcript length) versus Affymetrix oligonucleotide microarray hybridization for cerebellar cortex sample #41 (38.7 million reads). Microarray images were scaled to an average hybridization intensity of 200, and the threshold for expression was an average hybridization intensity of ≥50.

Figure 3. Comparison of Mahalanobis Distances of Gene Expression by Array Hybridization and Sequence Read Frequencies.

14 SCZ samples are indicated by blue circles, and 6 control samples by red circles. The Y-axis shows Mahalanobis distances of log transformed gene expression values. The dotted blue line indicates the cutoff value for outliers. Panel A: Log10 transformed Affymetrix array hybridization signals. Panel B: Log10 transformed genome-aligned read frequencies. Panel C: Log10 transformed transcript-aligned read frequencies. Log10 transformed array hybridization values (A) had a wider distribution of distances than Log10 transformed sequence read frequencies (B,C). Without log transformation, distances were greater and several samples represented outliers (data not shown).

Figure 4. Overlayed kernel density estimates of Gene Expression by Array Hybridization and Sequence Read Frequencies.

14 SCZ samples are indicated by blue lines, and 6 control samples by red lines. The X-axis shows log transformed gene expression values while the Y-axis shows kernel densities. Panel A: Log10 transformed Affymetrix array hybridization signals. Panel B: Log10 transformed genome-aligned read frequencies. Panel C: Log10 transformed transcript-aligned read frequencies. Log10 transformed array hybridization values less than 1.69 (equivalent to a calibrated hybridization signal of 50) are considered noise. Without log transformation, samples showed greater variability in kernel densities and sequence read frequencies showed a near exponential decay (Fig. 1C, data not shown).

Figure 5. Pairwise sample correlations of log10-transformed, genome-aligned read frequencies, showing pairwise correlation coefficients.

Figure 6. Pairwise sample correlations of log10-transformed, array hybridization signals, showing pairwise correlation coefficients.

Figure 7. Unsupervised principal component analysis of log10 transformed array hybridization signals (A) and log10 transformed read frequencies (B).

Three dimensional plots of principal component analysis by Pearson product-moment correlation. 14 SCZ samples are indicated by blue circles, and 6 control samples by red circles.

Figure 8. Principal components of variance of log10 transformed array hybridization signals (A) and log10 transformed read frequencies (B).

Variance components decomposition of principal components (with Pearson correlation), with partitioning of variability in terms of known effects. Patient (Diagnosis, Cause of death [Death_cause], Age, Race, Medication), sample (Post-mortem interval [PMI], brain pH [pH_2], RNA integrity number [RIN], RNA isolation date [RNA_isol_date]), and experimental (Year Sequenced, Average Read Quality, Average Read length, Cluster Station, Illumina 1G instrument [SNPster], Library Creator, % reads aligning, number of reads) parameters were examined to quantify sources of variability in read-frequency- and array hybridization-based gene expression.

Figure 9. Volcano plot of analysis of variance of log10-transformed, transcript- (A) and genome-aligned read frequencies (B), showing genes with FDR-corrected differences in LSMeans.

The x-axis shows the magnitude of the difference between 14 cases (S) and 6 controls (C), while the Y-axis shows the −log10(p value) for those differences. The red dashed line indicates the significance cutoff for differences with a control (−log10(p value) 3.17 for genome alignments and 2.94 for transcript alignments).

Figure 10. Cartoon illustrating functions and/or synaptic locations of 23 proteins corresponding to genes with altered expression in SCZ.

15 genes were upregulated (green), whereas 8 were downregulated (red). Underlined genes had >30% change in expression. Two genes involved in transport from the endoplasmic reticulum to the Golgi (GOLM1 and GPSN2) were downregulated, ten involved in transport from the trans-Golgi network to the synaptic vesicle were upregulated (GOLGA1, SLC35A3, COG6, TRIP11, AP1G1, ARFGEF2, USO1, ROCK1, RAB9B and VPS35) and two were downregulated (STX10 and ARFRP1), two genes involved with synaptic vesicle exocytosis (EEA1 and SYT1) were upregulated and two were downregulated (SV2A and NCDN), one gene involved in receptor-mediated endocytosis was upregulated (AAK1) and one involved in retrograde transport back to the Golgi apparatus was downregulated (SNX17). In addition, three post-synaptic membrane genes showed altered expression: GABRA1 (upregulated), ZACN (downregulated) and CACNG2 (upregulated).

Figure 11. Novel alternative splicing of the PRODH locus.

A. Exons (green boxes) and introns (yellow lines) of the PRODH locus are shown. A 100 bp shaded box covers part of intron 13, intron14 and exon 14 and part of exon 15. Two sSNPs are illustrated by vertical blue hash marks. B. Alignments of SBS reads to the introns and exons within the 1000 bp shaded box are shown. Sequence reads are shown as arrows pointing in the direction of their orientation relative to the genomic reference. Yellow arrows represent sequence reads that map to more than one region on the genomic reference, and green arrows represent uniquely mapping reads. Three cDNA reads that omit intron 14 and contain exon 14 and 15 sequences are highlighted. Also highlighted are cDNA reads that map within intron 14. C. The sequence of the PRODH genomic region shown in B. Green highlights represent exonic nucleotides with aligned sequence reads and yellow represent intronic nucleotides with aligned sequence reads. Reads aligning uniquely to intron 14 (75 nts) indicate the existence of an alternative splice isoform that reads through this intron. Similarly, the 3′ region of intron 13 appears to be included in a novel splice isoform(s).