Comparative genomic evidence for the involvement of schizophrenia risk genes in antipsychotic effects.

Genome-wide association studies (GWAS) for schizophrenia have identified over 100 loci encoding >500 genes. It is unclear whether any of these genes, other than dopamine receptor D2, are immediately relevant to antipsychotic effects or represent novel antipsychotic targets. We applied an in vivo molecular approach to this question by performing RNA sequencing of brain tissue from mice chronically treated with the antipsychotic haloperidol or vehicle. We observed significant enrichments of haloperidol-regulated genes in schizophrenia GWAS loci and in schizophrenia-associated biological pathways. Our findings provide empirical support for overlap between genetic variation underlying the pathophysiology of schizophrenia and the molecular effects of a prototypical antipsychotic.

Introduction

A major goal of human GWAS is to identify potential therapeutic targets for common, complex diseases. While genomic regions that reach genome-wide significance explain only a small fraction of disease risk, many nonetheless encode proteins that make effective drug targets (1). For example, a common genetic variant in HMGCR has a small (~5%) but significant (P=1×10−30) influence on low density lipoprotein levels but its inhibition by statins effectively treats hyperlipidemia (2). Another example comes from one of the earliest GWAS, which identified a common variant in the complement factor H gene for age-related macular degeneration (3–5) that led to targeting the complement cascade for the treatment of age-related macular degeneration (6).

Schizophrenia is a chronic, severe, and disabling brain disorder that has a median morbid lifetime risk of 0.72% (7, 8). The most recent schizophrenia GWAS identified 108 genome-wide significant loci encoding over 500 genes (9). Because screening these putative schizophrenia risk genes for those that are disease- and therapeutically-relevant is a substantial task, we adopted an alternative, in vivo molecular approach. We examined the overlap between schizophrenia risk genes and their orthologous mouse genes whose striatal expression was significantly altered following chronic haloperidol treatment. We improved upon prior studies of the effects of chronic antipsychotic exposure on gene expression in rodent brain (Supplementary Tables 1–2) by use of a better detection technology (RNA-sequencing), larger sample sizes to enable detection of more subtle effects (total N=38 mice), and examined several tissues including striatum, whole brain, and liver for comparative controls (Supplementary Table 3). We provide evidence that schizophrenia risk genes and pathways are relevant to haloperidol effects.

Materials and Methods

The goal of this study was to evaluate the effects of chronic administration of the antipsychotic haloperidol versus vehicle in mice. We have shown that we can reliably administer human-like steady-state concentrations of haloperidol (10–12). All experimental procedures were randomized to minimize batch artifacts (13) (e.g., assignment to haloperidol or vehicle, cage, order of dissection, RNA extraction, and assay batch). Experimenters were blind to treatment status.

We focused on striatum as it is relatively dense with dopaminergic neurons and a key site of action of the dopamine receptor antagonist haloperidol (14) (confirmation of choice of tissue is described in the Discussion). We also evaluated whole brain (to measure effects outside the striatum) and liver (to identify brain-independent and hepatic alterations consequent to chronic xenobiotic administration). For striatal samples, we use RNA sequencing (RNA-seq) to comprehensively identify differential gene expression resulting from chronic haloperidol exposure. For whole brain and liver, we used gene expression microarrays, which provide an inexpensive transcriptome evaluation (albeit with lesser dynamic range), and correlates well with RNA-seq (mean r=0.87 across 88 mouse brain samples assayed with both methods) (15).

Mice

All animal work was conducted in compliance with national guidelines (Institute of Laboratory Animal Resources, 1996) and approved by the UNC Institutional Animal Care and Use Committee. The study design is summarized in Supplementary Table 3. For striatal RNA-seq we chronically treated mice with haloperidol (N=16) or vehicle (N=12). All striatal samples were assayed using RNA-seq. Independent mice (N=20) were used to collect whole brain and liver (left lobe) from mice per treatment group for expression microarray analysis. To minimize the effects of the estrus cycle and other sources of heterogeneity, we evaluated male C57BL/6J mice (shipped at 6 weeks of age, Jackson Laboratory, Bar Harbor, ME, USA) for both RNA-seq and microarray experiments. Animals were maintained in standard environmental conditions (14-hour light/10-hour dark schedule, temperature 20–24°C, and 40–50% relative humidity). Mice were housed four per cage in standard 20 cm × 30 cm ventilated polysulfone cages with laboratory grade Bed-O-Cob bedding. Water and Purina Prolab RMH3000 were available ad libitum. Mice were housed in groups of four mice per cage (two haloperidol and two vehicle).

Haloperidol exposure

Eight week old mice were implanted with slow-release haloperidol pellets (3.0 mg/kg/day; Innovative Research of America; Sarasota, FL, USA)(16) or vehicle and treated for 30 days for a chronic haloperidol administration paradigm. Pellets were implanted subcutaneously, centrally above the scapulae, under isoflurane anesthesia and the incision sealed with VetBond (3M, St. Paul, MN, USA). We have previously demonstrated that this procedure reliably yields human-like steady-state concentrations of haloperidol in blood plasma and brain tissue, and that this results in vacuous chewing movements (an established model of extrapyramidal symptoms (17)) in C57BL/6J mice (10, 11, 18).

Tissue collection

After 30 days of exposure to haloperidol or vehicle (12 weeks of age), mice were sacrificed by cervical dislocation without anesthesia (to avoid effects on gene expression). All mice were euthanized between 8:00 AM and 12:00 PM, immediately after removal from the home cage. Tissues were dissected within 5 minutes of death, snap-frozen in liquid nitrogen, and pulverized using a BioPulverizer unit (BioSpec Products, Bartlesville, OK). Tissues collected were striatum, whole brain, or liver (left lobe). Striatum and whole brain were collected from separate animals. The striatum dissection consisted of capturing a 2 mm thick coronal section (Bregma coordinates +1.0 to −1.0) followed by manual isolation of the striatal region per a mouse brain atlas (Figure 1a) (19). Left and right striatum were pooled for each animal.

Figure 1

Schizophrenia risk genes and historical candidate genes are differentially expressed following chronic haloperidol exposure. (a) Overlap between schizophrenia GWAS loci from the Psychiatric Genomics Consortium (9) and one-to-one orthologous mouse genes with altered expression following haloperidol (q < 0.1). A significant enrichment was seen with striatum RNA-seq data (P=0.0004), but not with whole brain expression data (P=0.45). (b) Manhattan plot of schizophrenia GWAS results (9) showing 39 differentially regulated genes following chronic haloperidol exposure (red: increased, blue: decreased expression, q < 0.1). Underlined genes are located in a multi-genic GWAS locus. Genes of note are highlighted with a light green arrow. (c) Mouse genes with the most significant change in expression following chronic haloperidol exposure (all genes with q <0.05), ranked by fold change. This list is significantly enriched for historical schizophrenia candidate genes (◆ from SzGene database (43); ❖; from top 25 most studied candidate genes (42)).

RNA-seq

Total RNA was extracted from striatum using the Total RNA Purification 96-Well Kit (Norgen Biotek, Thorold, ON, Canada). RNA concentration was measured by fluorometry (Qubit 2.0 Fluorometer, Life Technologies, Carlsbad, CA) and RNA quality was verified using microfluidics (Bioanalyzer, Agilent Technologies, Santa Clara, CA). Barcoded libraries were created using Illumina (San Diego, CA) TruSeq Stranded mRNA Library Preparation Kit v2 with polyA selection using 1 μg total RNA as input. Equal amounts of all barcoded samples were pooled, to account for lane and machine effects. This pool was sequenced on eight lanes of the Illumina HiSeq 2000 (100 bp single-end reads). See Supplementary Table 4 for alignment summary of post-QC RNA-seq samples.

We mapped lane-level reads to the mouse genome (mm9) using Tophat (20) (v2.0.6, default parameters). Using samtools (21), we removed reads with quality score <10 or potential PCR duplicates, and a median of 87% of the reads mapped uniquely to the genome (range: 79–91%). The eight aligned BAM files for each sample (one per lane) were highly correlated and merged into one BAM file for each sample. Mapped reads were summarized into gene-level expression estimates of total read count (TReC). TReC is the number of reads that overlapped exonic regions of a gene (using the R package isoform (22)). Ensembl gene models were used (release 67, URLs). This yielded summarized read counts for 26,252 genes. We excluded genes with low expression levels (sum of TReC across all samples <50), resulting in 17,209 genes for analyses. TReC data were normalized using the weighted trimmed mean of M-values scale-normalization method in EdgeR (23). We tested for differential gene expression in striatum using the negative binomial generalized linear model approach in EdgeR (23) with gene-wise dispersion applied (23, 24). The drug effect was evaluated using log-likelihood statistics comparing null and alternative models. False discovery rate (FDR) correction was applied to gene-based P-values to account for multiple comparisons (R package qvalue) (25).

Gene expression arrays

Total RNA was extracted from ~25 mg of powdered tissue from whole brain using automated instrumentation (Maxwell 16 Tissue LEV Total RNA Purification Kit, Promega, Madison, WI). RNA concentration was measured and RNA quality verified as described above. Whole brain RNA from 20 male C57BL/6J mice (10 haloperidol, 10 vehicle) was hybridized to Affymetrix Mouse Gene 2.0 ST 96-Array Plate arrays using a GeneTitan instrument according to manufacturer instructions. As in our prior reports (15, 26), we used robust multiarray average method in the Affymetrix gene expression console to estimate normalized expression levels (default settings, median polish and sketch-quantile normalization). We excluded probes containing any known SNPs in C57BL/6J (27) resulting in 24,464 probesets for analysis. We searched for outliers using principal component analysis (PCA) and hierarchical clustering (R function hclust), and identified only one outlier in a liver sample. We evaluated potential confounding variables by examining the relationship between PC1-PC10 from the expression data and each variable. We found that all potential confounders had minimal impact on gene expression in whole brain and liver. For unmeasured confounders, we performed surrogate variable analysis (28). Surrogate variables 1 and 2 explained the majority of variation in the residual from a model including haloperidol treatment. We used the following model to identify genes displaying differential expression: y = β, where “drug” is an indicator of haloperidol exposure and “sv1” and “sv2” are the first two surrogate variables. FDR correction was applied to transcript-based P-values to account for multiple statistical comparisons (R package qvalue) (25).

Differentially expressed genes and GWAS results

We used INRICH (29) and MAGMA (30) to test for enrichment of GWAS signals in transcripts showing differential gene expression in haloperidol versus vehicle. INRICH evaluates whether a given gene set or pathway has an enrichment of smaller GWAS association P-values than expected by chance (accounting for gene size, SNP density, LD, and pathway size). MAGMA combines the GWAS P-values for each gene (10kb upstream to 1.5kb downstream) into a gene-level P-value and accounts for correlations between SNPs based on the LD (using 1000 Genomes Project European reference panel Phase 3 (31)). MAGMA applies a linear regression framework to test whether gene sets are significantly associated with a trait with respect to the rest of the genome or to some subset of genes. The latter was performed, by adding a covariate consisting of all genes located at the GWAS loci. To determine whether the significant enrichment of schizophrenia GWAS was specific to schizophrenia, we performed enrichment testing with INRICH using GWAS results from other psychiatric disorders (autism, bipolar disorder, major depressive disorder, Alzheimer’s disease; URLs)(32–34) and non-psychiatric traits (height, BMI, and type 2 diabetes; URLs)(35, 36). Because high LD and high gene density in the MHC may influence these analyses, we performed enrichment tests with inclusion and exclusion of the MHC region.

Functional enrichment analysis

We used ConsensusPathDB (release MM9 (11.10.2013) version for mouse) (37) to test differentially expressed genes (FDR q < 0.1) for enrichment in Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, WikiPathways, MouseCyc (38), and drug-gene interaction databases. Functional clustering analysis was performed separately on up- and down-regulated genes using a hypergeometric test to examine whether overlap between our list of genes and those present in each reference category was higher than expected by chance. The background gene list included genes assessed in our expression experiment.

Gene co-expression network analysis

To explore higher-order interactions in an unbiased fashion, we used WGCNA (39) (URLs). After removing transcripts with low expression levels, we applied blockwise Module function (power=6). For striatum RNA-seq data, we included genes with a TReC sum > 50 across all samples and then fit a negative binomal general linearized model for each gene using batch and RNA integrity number as covariates. We applied WGCNA to the residuals after subtracting the fitted values from original raw data. For the whole brain microarray data, we used 15,185 probesets that had a mean log2-transformed and normalized expression level > 6.5 across all samples.

Results

We performed RNA-seq of striatal brain tissue from adult male C57BL/6J mice chronically treated (30 days) with implanted haloperidol (N=16) or vehicle pellets (N=12). Additional mice were examined for differential RNA expression in whole brain and liver (N=20 mice, Supplementary Table 3). We first examined two positive control genes (Drd2 and Nts), and found that both genes showed a striatum-specific increase in expression after chronic haloperidol exposure, as expected from prior studies (40, 41) (Supplementary Table 5). We found that the transcriptional effects of chronic haloperidol exposure were brain-specific (Supplementary Figure 1) and most pronounced in striatum (Supplementary Table 6). Haloperidol-regulated genes were enriched for orthologous schizophrenia GWAS risk loci (Figure 1a, INRICH P=0.0004, MAGMA P=0.0003, 32 loci, 39 genes, Supplementary Tables 7–8). This enrichment remains when the MHC region is excluded (INRICH P=0.0006, MAGMA P=0.0001). These effects were not seen in whole brain (INRICH P=0.45) or liver (INRICH P=0.95), suggesting that these effects are anatomically specific. Although haloperidol-regulated and schizophrenia-associated genes tend to be expressed in striatum, enrichment for schizophrenia GWAS loci remained significant when restricting the analysis to genes expressed in mouse striatum (INRICH P=0.0009).

Consistent with our hypothesis that the overlap is largely specific to schizophrenia, only nominal associations were seen with a co-heritable condition, bipolar disorder, and no significant enrichment was observed for GWAS loci from human studies of autism, major depressive disorder, Alzheimer’s disease, height, BMI, or type 2 diabetes mellitus (Supplementary Table 9). Excluding the MHC region had a minimal effect on the significance of enrichment (no MHC: P=0.0006, 31 overlapping loci vs. included; with MHC P=0.0004, 32 overlapping loci). Most genes overlapping schizophrenia GWAS loci were down-regulated (Figure 1b, 27 of 39 genes), and 15 genes were located in a human multi-genic locus associated with schizophrenia. Therefore, these results provide support for plausible target genes within multi-genic loci along with hypotheses of direction of association.

We found that haloperidol-regulated genes are also enriched for historical schizophrenia candidate genes. Figure 1c shows genes with the greatest fold-change after chronic haloperidol exposure (full list in Supplementary Table 10). Three genes (CHRNA7, HTR2A, and SLC6A4) are among the top 25 most studied schizophrenia candidate genes (42) (P=0.00018) and eight (CARTPT, CHRNA7, GABRA6, GSTT2, HTR2A, NTS, PENK, SLC6A4) are among the 864 orthologous schizophrenia candidate genes in the SzGene database(43) (P=0.05, Figure 1c). Notably, genetic evidence for association with schizophrenia for most of these genes is currently lacking (42). Many of these candidate genes were proposed based on pharmacological properties (Supplementary Figure 2).

Discussion

We used an in vivo molecular approach to determine if there is significant overlap between genes and pathways involved in schizophrenia risk and those regulated by chronic antipsychotic treatment. We performed RNA sequencing of brain tissue from mice chronically treated with haloperidol and found that haloperidol-regulated genes are overrepresented within schizophrenia GWAS loci and schizophrenia-associated biological pathways. Our findings indicate overlap between genetic variation underlying the pathophysiology of schizophrenia and the molecular effects of a prototypical antipsychotic.

Our main analyses focused on striatum. We confirmed our choice of this tissue using single-cell RNA-seq in mouse brain (44, 45) where we demonstrated marked enrichment of antipsychotic drug targets (46) and differentially-expressed genes from this experiment: both analyses pointed at the dominant cell type in ventral striatum, medium spiny neurons (particularly those expressing Drd2) (47).

If such convergence occurs at the levels of biological pathways, it follows that some genes might be important for antipsychotic effects but do not harbor common variants that increase risk for schizophrenia. For example, the serotonin transporter and nicotinic acetylcholine receptor subunit alpha-7 are implicated in antipsychotic pharmacology (48, 49), show marked expression changes after chronic treatment with haloperidol but, at present, have no genetic association with schizophrenia (additional examples include HTR2A and NTS, Supplementary Figure 3–4). One might also expect genes within the same pathway to show gene expression changes and association with schizophrenia. Such is the case for the dopamine receptor D2 (DRD2), the direct target of all effective antipsychotics, and for synaptosomal-associated protein, 91kDa a novel synaptic vesicle protein (SNAP91; Supplementary Figure 3). Additional examples include CACNA1C and GRIN2A (Supplementary Figure 5).

Additional support for a role of putative schizophrenia risk genes in antipsychotic action is derived from gene co-expression network analyses (Supplementary Figures 6–8) and functional enrichment analyses (Supplementary Tables 11–12, Supplementary Figures 9–10). Overlap between haloperidol-regulation and schizophrenia risk extends to biological pathways and perhaps even to the composition of some multi-subunit receptors. For example, there is enrichment of the pathway representing nicotinic acetylcholine receptor regulation of dopaminergic synapses, an active process in the striatum (50) (Supplementary Figure 8a). The nicotinic receptor α4α6β2β3 (which is critical to striatal dopamine release) (51, 52) contains subunits encoded by a gene with differential expression but no genetic association (Chrna6) and a gene with the inverse pattern of findings (Chrna4, Supplementary Figure 8b). Furthermore, the net expression changes we observe in this pathway suggest that chronic haloperidol exposure likely decreases dopamine release in the striatum (50, 53). We also showed that the enrichment of haloperidol-regulate genes is specific to SCZ GWAS (Supplementary Table 9). A recent study used MAGMA to identify new drug targets using the PGC2 SCZ data (54). Interestingly, the authors found that targets of antipsychotics were enriched for association with the SCZ GWAS data. When they looked at druggable targets, they found that the PGC2 SCZ GWAS findings were associated with antipsychotics and anticonvulsants, as well as drugs targeting calcium channels and nicotinic acetylcholine receptors (54).

In summary, by integrating human genetic findings and mouse in vivo expression data, we provide evidence that some schizophrenia risk genes may be involved in chronic effects of haloperidol. Our findings suggest targets for antipsychotic drug development (e.g., genes highlighted in Figure 1b and pathways in Supplementary Figure 8a). Our results support the ongoing development of α7-nicotinic acetylcholine receptor agonists for cognitive enhancement in schizophrenia (55), and suggest that α6 is also a potential target, but this will require further experimental molecular data. We also show that the mouse can be a suitable and efficient model organism in which to use human GWAS results to learn more about drug mechanisms and to support compound development.