Transcriptome Sequencing of the Anterior Cingulate in Bipolar Disorder: Dysregulation of G Protein-Coupled Receptors

Bipolar disorder is an episodic and debilitating mood disorder that affects approximately 1% of the general population (1, 2). Extensive work has been done to understand the role of genes and regulation in the bipolar disorder brain. Candidate gene and genome-wide microarray expression analyses of postmortem human brains have shown that transcriptional dysregulation plays a role in the etiology of bipolar disorder (for a review, see reference 3). Some notable genes have been identified through these studies (4, 5); however, very few passed corrections for multiple testing, and findings from these studies have largely not been replicated (3, 5, 6). RNA sequencing, a technique that takes advantage of the recent development of high-throughput sequencing technologies, offers a number of advantages in comparison to previous methodologies, such as microarray studies, in that it provides direct estimates of transcript abundance, as well as nucleotide-level sequence. Differential expression can be measured both at the gene and transcript levels, thus providing unbiased and unparalleled evidence for novel RNAs (7) and regulatory mechanisms such as alternative splicing (8), which are not well represented on microarrays. Transcriptome analysis using RNA sequencing in bipolar disorder is timely and important, not only given the power of the technology but also given the need for greater understanding of the bipolar disorder pathophysiology and development of effective treatment options.

In order to understand the role that coding and noncoding RNAs play in brain regulation and how their potential dysregulation could affect brain function and ultimately onset of bipolar disorder, we investigated gene expression changes in postmortem brain tissue from bipolar disorder cases using RNA sequencing. While the precise bipolar disorder neuroanatomical circuits are not exactly known and there are data supporting the involvement of diverse brain regions, there is strong support for the role of the anterior cingulate cortex in the regulation of ideo-affective and mood functions and thus in the neurobiology of bipolar disorder (9, 10). Consequently, we focused this postmortem expression study on this region and found a global pattern of downregulation. Furthermore, we identified several differentially expressed genes and followed up these findings with an in vitro study that showed the mood stabilizers lithium, carbamazepine, and valproate to modulate the expression of these transcripts. By using RNA sequencing we hope to have achieved a more comprehensive level of understanding of the bipolar disorder brain and shed important light on the dysregulated mechanisms, as well as the potential implications for treatment.

Method

Postmortem Brain Samples and High-Throughput Transcriptome Sequencing

Post-mortem brain tissue was obtained from the Douglas-Bell Canada Brain Bank (www.douglasbrainbank.ca) (see the data supplement accompanying the online version of this article). Cases in this study were individuals who had a diagnosis of bipolar disorder type I or type II (N=13). Comparison subjects (N=13) had neither current nor past psychiatric diagnoses. Case and comparison subjects were matched for refrigeration delay, age, and brain pH (see Table S1 in the data supplement). RNA extraction and sequencing library preparation is described in detail in the methods section of the online data supplement. All sequencing was completed on the Illumina HiSeq2000 platform using 100-bp paired-end reads. Reads were aligned to the human genome reference (hg19) using TopHat v2.0.8b (11) (see Table S2 in the data supplement). For gene-level quantification, we employed HTSeq-count version 0.5.4p1 (12) (see Figure S1 in the data supplement). As validation, we also ran Cufflinks v2.1.1 (7) for gene-level counts, as well as for isoform-level counts. For differential expression analysis, fragment counts were normalized across libraries by using the weighted trimmed mean of log expression ratios from the edgeR v3.0.8 R package (13). Furthermore, counts were corrected for heteroscedasticity by employing voom from the limma v3.14.4 R package (14). The linear model used to fit the data included diagnosis, postmortem interval, and RNA integrity number as covariates.

Gene Ontology Analyses

Gene set enrichment analysis was performed using ermineJ v3.0.2 (http://erminej.chibi.ubc.ca/) (15) with a maximum gene set size of 300 and a minimum gene set size of 5, using the best scoring replicate. The precision-recall analysis was run for 10,000 iterations on all the transcripts from the differential expression analysis of the HTSeq genes.

Brain Region Expression Enrichment Analysis

We used the HBAset tool (http://www.chibi.ubc.ca/∼lfrench/HBAset/) (French and Pavlidis, in preparation). HBAset assembles the Allen Human Brain Atlas (16) expression data for each gene in an input set and computes an average expression level for each region. It also computes a probability reflecting the degree of enrichment of expression compared with random background genes.

Neural Progenitor Cell Line Chronic Drug Treatment

Human neural progenitor cells, derived from induced pluripotent stem cell line GM08330 obtained from a healthy male and previously characterized (17), were maintained in standard conditions (see the online data supplement). Chronic (1 week) treatments were performed with drugs commonly prescribed for bipolar disorder: lithium (1 mM), valproate (1 mM), and carbamazepine (50 μM), or no-drug control, after which cell pellets were collected and RNA was extracted. To validate the brain-like properties of neural progenitor cells, we performed immunohistochemistry with the neuron-specific and astrocyte-specific markers MAP2 and GFAP, respectively (see the data supplement).

Quantitative Real-Time Polymerase Chain Reaction (PCR)

Brain RNA for RNA sequencing and quantitative real-time PCR validation was used from the same original extraction. Complementary DNA synthesis and quantitative real-time PCR were performed as previously described (18) (see the data supplement). We investigated the stability of common endogenous genes in each sample set and determined the most suitable to be POLR2A (polymerase [RNA] II [DNA directed] polypeptide A) or ACTB (beta actin) using the NormFinder algorithm (19) (see Table S3 in the data supplement). All graphical data are presented as mean (standard error of the mean). Statistical differences between groups were analyzed by Student’s t tests, Mann-Whitney tests, one-way analysis of variance with Dunnet’s post hoc corrections, or Pearson’s correlation coefficients. Statistical significance was calculated using GraphPad Prism5 and SPSS 20. A p value ≤0.05 was considered statistically significant, and a value ≤0.1 was considered suggestive of a result that nearly reached statistical significance.

Results

Transcriptome Sequencing in the Anterior Cingulate Cortex of Bipolar Disorder Postmortem Brains

We used a directional library protocol that allows distinction of genes overlapping at the same chromosomal locus and a ribosomal-depletion transcriptome selection to allow for identification of non-polyA-tailed RNAs. On average, 391M 100-bp reads were mapped per individual, of which 276M reads had a mapping quality ≥50 based on TopHat (7). For gene-level quantification, we used HTSeq-count (12). We found on average 69.8M paired-end fragments mapping to 60,905 genes, while on average 67.5M fragments did not map to any genes from the reference annotation (see Table S2 in the online data supplement). This is expected, as many reads map to introns undergoing splicing or as-yet uncharacterized transcribed regions (20). We removed RNA transcripts with zero or aberrantly high counts (e.g., RN7SL2, RN7SK) and were left with a total of 27,706 genes. About 61% of the fragments were attributed to protein-coding genes, while the remaining fragments were attributed to other RNA classes, including long intergenic noncoding RNAs (lincRNAs), pseudogenes, antisense RNAs, etc. (see Figure S2 in the online data supplement).

Gene-Level Differential Expression: A Global Pattern for Downregulation

After quantifying gene-level expression using HTSeq-count (12), we identified all differentially expressed transcripts (see Table S4 in the data supplement) and by-and-large found a strikingly prominent global downregulation, with 70% of overall differentially expressed transcripts being downregulated (see Figure S3A in the data supplement) and a comparable enrichment of downregulated genes among the top 100 genes ranked by p value (72 of 100) (see Figure S3B in the data supplement). Of these, 10 were significantly differentially expressed at a false discovery rate ≤0.05 (see Table S4 in the data supplement). All 10 transcripts that passed the stringent false discovery rate cutoff were downregulated, and all were protein-coding (Table 1). In order to replicate our findings obtained with HTSeq, we ran Cufflinks (see Table S5 in the data supplement). Its advantage is that in addition to gene-level, it allows for isoform-level false discovery rate analysis (7). All 10 genes from the initial analysis were found to be differentially downregulated at a false discovery rate ≤0.05 (Table 1). Following validation by quantitative real-time PCR, we found that all 10 genes had fold changes in the expected direction; eight of the 10 gene expression values were correlated with the RNA sequencing data at a suggestive (p≤0.1) level, and eight were nominally significant (p≤0.05) with quantitative real-time PCR (Table 1, [also see Figure S4 in the data supplement]). These results support the accuracy of our expression quantification and the cutoffs we applied in representing the true transcriptomic landscape, as well as providing strong evidence for the importance of the identified genes in the etiology of bipolar disorder. While these genes are largely unstudied, it is worth noting that three (RXFP1, SSTR2, CHRM2) of the top genes belong to class A of the G protein-coupled receptor family of genes, which potentially suggests similar functions.

Replication of Results in External Bipolar Disorder RNA Sequencing Data Sets of the Anterior Cingulate Cortex and Prefrontal Cortex

We sought to replicate our findings in an independent cohort and obtained publicly available RNA sequencing data from the Stanley Neuropathology Consortium Integrative Database Collection (21), consisting of 61 samples (bipolar disorder group: N=26, comparison group: N=35) from the anterior cingulate cortex (22). In spite of the fact that this study had some methodological differences (see the methods section in the data supplement), querying the expression of our top 10 differentially expressed genes in this data set resulted in external validation of our results, with all genes showing consistent direction with our data (downregulation). Nine of the 10 genes were significant at a nominal level (p≤0.05) and one at a suggestive level (p≤0.1) (see Table S6 in the data supplement).

We also compared our findings with an RNA sequencing study of the dorsolateral prefrontal cortex (23), since different cortical regions have similar expression patterns (24–26). We performed a reanalysis of the Akula et al. (23) data (bipolar disorder group: N=11, comparison group: N=8) and found the significantly downregulated transcripts in our data set (1,761 at p<0.01) to be enriched at the top of the downregulated transcripts in the Akula et al. study (area under the curve=0.664) (see Figure S5A in the data supplement). Likewise, Akula et al.’s reanalyzed list of significantly downregulated transcripts (956 at p<0.01) were enriched in the top of the downregulated transcripts in our study (area under the curve=0.611) (see Figure S5B in the data supplement). Conversely, upregulated transcripts showed no enrichment between the two data sets (area under the curve=0.532 using our list of 324 upregulated transcripts and area under the curve=0.533 using the list of 529 upregulated transcripts in the Akula et al. study).

Enrichment of G Protein-Coupled Receptor Pathways in the Differentially Expressed Genes

In order to understand broader patterns of differentially expressed genes in the anterior cingulate cortex, we performed a gene set enrichment analysis on the HTSeq differentially expressed genes. Interestingly, the top two biological processes identified that also passed corrections for multiple testing were “G-protein coupled receptor signaling pathway, coupled to cyclic nucleotide second messenger” (GO: 0007187, corrected p=5.24E-09) and “adenylate cyclase-modulating G protein-coupled receptors signaling pathway” (GO: 0007188, corrected p=2.62E-09) (see Table S7 in the data supplement). We also investigated the top-10 differentially expressed genes to see whether they were enriched in expression in particular brain regions. Using an “expression enrichment” tool for the Allen Human Brain Atlas data, HBAset, we found that these gene expressions were significantly enriched compared with random genes in many cortical regions, including the anterior cingulate cortex (see Figure S6 in the data supplement), one of the top significant regions (p=0.0013; see Table S8 in the data supplement). This suggests that concerted dysregulation of these genes might have effects on multiple cortical regions, including the anterior cingulate cortex.

The Effect of Psychiatric Drugs on Identified Genes

In order to explore the role of medication commonly used to treat bipolar disorder on the dysregulated genes, we investigated the effect of lithium, valproate, and carbamazepine on the expression of the top differentially expressed genes. We performed an in vitro chronic treatment in neural progenitor cell lines that express both the neuron-specific marker MAP2 (microtubule-associated protein 2) and astrocyte-specific marker GFAP (glial fibrillary acidic protein) (Figure 1). All three drugs significantly decreased expression of CHRM2 and VWC2L, while increasing expression of DIRAS2. Only valproate had an effect on SLC7A14 and SSTR2, where both genes were upregulated (Figure 1).

Isoform-Level Analysis

We used Cufflinks to determine isoform-level differentially expressed genes, and after removing transcripts with very high or low expression, 120,845 remained. No transcripts were differentially expressed at a false discovery rate ≤0.05 (see Table S9 in the data supplement; transcripts with a p value ≤0.05). However, many of the top-ranked isoforms by p value belonged to the same genes identified by gene-level analysis. For proof-of-principle, we wanted to see whether we could detect any of these differences with quantitative real-time PCR. Since the gene B3GALT2 only has one isoform, we opted for the second-ranked transcript, ENST00000445907, isoform 1 of the gene CHRM2 (Figure 2A). We designed assays that would target all isoforms (Figure 2B), only ENST00000445907 (Figure 2C), and other isoforms besides ENST00000445907 (Figure 2D). We found that the total gene downregulation (p=0.0004) is driven by the ENST00000445907 isoform (p=0.0002) and not maintained when this isoform is not targeted by the quantitative real-time PCR assay (p=0.1065).

Analysis of Noncoding RNAs

Our approach of transcriptome selection by ribosomal depletion allowed us to also retain non-polyA-tailed RNAs. These include potentially interesting noncoding RNA classes, such as lincRNAs, antisense RNAs, and small nuclear and small nucleolar RNAs, etc. The relative abundance of these classes is lower than coding RNAs (see Figure S2 in the data supplement), though higher than previously suspected (20). A large number of reads was attributed to pseudogenes, a finding that has been previously reported (27). In the differentially expressed analysis, no noncoding RNA passed (false discovery rate ≤0.05) corrections for multiple testing; however, the top two differentially expressed noncoding RNAs, linc-KARS-3 (also known as TCONS_0024733) and linc-SFSWAP-3 (also known as TCONS_0021259), were ranked 14 and 16, respectively, by p value (see Table S4 in the data supplement). Interestingly, RP11–638F5.1 (also known as TCONS_0020164), ranked 23, mapped to the same genomic location on chromosome 12 as linc-SFSWAP-3 and appears to be a shorter isoform of the same locus, sharing two exons (see Figure S7B in the data supplement). Notably, both loci are highly conserved in chimp, gorilla, and rhesus monkeys (see Figure S7 in the data supplement), highlighting their potential functional importance. We found all three of these lincRNAs to be downregulated in the RNA sequencing data and succeeded in showing the same effect by quantitative real-time PCR analysis (linc-KARS-3, p=0.0487) (Figure 3). For the chromosome 12 locus, an assay targeting the two shared exons resulted in a significant downregulation (p=0.0503), while an assay querying just RP11–638F5.1 also resulted in a significant downregulation (p=0.0138) (Figure 3). For technical reasons, an assay specific to linc-SFSWAP-3 could not be designed (see Figure S7B in the data supplement). All quantitative real-time PCR results were significantly correlated with the RNA sequencing data. These results, along with the fact that there was no significant correlation between the two quantitative real-time PCR data sets, suggest that the RNA sequencing findings at this locus are independent of each other.

Discussion

In this study, we investigated the transcriptome of individuals with bipolar disorder in postmortem brains from the anterior cingulate cortex (Brodmann’s area 24) using RNA sequencing. Extensive evidence from microarray and candidate gene studies has demonstrated the role of transcriptional dysregulation in the etiology of bipolar disorder (3). Furthermore, accumulating evidence is starting to point to the key involvement of as-yet uncharacterized noncoding RNAs, in addition to the more commonly studied protein-coding RNAs in psychiatric disorders.

We showed excellent validation of our methods, both at the bioinformatics level across two different pipelines and at the molecular level with quantitative real-time PCR. Furthermore, we replicated our top findings in an external cohort consisting of anterior cingulate cortex samples from 61 postmortem brains from the Stanley Neuropathology Consortium Integrative Database collection, for which transcriptome sequencing findings were recently described in a study conducted by Zhao et al. (22). By-and-large, we found a strikingly prominent global downregulation, with all differentially expressed transcripts that passed multiple testing corrections (false discovery rate ≤0.05) being downregulated, as well as an enrichment of downregulated genes among the top 100 genes ranked by p value and a much more consistent expression pattern across subjects in downregulated genes compared with upregulated ones. The top gene identified by both gene-level bioinformatics pipelines and the isoform-level analysis was B3GALT2 (UDP-Gal:betaGlcNAc beta 1,3-galactosyltransferase, polypeptide 2) (28). Quantitative real-time PCR validation in our cohort was supportive of this finding, and in the Stanley Neuropathology Consortium Integrative Database replication cohort this gene was also downregulated, though at a suggestive significance level. This gene is a member of the beta−1,3-galactosyltransferase (beta3GalT) family, which encodes type-II membrane-bound glycoproteins. Though very little is known about B3GALT2 and associations with bipolar disorder have not yet been documented, genes from this family have been shown to be primarily brain-expressed in the mouse (29). Furthermore, B3galt2 knock-out mice have been shown to have decreased anxiety-related response and increased hyperactivity through the Open Field test, as well as decreased coping response through the Tail Suspension test, indicating a possible depression-like phenotype (30). It is also worth noting that Akula et al. (23), who conducted a bipolar disorder RNA sequencing study in the dorsolateral prefrontal cortex, also detected this gene to be downregulated.

One of our most interesting findings was the differential expression of G protein-coupled receptors, as suggested by the gene-set enrichment analysis that indicated an enrichment of G protein-coupled receptor pathways among the differentially expressed genes, supported by the fact that three of the top differentially expressed genes, SSTR2 (somatostatin receptor 2), CHRM2 (cholinergic receptor, muscarinic 2), and RXFP1 (relaxin/insulin-like family peptide receptor 1), are class A G protein-coupled receptors. These noted drug targets in many disorders, including those afflicting the brain (31), have previously been linked to mood disorders, including bipolar disorder (32). Of the top three G protein-coupled receptors that we have identified in this study, only CHRM2, a muscarinic receptor defined by the binding of acetylcholine and involved in adenylate cyclase inhibition, phosphoinositide degeneration, and potassium channel mediation, has been previously linked to bipolar disorder (33, 34). SSTR2 mRNA levels have been shown to be decreased in the prefrontal cortices of schizophrenia patients (35) and to decrease in response to stress in animal models (36). Interestingly, one recent study by Tomita et al. (32) that focused specifically on microarray profiling G protein-coupled receptor expression in psychiatric postmortem brains found SST (somatostatin) to be one of the top dysregulated genes in the anterior cingulate cortex of bipolar disorder patients. The specific G protein-coupled receptors that emerged in the top of our study were not mentioned in this study, though there was a great deal of overlap in the affected pathways, including inositol polyphosphate phosphatases, neuropeptide Y genes, and proteinphosphatase 1 genes (32). Our study, along with others, suggests that G protein-coupled receptors as a whole should be further investigated in bipolar disorder.

In order to further explore the role of the identified dysregulated genes, we investigated the effect of mood stabilizers lithium, valproate, and carbamazepine on the expression of the top differentially expressed genes through an in vitro chronic treatment experiment in cultured neural progenitor cells. These drugs were selected based on a long-standing history of documented efficacy in the clinical treatment of bipolar disorder (1, 2). Since the majority of differentially expressed transcripts were significantly downregulated in the bipolar disorder brain, we were interested in the possibility that the expression of these genes would be upregulated by mood stabilizers. DIRAS2 was upregulated by all the three drugs tested, while SSTR2 and SLC7A14 were significantly upregulated by valproate only. Very little is known about the implication of DIRAS2 (DIRAS family, GTP-binding RAS-like 2) in the brain, as it has only been studied in adult attention deficit hyperactivity disorder (37). Neither SLC7A14 (solute carrier family 7, member 14) nor SSTR2 has previously been investigated in bipolar disorder or valproate treatment, but the latter is a G protein-coupled receptor belonging to class 4A that has been implicated in adaptive response to stress (36, 38). Furthermore, other somatostatins have been linked to bipolar disorder genetics (39–41). Expression of CHRM2 and VWC2L was decreased by all three drugs, which is unexpected considering that they were also downregulated in the bipolar disorder brains. There are different possible explanations for these findings. Among these, it is possible that some of the brain expression differences observed were compensatory, thus in a similar direction as those observed with action of therapeutic drugs. Similarly, it is possible that substance and past treatment exposure of the bipolar disorder cases may explain some of our significant findings, including CHRM2 and VWC2L. While encouraging, our results provide but the first steps in the attempt to elucidate how commonly prescribed mood stabilizers may influence the expression of genes dysregulated in the bipolar disorder brain.

Finally, to our knowledge, this is the first study in psychiatry that used ribosomal depletion in the preparation of RNA sequencing libraries and thus can quantify all classes of RNA of 100 base pairs or longer, regardless of their possession of a poly-A tail. Furthermore, since we used a very high coverage, we were able to detect some very lowly expressed transcripts, including lincRNAs, antisense RNAs, small nuclear RNAs, and other noncoding classes that are not as abundant as protein-coding transcripts. This analysis is of interest because more and more reports have emerged in the last few years documenting the importance of noncoding RNAs in normal brain development, maintenance, and aging (42–44), as well as a variety of conditions including neurodevelopmental disorders such as autism (45, 46). Though they did not pass multiple testing corrections, the top three (ranked by p value) most significant lincRNAs were significantly downregulated when validated using quantitative real-time PCR, suggesting that these lincRNAs should be further investigated in bipolar disorder. Unfortunately, we have no information from the literature to help us understand how dysregulation of these lincRNAs could be connected to the expression of the coding genes identified. Since the exploration of noncoding RNA species is still in its infancy, undoubtedly computational tools will improve along with our understanding of these RNA classes, allowing us to extract even more valuable knowledge. Further work is warranted to fully exploit the abundance of information collected with total transcriptome sequencing analysis.