Our study provides the first high-resolution transcriptome profiling of the STG from subjects with schizophrenia. An average of 28 million random cDNA sequence reads for each of the 18 samples was generated from 76 bp reads. This delivered sufficient sequence coverage (median 21.0 fold) for transcriptome profiling 
, and provided an unparalleled assessment of the transcriptional complexity in the human STG. For each sample, more than 90% of the total mapped reads uniquely mapped to the reference genome allowing the identification of a comprehensive panel of DEGs and the incidence of alternative promoter usage and splicing variation associated with schizophrenia.
This analysis revealed that the expression of 772 genes was significantly different between cases and controls. These DEGs covered half of the previously highlighted genes with prominent and concordant expression change in schizophrenia, detected by at least two microarray studies across various brain regions 
. This level of consistency between genome wide expression studies was unexpected as previous comparisons between microarray studies typically show more heterogeneity. This concordance between mRNA sequencing and microarray technology is encouraging and suggests there is validity in this approach. Our sequencing data was also supported by qPCR analysis of three genes including DUSP1, BTG2, and EGR4. Interestingly, in rats both cortical BTG2 and DUSP1 expression has been shown to be sensitive to phencyclidine, a noncompetitive antagonist of NMDA receptors 
. Further comparison between mRNA sequencing analysis and existing microarray data both across the various brain regions and within the temporal cortex revealed three groups of genes consistently changed in schizophrenia including 14-3-3 family gene, oligodendrocyte- and myelin-related genes, and genes involved in signalling and synaptic function. Furthermore, for the remaining seven groups of genes ranging from GABA function to ubiquitin and stress genes derived from microarray studies (), at least one gene was found to be overlapped in our DEG list. This overlapping coverage in each functional group was broadly consistent with previous characterization of schizophrenia as a disease of synaptic, oligodendroglial, and metabolic origin. In addition to the comparisons at the discrete gene level, GO enrichment analysis at the gene set level was further carried out to explore the broad areas of schizophrenia dysfunction. This enrichment analysis revealed 55 and 113 over-represented GO terms in schizophrenia at a significance level of FDR<0.05 and FDR<0.4, respectively. These GO terms enriched in schizophrenia were organised into a network consisting of 6 functional clusters including microtubule related motility, ribosome and translation activity, mitochondrial function and ATP production, neurotransmission related functions, synaptic vesicle trafficking, and neural development.
Most striking was the simultaneous detection of 12 genes encoding SNARE complex associated proteins derived from the cluster of synaptic vesicle trafficking, which suggested a prominent role of SNARE complex in presynaptic dysfunction in schizophrenia. Two core SNARE proteins, the synaptosomal-associated protein of 25 kDa (SNAP25) and the vesicle-associated membrane protein (VAMP1, also known as synaptobrevin), were both significantly down-regulated in schizophrenia. The interaction of SNAP25 and syntaxin-1A with VAMP1 leads to the formation of the core SNARE complex between the vesicle and presynaptic terminal membranes 
. Furthermore, 10 SNARE-interacting proteins involved in regulation of vesicle release and recycling were differentially expressed. These genes included: (1) 2 complexins (CPLX1 and CPLX2) and 3 synaptotagmins (SYT1, SYT2, and SYT3), which together regulate calcium-dependent exocytosis 
; (2) one of the synapsins (SYN2) which are involved in vesicle tethering and controlling the number of vesicles available for release 
; (3) synaptophysin (SYP) which binds synaptobrevin (VAMP1) and plays a role in the control of exocytosis 
; (4) one of the RIM family members (RIMS3) which interacts with voltage-dependent Ca2+
channels and regulates neurotransmitter vesicle anchoring 
; and (5) NSF and NAPA(SNAP-α) which bind and then dissociate SNARE complexes for synaptic vesicle recycling 
. Though the impact of antipsychotic treatment on this altered gene expression is unclear, recent studies suggested that at least some of the changes (down-regulation of SNAP25, VAMP1, and SYT2; up-regulation of SYP and CPLX1) are unlikely to arise from the treatment, as animal models showed that haloperidol and/or clozapine increased SNAP25, VAMP, and SYT2 levels, while decreased SYP and CPLX1 levels 
Although altered expression of transcripts encoding presynaptic secretory machinery including SYN2, NSF, synaptotagmin, synaptobrevin, synaptophysin, SNAP25, and RIMS have been reported in various brain regions such as PFC, STG, and amygdala 
, our study is first to link this mechanism to a comprehensive panel of genes focused on the function of SNARE complex in the temporal cortex in schizophrenia. Consistent with our finding, a recent genomic convergence analysis of cerebellar cortices in schizophrenia using mRNA sequencing has identified altered expression of a group of genes involved in presynaptic vesicular transport 
. While this study emphasized genes involved in transport between the trans-Golgi network and synaptic vesicles, our data highlighted the dysregulation of the SNARE complex and associated proteins directly involved in synaptic vesicle release. Comparison between these two mRNA sequencing studies in schizophrenia suggests that different brain regions may share common pathways relating to synaptic abnormalities but involve different pathway members. The deficits in presynaptic gene function could result in compensatory enhancement of neurotransmission by down-regulating RGS4 expression, observed here and in previous microarray studies in PFC and STG 
. The detection of simultaneous alterations of presynaptic genes and RGS4 is supportive of the synaptic-neurodevelopmental model of schizophrenia, which suggests that these influences may affect postnatal synapse formation and pruning, which eventually leads to the disorder 
Further support for the synaptic dysfunction hypothesis of schizophrenia is suggested by the alteration of genes related to GABAergic and glutamatergic function, including GABAergic gene GAD1 (also known as GAD67; ), GABAA
receptor subunits alpha 5 and delta (GABRA5 and GABRD), and ionotropic NMDA type glutamate receptor subunit NR3B (GRIN3B). The reduction of GAD1 observed in our data is one of the most consistent changes in the schizophrenia neocortex 
. Interestingly, while antipsychotic treatment with haloperidol also impact on GABA receptor subunits and GAD1 expression in animal studies 
, these changes were in the opposite direction to our observations in schizophrenia.
The co-detection of eight genes involved in neuronal myelination was also consistent with previous analysis of cortical gene expression in schizophrenia, which reported dysregulation of MAG, MBP, MOBP, PLP1, PLLP, and CNP; transcription factor OLIG1; and oligodendroctye development and myelination associated gene CLDN11 
. This group of 8 myelin related genes, of which 6 significantly down-regulated, provided further evidence for oligodendrocyte-mediated myelination dysfunction.
Other interesting changes included enriched functional clusters related to mitochondria function, microtubule, ribosome and translation and altered expression levels of genes encoding creatine kinase (CKB) and 14-3-3 family members. It is likely that many of these genes are functionally intertwined. For example, our observation of global down-regulation of mitochondrial genes and genes related to energy metabolism in schizophrenia accords with previous studies 
, and suggests that altered energy metabolism may be an integral component underlying the cortical pathophysiology of schizophrenia. The decreased output of ATP may lead to further consequences such as a reduction in CKB involved in the storage of high-energy phosphates as observed in our dataset and bipolar disorder with mitochondrial dysfunction 
. Alternatively, decreased energy production will lead to adaptational changes such as the up-regulation of 14-3-3 family members (YWHAH and YWHAE), which are involved in mitochondrial function 
Coincidently, 9% of DEGs have also been indicated as putative schizophrenia susceptibility genes (Table S1
), including RGS4, NRGN, and APOE which are all ranked in the top 45 list of SZGene (http://www.szgene.org
). Furthermore, 5 out of 12 SNARE associated genes (SNAP25, CPLX1, CPLX2, SYN2, and SYP) and 7 myelin related genes, were all implicated in genetic association/linkage studies 
. These observations suggested that allelic variation might influence the expression of these genes as putative expression-quantitative trail loci, however, further investigation of the association between the differences in allele frequency and the gene expression would be required to establish a functional link.
Significant differences in alternative splicing and promoter usage between cases and controls were identified. Among these, DCLK1 is of particular interest as it is expressed in the CNS and plays important roles in brain development 
. DCLK1 contains a doublecortin-like (DCX) domain at the N-terminus and a kinase domain at C-terminus. The DCX domain at N-terminus was shown to be involved in the regulation of microtubule polymerization and neuronal migration 
. One of the C-terminal isoforms has also been suggested to be a candidate neural plasticity gene with a potential role in synaptic remodelling 
. DCLK1 has a few transcript isoforms that include either its full-length, or only the N-terminal DCX-like region or the kinase-encoding C terminus alone 
. Our data suggests some bias for TSSIII usage in schizophrenia, producing both D02 (full length DCLK1) and D04 (N-terminal DCL-like region) variants, which have been shown to prevent apoptosis in neuroblastoma cells 
. Under TSSIII, schizophrenia subjects had a higher percentage of D02 variants (case 34.4%) than the controls (6.3%). Interestingly, DCLK1 expression was shown to be increased in the mouse brain by clozapine and olanzapine, suggesting that reduction observed in the STG was not likely to be a result of antipsychotic medication 
In addition, we also observed the differential splicing of PLP1 between cases and controls. This gene is noteworthy as it encodes one of the major components of myelin protein, the proteolipid protein (PLP) and has been reported to be down-regulated in a number of previous studies of schizophrenia 
. These two isoforms of PLP1 arise from differential use of two 5′ splice sites in exon 3, with the upstream site used for production of DM20 and the downstream splice site for PLP mRNA 
. Although the functional difference between PLP and DM20 for the development of compact myelin remains to be elucidated, it has been reported that both forms are coordinately expressed during early and late stages of myelination. During embryonic stages of development, DM20 is the predominant product 
, but is overtaken in postnatal stage by PLP, which is highly expressed in oligodendrocytes 
. This observation is consistent with our detection of 78.5% and 83.4% of mRNA isoform encoding PLP in both control and case, respectively. It has also been suggested that DM20 and PLP may be identical in topology, but differ in the conformation due to the extra 35 amino acid cytoplasmic peptide in PLP 
. This peptide with the normal function of detecting conformational changes and triggering intracellular events such as membrane compaction, imparts a susceptibility to mutations in extracellular domains. The functional significance of the two distinct isoforms in relation to schizophrenia warrants further investigation.
A few limitations concerning the use of postmortem brain tissue should be noted here. First, even though all the variables thought to influence RNA integrity and gene expression were controlled for as much as possible in this study, we cannot completely rule out the impact of antipsychotic medication on changes in gene expression in post-mortem brain tissue. Furthermore, two of the schizophrenia subjects died as a result of suicide, which may potentially contribute other changes to the transcriptome that are independent of the underlying neuropsychiatric diagnosis. It is also possible that some of the results reflect differences in the ratios of various cell types in the tissue as a consequence of differential shrinkage in specific cell types.
In summary, this is the first study to comprehensively profile mRNA expression in STG in schizophrenia using ultra high-throughput sequencing technology. The predominant change in gene expression pattern accorded well with previous microarray studies in its characterisation of schizophrenia as a disease of synaptic, oligodendroglial, and metabolic origin 
. In contrast to expression by hybridisation, the sequencing approach demonstrated greater sensitivity and specificity. This enabled us to reveal the significance of the SNARE complex in synaptic abnormalities and indicated a schizophrenia candidate gene DCLK1 believed to be involved in neuronal development. It also enabled the detection of novel disease-associated alternative splicing and promoter usage, suggesting that aberrant RNA processing could play a significant role in the complex pathophysiology of schizophrenia. Systematic investigation and functional analysis of alternative splicing and promoter usage in relation to the neuropathology of schizophrenia is needed to advance our current understanding of the disease pathophysiology and will no doubt comprise an essential part to future drug discovery efforts.