|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: GPS SFK. Performed the experiments: JM NAM IK IEL GM JJH LS LZ JCvV ADF RJL RWK. Analyzed the data: JM NAM IK IEL GM JJH JCvV ADF SL WB FDS SMV CFB MKM LCM RJL SKK MG VAG REH WC SM RDW NPB GPS SFK. Contributed reagents/materials/analysis tools: NAM GM ADF SL WB FDS SMV CFB MKM LCM JPU RWK RCR WC SM NPB GPS. Wrote the paper: JM SKK MG REH RDW NPB GPS SFK.
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.
Schizophrenia (SCZ [MIM #181500]) is a common, severe psychiatric disorder with a strong genetic component and complex inheritance –. Cytogenetic, linkage, positional cloning, candidate gene association and genome-wide studies have identified several credible SCZ risk genes. However, few have yet been replicated or translated into causal alleles, in vitro diagnostics or therapeutics –. In many studies of SCZ, genetic analysis has been impeded by phenotypic definition based upon multiple, subjectively ascertained, behavioral parameters that lack reference to biological mechanism –. In addition, more than 20 whole-genome linkage scans have demonstrated heterogeneity of linkage , suggesting the existence of genocopies (similar phenotypes that are determined by distinct risk loci) . Evidence exists that the genetic architecture of SCZ may be further obscured by allelic heterogeneity (additive genetic variance segregating in the population at causative loci), epistasis (different combinations of loci producing a phenotype in different pedigrees), pleiotropy (loci that affect more than one phenotype) and phenocopies (heterogeneous environmental factors that mimic allelic effects) . In addition, contributions of risk alleles to complex traits may not fit basal multiplicative and/or threshold models, and studies performed to date may not have had sufficient power, or appropriate theory, to assess non-linear (i.e., epistatic and genotype-by-environment) models. As a result, case-control association studies have identified numerous significantly associated susceptibility loci, but lack of replication among studies is widespread , , , –. Furthermore, the most validated loci were largely selected based on involvement in networks implicated in SCZ (such as dopaminergic and glutamatergic neurotransmission), introducing bias and limiting identification of novel risk factors.
Absence of a clear understanding of the molecular basis of SCZ imposes significant challenges to timely diagnosis and prognostic or therapeutic categorization , . Supplementation of diagnostic criteria with biomarkers that are causally related to SCZ or endophenotypes may allow definition of homogeneous subgroups that are predictive of progression and therapeutic response in individual patients  and would serve as a starting point for development of therapeutics directed at causal variants. An alternative approach for molecular dissection of SCZ is identification of altered gene expression in affected tissues. Because gene expression reflects both genetic and environmental influences, it may be particularly useful for identifying risk factors for a complex disorder such as SCZ, which is believed to have a multifactorial etiology . Two factors have hitherto limited the effectiveness of gene expression analysis in SCZ: Firstly, mRNA analyses in post-mortem brains in SCZ is challenging due to type I and type II errors resulting from variation in cause of death (affecting agonal gene expression), postmortem interval (affecting RNA quality), concurrent medication, substance abuse, age, sex, race and duration of illness . Secondly, in common with genome-wide association studies, gene expression comparisons employing available cohort sizes sail between the Scylla of many false-positives due to multiple comparisons and the Charybdis of insufficient power to detect true-positives following statistical correction . Recently, however, studies of mRNA expression in post mortem brains in SCZ that account for these variables have started to be reported , , .
An elegant, new approach to navigate Scylla and Charybdis and, thereby, accomplish molecular definition of SCZ may be genomic convergence analysis . Predicated on an implication of the central dogma of molecular biology, genomic convergence analysis posits that clinically relevant nucleotide variation should result in detectable cis- and trans-effects in messenger RNA (mRNA) that amalgamate into functional changes in networks and pathways. Importantly, genomic convergence analysis provides a strategy to collectively interpret and employ the massive, disease-related data sets produced by unbiased (i.e. non-hypothesis driven) linkage and expression studies. Indeed, integration of gene expression and genetic linkage data has shown promise in several neurologic disorders – and has started to be applied to SCZ , –.
mRNA sequencing with shotgun, massively parallel sequencing platforms has recently shown utility for measurement of transcript abundance, splice isoforms and allelic influence on gene expression –. mRNA abundance is determined by sequencing either 3′ end tags or random cDNA fragments (digital transcript expression, DTE), followed by read alignment to reference databases and calculation of aligned read frequencies. Potential advantages of DTE in comparison to array hybridization include: single molecule sensitivity (corresponding to approximately 1 mRNA molecule per 30 cells; Hayashizake, personal communication); absolute, rather than relative, measurement of transcript abundance; sequence verification for each measurement; comprehensive detection of both known and novel, unannotated transcripts and isoforms; applicability to any eukaryotic species; very little technical imprecision; absence of interference from abundant transcripts (e.g. globin); and extensibility to concomitant measurement of non-coding RNA and to detection of nucleotide and structural variation.
As the first stage of a genomic convergence analysis of SCZ, we describe shotgun mRNA sequencing of an affected tissue (post-mortem cerebellar cortex) of patients and controls, together with a rigorous and systematic approach to DTE. The approach presented is conservative due to application of statistical and bioinformatic methods that substantially reduce type I error rates by using statistical significance criteria rather than fold-change values and by incorporation of potential confounding variables, such as psychotropic medication and cause of death. It is also more comprehensive than prior microarray studies of SCZ. The stringency of this approach, in combination with genome-wide genotypes from the same individuals, is anticipated to enable a genomic convergence analysis that advances understanding of the molecular basis of SCZ.
16.7 billion nucleotides of shotgun, full length cDNA sequence data were generated using Illumina Genome Analyzer platforms with sequencing-by-synthesis (SBS) chemistry from 20 mRNA samples , , . mRNA samples were isolated post-mortem from the lateral hemispheres of the cerebellar cortices of 14 patients with SCZ and 6 control individuals ,  (Table 1). Unrelated subjects were chosen to facilitate sampling of genetic heterogeneity , . Cases and controls were approximately matched for age (cases, 45.2±11.8 years; controls 41.3±9.2 years), sex (all male), race, post-mortem interval (cases, 12.2±5.0 hours; controls 17.7±3.3 hours), cause of death, autopsy brain pH (cases, 6.54±0.19; controls 6.46±0.10) and RNA integrity number (cases, 8.06±0.53; controls, 7.87±0.41) (Table 1). 12.5–38.7 million, high quality sequences of length 32–36 bp were generated per sample (Table 2). Sequences were aligned to the human genome and RefSeq transcript databases using the algorithm GMAP, which allowed <2 (≤6%) mismatches . There was little intra-sample variability in the number of sequences aligned to each locus from run-to-run or instrument-to-instrument (Fig. 1A; all source coefficient of variation 3.4%). 43.5±6.7% of sequences aligned to a transcript and 69.4±9.6% to the genome, evidence that annotation of mRNA isoforms in Homo sapiens is incomplete  (Table 2). 91% of alignments were unique (Table 2). Reads aligning to more than one location contained repetitive, paralogous, polymorphic or low complexity sequences and primarily mapped to untranslated regions or highly polymorphic gene families, such as major histocompatibility genes . Unmapped sequences did not align to mitochondrial or 1879 viral genomes, offering negative evidence of chronic viral etiology for SCZ in these patients.
In order to further investigate differences in proportions of sequences aligning to different reference databases, an additional 2.2 billion nucleotides of shotgun, full length, paired cDNA sequence data were generated from mRNA sample 3S using an Illumina Genome Analyzer II platform and aligned to three reference databases (the human genome, RefSeq transcript and UniGene transcript databases): 46.1% of reads aligned to all three databases, representing known exons. 31.4% of reads aligned to the genome alone, representing novel exons, splice isoforms  or transcripts . 5.5% mapped to RefSeq transcript alone, and 6.5% to UniGene transcript alone, representing sequences that span known exon junctions . 10.6% of reads did not align to any of these reference sets, representing reads that span boundaries of novel exons or splice isoforms, novel genomic sequence, or poor quality sequence.
The number of transcripts detected in cerebellar cortex mRNA differed little between samples (33,200±1,000; Table 2), corresponding to 85±3% of RefSeq transcript entries. While 12.5 million sequences per sample was sufficient to reach a plateau in the number of transcripts detected, deeper sequence generation resulted in linear increase in average depth of coverage (Fig. 1D). As anticipated with hexamer-primed cDNA synthesis, the distribution of sequence alignments along transcripts appeared random. Transcript coverage showed a moderate 3′ bias, particularly among long transcripts, attributable to only very slight levels of degradation in these samples together with two rounds of poly-A+ RNA selection (Fig. 1B). Decreased coverage was observed at 5′ and 3′ termini due to edge effects of random priming very close (<30 nts) to the ends of the mRNA. No compositional bias was detected. The average depth of coverage achieved was 8.0-fold (Table 2). The distribution of transcript abundance in cerebellar cortex was interesting: 76–81% of expressed transcripts had an abundance of at least 1 read per million. However, as abundance increased, the number of transcripts declined greater than exponentially (Fig. 1C, r2=0.97). Thus, only 0.3–0.5% of expressed transcripts had an abundance of ≥1000 reads per million.
Gene expression was evaluated in mRNA samples from the 20 cerebellar cortices using both oligonucleotide array hybridization, the current standard, and tag frequencies from the Illumina cDNA sequencing assay (aligned reads per million). The dynamic range of read frequencies was two orders of magnitude greater than array hybridization (log10 dynamic range DTE–4.37; log10 dynamic range array hybridization–2.33; Fig. 2). While 41% of arrayed oligonucleotide probes had hybridization intensities greater than the conventional signal: noise threshold (Fig. 2), 85±3% transcripts had aligned reads. Array hybridization signals are transformed and often normalized prior to evaluation of inter-sample differences. For read frequencies, log transformation, but not normalization, improved overlaid kernel density estimates, univariate distribution results, and Mahalanobis distances (Fig. 3 and and4).4). Therefore, log transformed values were used in all subsequent expression analyses. Aligned read frequencies, adjusted for transcript length, correlated weakly with oligonucleotide array hybridization signals (r2=0.35; Fig. 2). Read frequencies exhibited much higher correlation coefficients in pair wise sample comparisons (r2=0.93–0.99) than array hybridization (r2=0.83–0.88; Fig. 5 and and6).6). Furthermore, pair wise sample correlations of genome- and transcript-aligned read frequencies had very similar correlation coefficients (data not shown). SCZ patients could readily be distinguished from controls by unsupervised principal component analysis (by Pearson product-moment correlation, Fig. 7B) or Ward hierarchical clustering of Pearson product-moment correlations of read frequencies (data not shown). In contrast, patients and controls could not be separated with these methods on the basis of array hybridization signals (Fig. 7A and data not shown). Log transformed genome- and transcript-aligned read frequencies showed identical hierarchical clustering and Pearson product-moment correlations of samples (data not shown).
Patient, sample, and experimental parameters were examined to quantify sources of variability in read-frequency-based DTE. Decomposition of principal components of variance showed that DTE attributed a greater proportion of the total variance to diagnosis (SCZ versus control, 45.3%) than array hybridization (14.1%) (Fig 8). The largest sources of variability of array hybridization results were mRNA quality metrics (Fig. 8A), whereas DTE results were only marginally influenced by mRNA quality. In contrast, large sources of variability in DTE results were year of sequence generation, sequencing instrument-to-instrument variation and cause of death (Fig. 8B). Variance component categories were identical in genome- and transcript-aligned read frequencies (data not shown).
Differences between SCZ patients and controls in gene expression in cerebellar cortex were identified with analysis of variance with diagnosis as the discriminatory effect and the major non-diagnosis components of variance as fixed effects (brain pH, post-mortem interval, RNA integrity number and RNA isolation date for array hybridization, and cause of death, post-mortem interval, sequencing instrument and year sequenced for read frequencies). Following FDR correction, no differences were identified between SCZ patients and controls by array hybridization, due to the magnitude of non-diagnosis components of variance. In contrast, 88 genes exhibited FDR-corrected differences in cerebellar expression between SCZ patients and controls in genome-aligned read frequencies and 152 genes differed significantly in transcript-aligned reads (Fig. 9; Table 3). Between these two sets of genes, 25 were identified in common to both genome- and transcript-aligned reads. 24 of these 25 genes showed congruent direction of change in genome- and transcript-aligned reads. 95% of genes exhibiting significant change in expression with only one alignment had congruent but non-significant change in the second alignment or were absent in one of the reference sets. The correlation coefficient of the case-control ratio of the 215 differentially expressed genes between genome- and transcript-aligned reads was 0.73 and between transcript-aligned reads and array hybridization was 0.46 (data not shown). A majority of genes with disparity in the magnitude of gene expression change between alignments had greater read counts for genome- than transcript-alignments, which, in some cases, appeared to be due to unannotated exons and novel splice isoforms (see below).
Genes with the largest increase in expression in SCZ samples were the ISL2 transcription factor, LIM/homeodomain, >800% in genome and transcript alignments; AMME complex, gene 1, 580% in transcript alignments; 3-hydroxy-3-methylglutaryl-Coenzyme A synthase, 372% in transcript alignments; kelch domain containing 1, 360% in transcript alignments; interphotoreceptor matrix proteoglycan 2, 321% in transcript alignments; and NK2 transcription factor related, locus 3, 320% in genome alignments. Genes with the largest decrease in expression in SCZ samples were the solute carrier family 25 (ornithine transporter), member 2, 59% in transcript alignments and LOC126170, similar to peptidylprolyl isomerase A isoform 1, 73% in transcript alignments.
An SCZ candidate gene showing differential expression by both measures was the gamma-aminobutyric acid (GABA)-mediated neuroinhibitory receptor, GABRA1 (increased expression in SCZ by 52% [transcript alignment] and 46% [genome alignment], Table 3). Eleven other genes involved in GABA neurotransmission showed non-significant increases in expression: GABRA2 (59%), GABRA4 (72%), GABRR1 (69%), GABRB1 (168%), GABRG2 (42%), GABRE (35%), GABRB2 (33%), GABBR2 (27%), GABARAPL2 (35%), GABRA6 (14%) and GABRB3 (8%). Four GABAergic genes showed non-significant decreases in expression: GABBR1 (22%), GABRG3 (54%), GABRR2 (26%) and SLC6A1 (GAT-1; 16%).
GO annotation of genes with significantly altered expression revealed over-representation of membrane associated genes (Χ2 test, p<0.01), genes involved in Zinc binding or transport (p<0.02), regulation of transcription (p<e−6), Golgi apparatus (p<e−6) and vesicle-mediated transport (p<0.01). 13 of 20 genes involved in vesicular transport or Golgi apparatus were up-regulated.
Further annotation of these 20 genes revealed up-regulation of:
Down-regulated transcripts corresponded to:
In addition, three post-synaptic membrane genes showed altered expression: GABRA1 (see above), ZACN (ligand-gated ion channel, zinc activated 1; down-regulated 21% in transcript alignments and 5% in genome alignments) and CACNG2 (calcium channel, voltage-dependent, gamma subunit 2; up-regulated 169% in transcript alignments and 10% in genome alignments) (Fig 10, Table 3).
10,022 genes exhibited ≥20% more reads aligned to genomic loci than to the corresponding, annotated transcripts, evidence for novel, unannotated exons. Putative use of alternative 5′ and 3′ terminal exons, novel internal exons and read-through of introns were all detected . A SCZ susceptibility locus that exhibited intron-read-through in cerebellar mRNA samples was proline oxidase (PRODH, EC 18.104.22.168, SCZ4): 421 cerebellar cortex 41 reads aligned to the sole annotated PRODH transcript, while 838 aligned to the genomic locus. 33 of the latter mapped to intron 14 (nucleotides 17280876–17280950 on Chromosome 22), resulting in insertion of 25 new, in-frame amino acids adjacent to the N-terminus (Fig. 11). In addition to reads mapping within this intron, reads showing splicing out of the intron were also observed (Fig. 11). A further 112 reads aligned to 660 nucleotides of intron 13, generating an alternative N-terminus and premature stop codon (Fig. 11). Presence of these intronic PRODH sequences was observed in all cerebellum mRNA samples and is supported by independently isolated Sanger EST sequences (Genbank accession numbers CN362766, AA300535, AA322439, CD671570, CN362770).
As noted above, the presence of unannotated exons and novel splice isoforms appeared to explain disparity in the magnitude of gene expression change between cases and controls when using genome- and transcript-alignments. Dedicator of cytokinesis 3 (DOCK3), for example, was significantly up-regulated in SCZ using genome alignments (225 reads/million in cases versus 200 reads per million in controls, −log10(p-value)=3.82) but non-significantly decreased using transcript alignments (138 reads/million in cases versus 147 reads per million in controls). Inspection of genomic read alignments revealed the existence of several, putative, unannotated exons. For example 26 reads in SCZ case 5S and 17 reads in SCZ case 5 aligned uniquely to DOCK3 intronic sequences corresponding to nucleotides 5,130,7158–51307328 on Chromosome 3 (data not shown).
As the first stage of a genomic convergence analysis of SCZ, we sought gene expression changes between post-mortem cerebellar cortices of 14 patients and 6 controls using shotgun, clonal mRNA sequencing, together with a rigorous and systematic analysis approach. Cerebellum was chosen for study since it has been shown to be affected in SCZ: Patients exhibit impaired sensory-motor integration, reduced cerebellar volume and grey matter, increased cerebellar blood flow and glucose uptake, both at rest and following cognitive challenge, greatest approximate conditional likelihood score in cerebellar vermis, and gene expression changes in the lateral hemispheres of the cerebellar cortex consistent with decreased GABAergic and increased glutamatergic neurotransmission , , , –. In addition, the simple organization and lack of perturbation by antipsychotic medication of cerebellar cortical circuits relative to other brain regions affected by SCZ was appealing for the current study. Specifically, granule cells, which are important in processing sensory information, are by far the most abundant cerebellar glutamatergic neuron, comprising the majority of neurons in the brain .
Cases and controls were all male and were matched reasonably well for age, race, cause of death and RNA quality metrics. Since only cases had received psychotropic medication, this variable could not be matched. There is currently considerable excitement about stratification of patients with SCZ on the basis of clinical or neuroimaging endophenotypes . Unfortunately, that was not possible in this series.
Approximately 22 million random, cDNA sequence tags were generated from each of 20 samples, providing the deepest mRNA coverage reported to date and providing an unparalleled assessment of the transcriptional complexity of the human cerebellar cortex. Sequences did not align to 1879 viral genomes, offering negative evidence of chronic viral etiology for SCZ in these patients. 85% of annotated transcripts were expressed (based on molecular counting, Table 2), of which 63% were detected at a level of at least 1-fold coverage (equivalent to 56 reads per million or 1.1 copies per million). This pervasive transcription accords with other recent studies , , , . In contrast to 3′ end tag sequencing, random-primed mRNA sequencing provided good coverage of coding domains (Fig. 1B), allowing assessment of exon utilization in cerebellar cortex. Alternative splicing is thought to generate more than 5 transcripts per human locus, most of which have not yet been annotated . 85% of loci expressed at a level of at least 1-fold coverage showed evidence for utilization of unannotated exons, based on excess genomic read alignments. Proline oxidase, for example, is a very well studied gene that appeared to utilize hitherto unannotated exons in human cerebellum (Fig. 11). For this reason, sequences were aligned both to a reference transcript database (for estimation of expression of well-annotated transcripts) and to the human genome (for estimation of total expression of all exons, whether annotated or not). With additional improvements in library preparation and analysis software it should be possible to measure digital gene expression on an exon-by-exon basis in a strand-specific manner –, . The number of transcripts with aligned reads declined greater than exponentially as a function of level of expression (Fig 1B). As a result, 29 million aligned reads provided an impressive linear dynamic range (4.37×log10, Fig. 2). Despite lower read lengths and quality scores in this study (Table 2) than currently obtained with the Illumina GA II instrument, only 9% of reads mapped non-uniquely. Finally, technical imprecision of DTE measurements was paltry (CV 3.4%). In summary, our technical assessment indicated deep mRNA sequencing to be a powerful tool for digital gene expression and empiric annotation of transcribed elements in genomes, in accord with a rapidly growing literature –, .
Gratifyingly, digital transcript expression values, expressed as aligned read frequencies, were amenable to very similar data transformation, quality control, pattern discovery, row-by-row modeling and annotation analyses as array hybridization datasets. This enabled direct comparison of digital transcript expression and array hybridization data metrics (Figs. 4–9).9). Like array hybridization signals, aligned read frequencies benefited from log transformation as evidenced by improved Mahalanobis distances, overlaid kernel density estimates, univariate distribution results, and unsupervised principal component analysis. Aligned read frequencies had greater correlation coefficients in pair wise sample comparisons than array hybridization, and much greater ability to distinguish SCZ cases from controls, as evaluated by unsupervised principal component analysis (PCA), Ward hierarchical clustering of Pearson product-moment correlations of read frequencies, and the magnitude of the component of variance attributable to diagnosis. DTE was clearly better than array hybridization for PCA-based segregation of cases and controls. Read frequencies based upon alignment to a reference transcript database and to the human genome were almost indistinguishable by these metrics. Partitioning of variance components of clinical, sample and experimental metadata revealed the largest components of array hybridization variability to be metrics of mRNA quality, in accord with recent studies of post-mortem brain tissue . DTE was much less influenced by RNA quality. Non-diagnosis DTE variability included year of sequence generation, sequencing instrument variation and cause of death (Fig. 8). This is not surprising since substantial improvements occurred in sequencing instrument specifications, reagents and base-calling software during the project. Given continued, rapid evolution of generation II sequencing technologies, it is important to keep instrument specification, reagents and base-calling software relatively stable during mRNA sequencing projects to minimize non-hypothesis-related variability. Importantly, variance decomposition allowed incorporation of large clinical, sample and experimental effects as fixed effects in analyses of variance, minimizing detection of gene expression changes associated with confounding variables. Incorporation of fixed effects in array hybridization, but not DTE, analyses of variance resulted in no expression changes meeting significance cutoffs. Coefficients of correlation between aligned read frequencies, adjusted for transcript length, and array hybridization signals were 0.35 for all values and 0.46 for fold-change in genes with significant differences in expression. These values are somewhat lower than similar comparisons between array hybridization platforms and are probably reflective of the exacting RNA source . It should be noted that aligned read frequencies correlate very well with quantitative PCR results , . In conclusion, mRNA sequencing appears to be superior to array hybridization for gene expression analysis and is amenable to standard statistical procedures for quality assessment and gene expression analysis.
Expression of 215 genes differed significantly between cases and controls in post-mortem cerebellar cortex by analysis of variance (Table 3). Of these, 88 were identified in genome-aligned reads and 152 in transcript-aligned reads. Only 25 genes were common to both alignments, but 96% of the 215 had congruent direction of change and correlation between genome- and transcript-aligned fold-change was 0.73. Major contributors to disparity between alignments appeared to be detection of unannotated exons in genome alignments and forced mis-alignments in incomplete transcript datasets. The latter may largely be avoided by removal of non-unique alignments from calculations. As RefSeq transcript becomes more complete, this disparity should decrease. In the interim, it is wise to align both to a well-annotated transcript dataset and to the genome.
17 of the 215 genes with altered cerebellar expression have previously shown changes in SCZ in the dorsolateral prefrontal cortex or superior temporal gyrus (APBA2, BTG1, CACNG2, CAP2, DKFZP434A0131, GABRA1, GOLGA1, HSPBP1, KIAA0256, KPNA1, RPL3, SLC35A3, SLC39A7, SRRM2, TCF4, ZNF148 and ZNF195). Despite substantial differences in gene expression in cerebellar cortex and cerebral cortex , 13 genes were congruent with the direction of change reported previously. Three genes involved in GABAergic neurotransmission (GABRA6, GABRB3 and SLC6A1) showed changes in cerebellar cortex in SCZ that agreed with a previous report . In the context of cerebellar cortical function in SCZ, alteration in expression of genes may be causal, consequential, compensatory, or the result of confounding factors. Because the substantial heritability of schizophrenia appears to be polygenic –, expression differences might reflect SCZ-associated nucleotide variants that alter expression in cis (e.g. eSNPs). Indeed, seven of these genes have previously shown association with SCZ in genetic association studies (CACNG2, GABRA1, GPSN2, HIRA, PSAP, RANBP5 and TCF4) . Alternatively, expression differences may represent compensatory changes in pathways and networks. GO annotation of genes with altered expression revealed over-representation of membrane association, Zinc binding or transport, regulation of transcription, Golgi apparatus and vesicle-mediated transport.
Most striking were 23 genes involved in presynaptic vesicular transport / Golgi apparatus or post-synaptic neurotransmission, 15 of which were up-regulated and 8 down-regulated. Up-regulated genes included nine involved in transport from the trans-Golgi network to the synaptic vesicle (golgi autoantigen A1, UDP-N-acetylglucosamine transporter (SLC35A3), component of oligomeric golgi complex 6, thyroid hormone receptor interactor 11, adaptor-related protein complex 1G1, ADP-ribosylation factor guanine nucleotide-exchange factor 2, vesicle docking protein p115, Rho-associated, coiled-coil containing protein kinase 1 (ROCK1) and RAB9B), one involved in synaptic vesicle exocytosis (synaptotagmin I), one involved in clathrin-mediated endocytosis (AP2 associated kinase 1 (AAK1)), one involved in homotypic fusion of early endosomes (early endosome antigen 1 (EEA1), and one involved in retrograde transport to the Golgi apparatus (vacuolar protein sorting 35 (VPS35)). Down-regulated genes included two involved in transport from the trans-Golgi network to the synaptic vesicle (syntaxin 10 and ADP-ribosylation factor related protein 1), two involved with synaptic vesicle exocytosis (synaptic vesicle glycoprotein 2A and neurochondrin), two involved in transport from the endoplasmic reticulum to the Golgi (golgi phosphoprotein 2 and synaptic glycoprotein 2) and one involved in retrograde transport to the Golgi apparatus (sorting nexin 17) (Fig. 10, Table 3). Up-regulated post-synaptic membrane genes included GABRA1 (see above) and CACNG2 (calcium channel, voltage-dependent, gamma subunit 2). Several other genes involved in GABAergic neurotransmission showed non-significant changes in expression in cerebellar cortex in SCZ, in agreement with previous reports . One down-regulated post-synaptic membrane gene was ligand-gated ion channel, zinc activated 1. Most of the up-regulated genes, but none of the down-regulated genes in this set, had change in expression of greater than 30%. The largest alterations were CACNG2 (+169%), EEA1 (+88%), ROCK1 (+85%), SLC35A3 (+69%), AAK1 (64%) and GABRA1 and VPS35 (both +52%). Several previous studies have shown decreased expression of genes involved in presynaptic vesicular transport in the dorsolateral prefrontal cortex and hippocampus in SCZ –. While the molecular determinants of synaptic vesicular transport have been intensively studied, it is not yet possible to integrate these transcriptional changes into a single biological outcome. Isolated elevation of AAK1, for example, should decrease neurotransmitter endocytosis , while elevated EEA1 should increase early endosomal fusion . However, activity of most of these proteins is dependent either upon post-translational modification or presence of interacting proteins. Furthermore, there is likely to be heterogeneous expression of these proteins among cerebellar cortical neurons. Nevertheless, it appears clear that presynaptic vesicular transport, Golgi function and GABAergic neurotransmission are perturbed in cerebellar cortex in SCZ, in agreement with previous analyses of gene expression in post-mortem brain in SCZ , ,  and genetic associations of nucleotide variants in synaptic vesicular transport genes with SCZ (namely CACNG2, COG2, STX1A, SNAP29, MUTED, NRG1, PLDN, NRG3, HOMER3, RTN4R, NRG2, RTN4, CPLX2, BLOC1S1, SYNJ1, HOMER1, CLINT1, SYP, BLOC1S2, SYN3, SYT11, SYNGR1, SNAPAP, STX7, BLOC1S3, DISC1, SNAP25, DLG4, CPLX1, CNO, SYN2, DTNBP1 and DAOA) . Replication of alterations in presynaptic vesicular transport and Golgi function, particularly as affecting GABAergic neurotransmission, in additional cases or other brain regions–such as prefrontal cortex–would substantiate novel molecular mechanisms underpinning SCZ and establish novel targets for therapeutic intervention.
The current study reports the first phase of a genomic convergence analysis of SCZ, which is intended to integrate linkage and expression studies. In a subsequent manuscript, we will present results of expressed nucleotide variants that differ in allele frequency between cases and controls, and integrate nucleotide variant and expression analyses. In addition, a recent study has shown substantial variation in alternative splice isoform expression and alternative polyadenylation in cerebellar cortex between normal individuals , and it will be interesting to ascertain whether specific examples of alternative isoform expression show association with SCZ. Confirmation of causal components of the molecular basis of SCZ is anticipated to have significant impact on clinical practice, particularly with respect to timely diagnosis and prognostic or therapeutic categorization , .
was as previously described . Anonymized cerebellar samples were provided from the Maryland Brain Collection with permission from the Maryland Brain Collection Steering Committee. Permission for the study of brain tissues was provided by families post-mortem in accordance with the guidelines of the Uniform Anatomical Gift Act. Fourteen samples were from patients with a diagnosis of SCZ according to DSM-IV criteria and six were controls. Cortical areas corresponding to crus I/VIIa of cerebellar hemispheres were dissected at −20°C and frozen at −80°C, as described . The average pH of the samples was 6.6±0.17. Cerebellum specimens were obtained according to NIH guidelines for confidentially and privacy. Genomic DNA and total RNA were isolated from samples using standard techniques (Qiagen, Valencia, CA).
Cerebellum samples were sequenced using Illumina Genome Analyzers with modifications for mRNA samples , , , , : Following quality assessment using a Bioanalyzer 2100 (Agilent Inc., Santa Clara, CA; Table 1), poly A+ RNA was isolated from 5–10 µg total RNA by two rounds of oligo-dT selection (Invitrogen Inc., Santa Clara, CA). mRNA was annealed to high concentrations of random hexamers and reverse transcribed. Following second strand synthesis, end repair, and A-tailing, adapters complementary to sequencing primers were ligated to cDNA fragment ends. Resultant cDNA libraries were size fractionated on agarose gels, 200 bp fragments excised, and amplified by 15 cycles of polymerase chain reaction. Following quality assessment using a Bioanalyzer 2100, single-stranded cDNA-adapter fragments were randomly annealed to the surface of a flow cell in a cluster station (Illumina Inc., San Diego, CA) via primers complementary to the adapters and incubated under conditions fostering annealing of the ends of cDNA-adapter fragments to adjacent complementary primers. Primers with annealed cDNA-adapter fragments were extended with DNA polymerase and unlabelled dNTPs in a solid-phase “bridge amplification”. Resultant double-stranded products were denatured, yielding 2 single stranded fragments. The latter steps were repeated for 35 cycles, generating ~40 million clusters of clonal amplicons. Subsequently, 32–36 cycles of sequencing-by-synthesis chemistry were performed in Illumina Genome Analyzer instruments with 4 dNTPs featuring cleavable dyes and reversible terminators. Following each base extension, 4 images (one for each nucleotide) are taken upon laser excitation. Incorporation of the next base occurred after removal of the blocked 3′ terminus and fluorescent tag of the previously incorporated nucleotide. Sequences were retained if of high quality (defined as those passing the default quality filtering parameters used in the Illumina GA Pipeline GERALD stage, i.e. clusters with intensities greater than 0.6-times the average of the highest and the sum of the two highest intensities for the first 12 cycles). Furthermore, sequencing runs with poor average nucleotide quality score graphs (associated with a minority of reads passing the default quality filtering parameters) were discarded.
High-quality reads were aligned to the human genome, Build 36.2, RefSeq transcript database, Release 22 , Unigene Hs, build 215  using the algorithm GMAP  and the software system Alpheus® , with adjustments for short SBS reads (oligomer overlap interval=3 nt, identities ≥34/36 or 94%). A read was denoted ‘aligned’ to a locus if its sequence alignment to the genomic reference sequence (NCBI build 36.2) fell within the boundaries of the locus coordinates on the chromosome. Locus boundaries on the genome were defined by NCBI annotations, as reported through the Nucleotide database. Annotation data was downloaded on 2/2/2007 using the Batch Entrez query tools (http://www.ncbi.nlm.nih.gov/sites/batchentrez?dbNucleotide). Only the highest scoring alignment(s) was retained. Reads with a single best alignment or with equally good alignments to alternative transcripts of the same locus were considered uniquely aligned. Aligned read frequencies (per million aligned reads) were calculated for each sample and locus using Alpheus® .
Array-based gene expression profiling was performed using Affymetrix Human Genome U133 Plus 2.0 oligonucleotide arrays and standard procedures (Affymetrix, Santa Clara, CA). Array-based gene expression profiling of cerebellar cortex of unaffected individuals was performed using Affymetrix Human U95A oligonucleotide arrays as previously reported . Primary image analysis was performed with Affymetrix GENECHIP 3.2. Images were scaled to an average hybridization intensity of 200, and the threshold for expression was an average hybridization intensity of ≥50.
Array hybridization signals and read frequencies were log10 transformed prior to evaluation of inter-sample differences. Overlaid kernel density estimates, univariate distribution results, Mahalanobis distances, correlation coefficients of pair wise sample comparisons, unsupervised principal component analysis (by Pearson product-moment correlation) and Ward hierarchical clustering of Pearson product-moment correlations of read frequencies were performed with JMP Genomics, Version 3.2 (SAS Institute, Cary, NC). Decomposition of principal components of variance of array signals and read frequencies, FDR corrected analysis of variance (with inclusion of the non-diagnosis components of variance discussed above as fixed effects) and chi squared comparisons of GO annotations were also performed with JMP Genomics, Version 3.2. Patient, sample, and experimental parameters examined to quantify sources of variability in measurements were age, alignment percentage, brain pH, cause of death, cluster station, diagnosis, library creator, number of reads, post-mortem interval, psychotropic medication, race, read length, read quality score, RNA isolation date, RNA integrity number, sample storage duration, sequencing instrument and year sequenced.
A Deo lumen, ab amicis auxilium. We thank Steven Day, Dan Weems, Melodie J. Rice, Terri L. Gomez, Linda Julien, Kamal Gajendran for assistance with Alpheus software development. We thank Becky Archuleta, Dawn Doyle, Pauline Brito, and Leisha Pak for assistance with variant annotation. We thank Eric Vermaas, Victor Quijano, Juying Yan, Leonda Clendenen and Jimmy Woodward for assistance in generating some of the sequencing data. We thank the NIH Neuroscience DNA Microarray Consortium and Translational Genomics Research Institute (Phoenix, AZ) for help with gene expression profiling of human cerebellar tissue. Irina Khrebtukova, Shujun Luo, Lu Zhang and Gary P. Schroth are employees of Illumina Inc. Study data is available at NCBI Gene Expression Omnibus with accession number GSE12297 and at the Short Read Archive with accession number SRA002355.1 and is freely available for query at http://schizophrenia.ncgr.org.
Competing Interests: Irina Khrebtukova, Shujun Luo, Lu Zhang and Gary P. Schroth are employees of Illumina Inc. Wendy Czika, Stanton Martin and Russell D. Wolfinger are employees of SAS Institute, Cary, NC.
Funding: This work was supported by the National Science Foundation, State of New Mexico, National Institutes of Mental Health, National Center for Research Resources. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.