A whole genome SNP database and computational genetic mapping were used to analyze the murine genetic model of HIT. Guided by the mouse genetic analysis, we demonstrate that genetic variation within an ABC-drug efflux transporter (Abcb5) affected susceptibility to HIT. In situ hybridization results reveal that Abcb5 is expressed in brain capillaries, and by cerebellar Purkinje cells. We also analyzed chromosome substitution strains, imaged haloperidol abundance in brain tissue sections and directly measured haloperidol (and its metabolite) levels in brain, and characterized Abcb5 knockout mice. Our results demonstrate that Abcb5 is part of the blood-brain barrier; it affects susceptibility to HIT by altering the brain concentration of haloperidol. Moreover, a genetic association study in a haloperidol-treated human cohort indicates that human ABCB5 alleles had a time-dependent effect on susceptibility to individual and combined measures of HIT. Abcb5 alleles are pharmacogenetic factors that affect susceptibility to HIT, but it is likely that additional pharmacogenetic susceptibility factors will be discovered.

We know very little about the genetic factors affecting susceptibility to drug-induced central nervous system (CNS) toxicities, and this has limited our ability to optimally utilize existing drugs or to develop new drugs for CNS disorders. For example, haloperidol is a potent dopamine antagonist that is used to treat psychotic disorders, but 50% of treated patients develop characteristic extrapyramidal symptoms caused by haloperidol-induced toxicity (HIT), which limits its clinical utility. We do not have any information about the genetic factors affecting this drug-induced toxicity. HIT in humans is directly mirrored in a murine genetic model, where inbred mouse strains are differentially susceptible to HIT. Therefore, we genetically analyzed this murine model and performed a translational human genetic association study.

These findings indicate that Abcb5 is a component of the blood-brain barrier in mice and suggest that genetic variants in the gene encoding this protein underlie, at least in part, the differences in susceptibility to haloperidol-induced toxicity seen among inbred mice strains. Moreover, the human genetic association study indicates that a specific ABCB5 allele also affects the susceptibility of people to haloperidol-induced toxicity. The researchers note that other ABCB5 alleles or other genetic factors that affect haloperidol-induced toxicity in people might emerge if larger groups of patients were studied. However, based on their findings, the researchers propose a new model for the genetic mechanisms that underlie inter-individual and cell type-specific differences in susceptibility to haloperidol-induced brain toxicity. If confirmed in future studies, this model might facilitate the development of more effective and less toxic drugs to treat a range of brain disorders.

The researchers used a database of genetic variants called single nucleotide polymorphisms (SNPs) and a computational genetic mapping approach to show first that variations within the gene encoding Abcb5 affected susceptibility to haloperidol-induced toxicity (indicated by changes in the length of time taken by mice to move their paws when placed on an inclined wire-mesh screen) among inbred mouse strains. Abcb5 is an ATP-binding cassette transporter, a type of protein that moves molecules across cell membranes. The researchers next showed that Abcb5 is expressed in brain capillaries, which is the location of the blood-brain barrier. Abcb5 was also expressed in cerebellar Purkinje cells, which help to control motor (intentional) movements. They also measured the measured the effect of haloperidol and the haloperidol concentration in brain tissue sections in mice that were genetically engineered to make no Abcb5 (Abcb5 knockout mice). Finally, the researchers investigated whether specific alleles (alternative versions) of ABCB5 are associated with haloperidol-induced toxicity in people. Among a group of 85 patients treated with haloperidol for a psychotic illness, one specific ABCB5 allele was associated with haloperidol-induced toxicity during the first few days of treatment.

Although drugs have been developed to treat various brain disorders, more active and less toxic drugs are needed to improve the treatment of many if not most of these conditions. Unfortunately, relatively little is known about how the blood-brain barrier regulates the entry of drugs into the brain or about the genetic factors that affect the brain’s susceptibility to drug-induced toxicities. It is not known, for example, why about half of patients given haloperidol—a drug used to treat psychotic disorders (conditions that affect how people think, feel, or behave)—develop tremors and other symptoms caused by alterations in the brain region that controls voluntary movements. Here, to improve our understanding of how drugs enter the brain and impact its function, the researchers investigate the genetic factors that affect haloperidol-induced toxicity by genetically analyzing several inbred mouse strains (every individual in an inbred mouse strain is genetically identical) with different susceptibilities to haloperidol-induced toxicity and by undertaking a human genetic association study (a study that looks for non-chance associations between specific traits and genetic variants).

The brain is the control center of the human body. This complex organ controls thoughts, memory, speech, and movement, it is the seat of intelligence, and it regulates the function of many organs. The brain comprises many different parts, all of which work together but all of which have their own special functions. For example, the forebrain is involved in intellectual activities such as thinking whereas the hindbrain controls the body’s vital functions and movements. Messages are passed between the various regions of the brain and to other parts of the body by specialized cells called neurons, which release and receive signal molecules known as neurotransmitters. Like all the organs in the body, blood vessels supply the brain with the oxygen, water, and nutrients it needs to function. Importantly, however, the brain is protected from infectious agents and other potentially dangerous substances circulating in the blood by the “blood-brain barrier,” a highly selective permeability barrier that is formed by the cells lining the fine blood vessels (capillaries) within the brain.

Funding: GP and MZ were partially supported by funding from a transformative RO1 award (1R01DK090992-01) DLD was supported by a King Abdullah University of Science and Technology (KAUST) research grant under the KAUST Stanford Academic Excellence Alliance program. ST was supported by Stanford’s CURIS program for summer research interns. LSE was supported by a postdoctoral fellowship from the Center for Molecular Analysis and Design (CMAD). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Because no pharmacogenetic factors for HIT have previously been identified in humans or mice, we wanted to analyze this murine genetic model of HIT. Of concern, mouse genetic studies have, overall, produced rather disappointing results; they have definitively identified <1% of the genes affecting the studied traits [ 7 ]. Beyond the well described problems associated with murine linkage studies [ 8 ], several analyses have predicted that the conventional methods used to analyze genome wide association studies (GWAS) using inbred mouse strains will produce even more disappointing results, owing to their low power and a high false positive rate [ 7 , 9 , 10 ]. However, the availability of next generation sequence (NGS) data, which can provide complete genome-wide variant information for an expanded number of inbred strains [ 11 ], could enable murine GWAS to identify causative genetic factors for many biomedical traits [ 12 ]. The detailed genetic variation information obtained using NGS has increased power for pinpointing underlying causal variants relative to studies using SNP array data [ 13 ]. Unlike the conventional murine GWAS methods, haplotype-based computational genetic mapping (HBCGM) [ 14 ] has provided an alternative method for analysis of murine GWAS data that has identified causative genetic factors for many biomedical traits (reviewed in [ 12 ]). HBCGM first organizes the allelic information into haplotype blocks, which facilitates the identification of the causative genetic factors through correlation of the phenotypic data with the haplotype-based pattern of genetic variation. Since our prior studies analyzed much less complex datasets, which had fewer SNPs and haplotype blocks, we did not know if HBCGM could analyze the far more complex datasets produced using NGS data. To determine if HBCGM could be successfully combined with an SNP database generated from whole genome sequence data, we analyzed a murine genetic model of HIT in order to identify causative genetic factors. A translational human genetic association study was then performed to determine if similar genetic factors affected human susceptibility to HIT.

Relatively little is known about the mechanisms regulating the entry of many drugs into the brain, or about the genetic mechanisms affecting susceptibility to many drug-induced central nervous system (CNS) toxicities. Because of this lack of information, candidate drugs affecting the CNS have the lowest rate of survival during clinical development (<2%), and drug distribution across the blood-brain barrier has been identified as a contributing cause for drug candidate failure [ 1 ]. The blood-brain barrier is a very complex multicellular vascular structure, which separates the brain from the peripheral blood by actively controlling the entry and efflux of drugs [ 2 ]. Therefore, the mechanistic insight obtained by identification of genetic factors affecting a drug-induced CNS toxicity could aid in the development of new medicines and could improve our ability to optimally utilize existing treatments. For example, haloperidol is a potent dopamine receptor antagonist that is used to treat psychotic disorders [ 3 ]. However, haloperidol-induced alterations in the extrapyramidal motor system depress the ability to initiate voluntary movements. This effect triggers the characteristic extrapyramidal symptoms (EPS) that develop in 40%–76% of chronically treated human patients, which include tremors, Parkinsonian rigidity, and decreased spontaneous movement [ 4 , 5 ]. The same symptoms appear in some (but not all) of the 27 inbred strains that were treated with haloperidol [ 6 ]. In this murine genetic model, the time required (latency) for a mouse to move all four paws after being placed on an inclined wire-mesh screen was measured as an index of the haloperidol-induced Parkinsonian rigidity that is observed in haloperidol-treated human patients. Although all strains had plasma drug levels that were similar to those measured in human patients, the extent of the haloperidol-induced toxicity (HIT) varied in a highly (85%) heritable manner [ 6 ] among the inbred strains.

Methods

Next-Generation Sequencing of the Genome of Nine Inbred Mouse Strains Genomic DNA was obtained from the Jackson Laboratory (http://www.jax.org). A total of 3 μg of DNA from each of nine inbred mouse strains was sonicated using a Covaris instrument with the following settings: duty cycle, 10%; intensity, 4; cycles per burst, 200; and time, 60 seconds. After end repair of the DNA fragments and addition of “A” base to the 3′ ends, standard Illumina genomic DNA sequencing adaptors were ligated. The samples were run on a 3% agarose gel, gel slices between 300–500 bp were excised, and the DNA was purified with a Qiagen column. The DNA was divided into three aliquots of equal quantity, and amplified separately using one of the following three DNA polymerases: Phusion HF (Finnzymes), Kapa HiFi (Kapa Biosystems), and Herculase (Agilent). For analysis of the SJL strain, Phusion GC (Finnzymes) was also used to amplify a separate aliquot. The PCR conditions for amplification were as follows: 98°C for 3 minutes, followed by nine cycles of 98°C for 20 seconds, 62°C for 30 seconds, and 72°C for 45 seconds, and then a final 72°C for 5 minutes. The amplified products were separated on an agarose gel, slices between 400 and 600 bp were excised, and the DNA was eluted from the gel slices and sequenced on an Illumina HiSeq 2000 instrument.

Sequence Analysis and Variant Calling The quality statistics for the allele calls at each position were evaluated using the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). A rapid decay in sequence quality was observed toward the end of each read, usually within the last 10 (or sometimes 20) bps. To avoid artifacts caused by low quality sequence data, the last 10 or 20 bps of each read were trimmed, based upon a visual inspection of the quality statistics plot for the reads in each lane. The extent of trimming was designed to ensure that the quality scores should be ≥30 in the 1st quartile and ≥20 in the 5th percentile. All sequences were trimmed by 10 bps; with the exceptions being one of the three lanes for the MA/MyJ and MRL, and two of the three SM/J lanes. Since a much larger number of sequences were obtained for SJL, SJL sequences were trimmed by 20 bps to ensure the overall quality of the base calls. The read sequences were then aligned to the reference C57BL/6J mouse genome sequence (MGSCv37 assembly) using BWA [15]. Optical and PCR duplicates were then marked using the Picard tool (http://picard.sourceforge.net/). The sequence variants (SNPs and INDELs) for each strain were identified using Samtools [16], and the default settings were applied. The raw sequence data for the 17 inbred mouse strains described in [11] were obtained from the European Nucleotide Archive (http://www.ebi.ac.uk/ena/), and the same data processing procedures were applied, but these sequences were not trimmed. To ensure high quality variant calls, the variants were further filtered using two additional criteria. (i) To ensure a high level of confidence in the existence of the variant, we require the QUAL score, which represents the (Phred-scaled) probability that all observed samples are homozygous reference alleles, should be ≥50. (ii) To ensure a high level of confidence that the variant allele is indeed homozygous, the genotype call (the “GT” value) must be a homozygous alternative call, and that the Phred-scaled genotype likelihood (the “PL” value) for the homozygous alternative call is at least 20 units above that of any other PL value. If the variant allele call for a strain does not meet these criteria, the call is converted to an “N” (undetermined) instead of being disregarded, which ensures that we don’t mistakenly assume that it is the reference allele. Information about the fold-coverage and other aspects of the sequencing methods are provided in S1 Text.

Haplotype Block Construction and Genetic Mapping in Mice HBCGM utilizes bi-allelic SNPs that are polymorphic among the strains. Only bi-allelic SNPs with at least one definitive homozygous alternative call among the 23 strains were included, while all other variants (including INDELs) were marked and excluded from HBCGM analysis. The locations and potential codon-changes caused by the SNPs are then annotated relative to the encoded genes using predictive gene models from Ensembl version 65. Only SNPs meeting the following criteria were used for haplotype block construction: (i) polymorphic among the strains with input trait data; and (ii) there were at least eight strains (or one-half of the number of input strains) with unambiguous allele calls. Haplotype blocks with 2, 3, 4, or haplotypes were then dynamically produced as described in [17], and the correlation between the input phenotypic data and the haplotype pattern within each identified block was evaluated as described [14]. The genes are then sorted based upon the ANOVA p-value (in increasing order) for numeric data or by the F statistic (in decreasing order) for categorical data. For a gene that is covered by multiple blocks, the smallest p-value (or the largest F statistic) obtained for all blocks within that gene was used. The genetic effect size is calculated using our previously described method [18]: Genetic effect size = (SSB − (k − 1) * MSE)/(SStotal + MSE), where SSB is the between-group sum-of-squares of the ANOVA model, SStotal is the total sum-of-squares, k is the number of groups, and MSE is within-group sum-of-squares divided by (n − k), where n is the total number of objects in the model. Given the very large number of haplotype blocks produced by the HBCGM method, all results were automatically filtered according to pre-determined criteria using software that selected the correlated genes that were expressed in a designated target organ and had missense SNPs.

Robustness Assessment To assess the robustness of the HBCGM output obtained using the data obtained from different sets of input strains, we developed software that would automatically perform a re-sampling analysis. For each iteration, a number (1, 2, or 3) was randomly selected to indicate the number of strains whose data would be randomly excluded from the 17 strains with available latency data. Next, the selected number of strains (and each of the individual strains to be deleted were also randomly selected) was excluded for this iteration. Then, HBCGM were run on this re-sampled data and the result was recorded (using p = 0.01 as the cutoff). This process was repeated 100 times. If a re-sampled strain panel was identical to one that had already been re-sampled, this iteration was disregarded, and another subset was randomly chosen using the same procedure. To measure the overall association of a gene with the phenotype among the 100 random re-sampled sets of output, we defined a robustness score based on the p-values obtained, which was calculated as: , where p i was the resulting p-value for that gene in the ith re-sampled experiment (if the p-value was bigger than the cutoff and therefore censored in the result, p i was set to 0.01, i.e., the cutoff). Genes with a larger robustness index have stronger overall correlation among the different subsets analyzed. Their calculated scores then were used to rank the genes, and those with the highest scores were output.

Statistical Analysis of Haloperidol Latency Data in Mice The haloperidol-induced latency data measured on days 0, 3, 7, 30, 60, and 120 (MPD datasets: 39407, 39408, 39409, 39410, 39411, and 39445, respectively) and the plasma haloperidol concentrations measured on days 30, 60, 90, and 120 (MPD numbers: 39403–39406, respectively) were obtained from the Mouse PHENOME Database (http://phenome.jax.org/). The correlation analyses were performed using the phenotypic data for all 27 strains that were evaluated. The non-parametric Spearman’s correlation coefficient (rho) was calculated [19], and the statistical significance of the correlation between the tested datasets was evaluated using this result. This non-parametric test was used to avoid producing spurious results that could be caused by the non-normal distribution of the data and by the presence of extreme outlying values. The p-values were then adjusted for multiple testing using the Benjamin-Hochberg method [20]. To minimize the effect produced by the strain with an extreme outlying value and to better visualize the relationship of the variables, the haloperidol-induced latency data from day 7 and day 30 were log-transformed (10-based). (It should be noted that the log-transformation has no effect on the Spearman’s rho, nor did it affect the corresponding test result.)

Haloperidol Toxicity Measurements All animal experiments were performed according to protocols approved by the Stanford Institutional Animal Care and Use Committee. All mice were obtained from Jackson Labs, and were used at 10–12 weeks of age, and the results are reported according to the ARRIVE guidelines (S1 Text) [21]. Haloperidol was administered to the mice by one of two methods: (i) a haloperidol pellet (Innovative Research of America) was subcutaneously implanted, which continuously released a dosage equivalent to 3 mg/kg/day, using published methods [22]; or (ii) with 10 mg/kg haloperidol (Sigma) IP qd for measurement of drug and metabolite concentrations. The haloperidol-induced latency was measured as the time required for a mouse to move all four paws after being placed on a vertical metal mesh screen as previously described [22].

Haloperidol Quantitation Haloperidol was purchased from Sigma-Aldrich; haloperidol-d4 was purchased from Toronto Research Chemicals Inc.; HPLC grade acetonitrile and water were from Honeywell Burdick and Jackson International Inc.; and formic acid was purchased from Pierce Thermo Scientific. Mouse brain tissues of known weight were homogenized in methanol using a Precelly 24 (Bertin Technologies) homogenizer. 50 μl of homogenate was aliquoted, and haloperidol-d4 was added as internal standard to a final concentration of 100 ng/ml. After the mixture was incubated with three volumes of cold acetonitrile at −20°, it was centrifuged at 14,000 rpm for 10 minutes. The supernatant was dried in a speed vacuum and re-suspended in the original volume. A standard curve was prepared using brain homogenate from an untreated mouse in a linear range from 0.5 ng/ml to 250 ng/ml. Quantitative analysis of haloperidol and its metabolite (HPP+) was performed on an Agilent accurate mass QTOF 6520A coupled with an UHPLC Infinity 1290. The analytes were separated on a Phenomenex Kinetex XB-C18 column 2.1 × 100 mm. The flow rate was 0.5 ml/min with a gradient of solvent B (acetonitrile with 0.1% formic acid; solvent A is 0.1% formic acid) from 5% to 35% in 10 minutes, then 35% to 95% in 5 minutes. Positive electrospray full scan in the range of m/z 110–1,000 spectra were collected. Data analysis was done using Agilent quantitative analysis software.

Desorption Electrospray Ionization Mass Spectrometroscopic Imaging High mass resolution/high mass accuracy mass spectrometry tissue imaging was performed using a lab-built desorption electrospray ionization mass spectrometry imaging (DESI-MSI) source coupled to an LTQ-Orbitrap XL mass spectrometer (Thermo Scientific). DESI-MSI was performed in the positive ion mode from m/z 200–500, using the orbitrap as the mass analyzer at a resolving power of 60,000. The spatial resolution of the imaging experiments was of 150 μm. A histologically compatible solvent system dimethylformamide:acetonitrile (DMF:ACN) 1:1 (v/v) was used for analysis [23] at a flow rate of 0.8 μl/min. The N 2 pressure was set to 175 psi. After DESI-MSI, the same tissue section was subjected to hematoxylin–eosin (HE) staining. The software ImgGenerator (freeware, http://www.msimaging.net/) was used for converting raw files into 2-D images. Spatially accurate ion images were assembled using BioMap software (freeware, http://www.maldi-msi.org/). The protonated form of haloperidol was detected in the tissue sections at m/z 376.14679 (mass error of −1.7 ppm). The isotopic distribution of the ion at m/z 376.14679 was also found to agree with the molecular formula of haloperidol. DESI-MS ion images of haloperidol were compared to optical images of the same tissue in HE-stained tissue sections. The total abundance of haloperidol in each pixel from a DESI-MS ion image was extracted and averaged to obtain the average abundance of haloperidol per tissue section. This procedure was performed for each tissue section analyzed for each strain of mice. Based on the average total abundance of haloperidol for each strain, the fold change was then calculated.

RNA In Situ Hybridization Abcb5 and Drd2 probes were designed to hybridize to their corresponding mRNAs. The 741 bp Abcb5 probe was PCR amplified from RIKEN_cDNA_clone B020023O04 using the primer set: AGGCCAGAAACAGAGGATTG; and CCTGGTAGAGCATGGCTTTG. The 1,126 bp Drd2 probe was PCR amplified from total cDNA from a C57BL/6J mouse using the primer set: CCGTGAACCCCATCATCTAT and GGTTTGGTGCATGTATGGTG. PCR reactions were performed using Phusion High-fidelity DNA polymerase (Thermo Scientific) and the products were cloned into pCR-BluntII-TOPO vector (Life Technology). Clones with correct insertions were linearized using either BamHI or NotI (New England Biolabs) for in vitro transcription. Anti-sense and sense RNA probes were then in vitro transcribed using Riboprobe System-SP6/T7 (Promega) and DIG RNA labeling mix (Roche), DNase treated, precipitated in lithium chloride (Ambion), and stored at −80°C in aliquots. Freshly dissected mouse brain was immediately bisected sagittally, and fixed in 4% paraformaldehyde in DEPC treated phosphate buffered saline at 4°C overnight. Fixed tissues were then saturated in DEPC-treated 20% sucrose at 4°C overnight, then flash frozen in OCT compound (Tissue-Tek) and stored at −80°C. The brain sections were cut at 20 μm thickness, and serial sections were collected on Superfrost Plus glass slides (Fisher Scientific). Sets of serial sections across the whole brain were then hybridized anti-sense probes, and adjacent sections were hybridized with the sense probe (as negative controls), at 63°C in a moist chamber overnight. After washing and blocking, tissue sections were then incubated with AP-conjugated anti-DIG antibody (Roche) at 4°C overnight. Tissue sections were then washed and chromogenic reactions were then performed by incubation with the AP substrate BM Purple (Roche) at room temperature. The sections were then examined using a light microscope, and a purple/blue color indicates the localization of the targeted mRNA within the tissue.

Statistical Analysis of Data Obtained from Inbred and Chromosome Substitution Strains The latency measurements from chromosome substitution strain (CSS) mice were log-transformed, and a two-sample t-test was applied to compare the latency of CSS12 to that of C57BL/6. The differences between the latency on the log-scale and the corresponding 95% confidence intervals were evaluated, and then transformed back to the original scale. Similarly, the measured haloperidol and HPP+ concentrations were log-transformed and a two-sample t-test was applied to compare the concentration in A/J to that in C57BL/6. The differences between the latency on the log-scale and the corresponding 95% confidence intervals were evaluated and then transformed back to the original scale. For comparison of the HPP+ levels in A/J, CSS12, and C57BL/6 strains, the concentrations were log-transformed and an F-test was applied. The differences between the latency on the log-scale and the corresponding 95% confidence intervals were evaluated using SAS (version 9.3), and were then transformed back to the original scale.

Analysis of Abcb5 Knockout Mice Mice with a homozygous Abcb5 gene deletion were prepared and bred onto the C57BL/6 by nine backcrosses, and these mice did not express Abcb5 mRNA (JPG, MMG, and colleagues, personal communication). The latency in control wild-type littermate and Abcb5 knockout mice were analyzed after 0 to 7 days of haloperidol 10 mg/kg/day IP administration. The statistical significance of the difference in the latency measurements between wild-type and Abcb5 knockout mice was determined using a two-sample two-sided t-test on log-transformed data for each day as described in the preceding paragraph. Since gender could affect the haloperidol response, we performed a two-way ANOVA to assess the effect of gender. However, the gender effect was only significant on day 4 (p = 0.046) and was not significant for all of the other days. Therefore, an adjustment for the gender effect was not performed.