Genetic Modulation of Initial Sensitivity to Δ9-Tetrahydrocannabinol (THC) Among the BXD Family of Mice.

Cannabinoid receptor 1 activation by the major psychoactive component in cannabis, Δ9-tetrahydrocannabinol (THC), produces motor impairments, hypothermia, and analgesia upon acute exposure. In previous work, we demonstrated significant sex and strain differences in acute responses to THC following administration of a single dose (10 mg/kg, i.p.) in C57BL/6J (B6) and DBA/2J (D2) inbred mice. To determine the extent to which these differences are heritable, we quantified acute responses to a single dose of THC (10 mg/kg, i.p.) in males and females from 20 members of the BXD family of inbred strains derived by crossing and inbreeding B6 and D2 mice. Acute THC responses (initial sensitivity) were quantified as changes from baseline for: 1. spontaneous activity in the open field (mobility), 2. body temperature (hypothermia), and 3. tail withdrawal latency to a thermal stimulus (antinociception). Initial sensitivity to the immobilizing, hypothermic, and antinociceptive effects of THC varied substantially across the BXD family. Heritability was highest for mobility and hypothermia traits, indicating that segregating genetic variants modulate initial sensitivity to THC. We identified genomic loci and candidate genes, including Ndufs2, Scp2, Rps6kb1 or P70S6K, Pde4d, and Pten, that may control variation in THC initial sensitivity. We also detected strong correlations between initial responses to THC and legacy phenotypes related to intake or response to other drugs of abuse (cocaine, ethanol, and morphine). Our study demonstrates the feasibility of mapping genes and variants modulating THC responses in the BXDs to systematically define biological processes and liabilities associated with drug use and abuse.

Introduction

Recombinant inbred (RI) rodent populations are a valuable resource for forward genetic mapping and systems genetics analysis. Individual RI lines are stable and each genotype can be resampled to boost mapping power and improve quantitative trait loci (QTL) detection (Phillips et al., 1991). In addition, RI population genotypes and trait data collected across time and laboratory settings can be combined into a powerful database of legacy phenotypes amenable to systems scale analyses both across and between trait levels. Here, we leverage the BXD RI family of strains originally derived by crossing C57BL/6J (B6) and DBA/2J (D2) mice followed by intercrossing and inbreeding of their progeny to yield stable RI lines (Figure 1A). Since its inception in the late 1970s, the BXD RI panel continues to be a valuable resource for characterization of the genetic architecture of addiction-related behavior and identification of genes and variants modulating the response to drugs of abuse (Phillips et al., 1991; Rodriguez et al., 1994, 1995; Hitzemann et al., 2003; Philip et al., 2010; Jackson et al., 2011; Dickson et al., 2016).

FIGURE 1

Experimental Overview. (A) Overview of the generation of the BXD RI population. This family of strains has been derived at several points in time from both F2 and advanced intercrosses and there are now ∼150 strains available. (B) Description of daily treatment regime. Initial sensitivity is measured by subtracting Day 1 trait values from those of Day 0 for each individual. (C) Timing of daily trait measurements for body temperature (Temp), time mobile in the open field (OF), and latency to withdrawal the tail in response to a thermal stimulus (TF) are shown as minutes post-injection.

With well over 100 BXD lines currently available (Ashbrook et al., 2021), the BXD RI family is the largest and best characterized rodent genetic reference population available (Mulligan et al., 2017; Parker et al., 2017; Ashbrook et al., 2021). This family contains over 6 million segregating variants and ∼10,000 recombinations, resulting in comparatively high mapping precision (Ashbrook et al., 2021). The BXD family has also been deeply phenotyped for molecular, behavioral, physiological, and pharmacological traits since the late 1970’s and includes data related to cognitive function, anxiety and stress, social interactions, and response to drugs of abuse (Philip et al., 2010; Laughlin et al., 2011; Newbury and Rosen, 2012; Harenza et al., 2014; Loos et al., 2014; Neuner et al., 2016; Ashbrook et al., 2018; Knoll et al., 2018).

Previously, we established that the parents of the BXDs (strains B6 and D2) differ greatly in their initial responses to a single dose of THC (10 mg/kg, i.p.), the major psychoactive component in cultivars of cannabis (Parks et al., 2020). Activation of the cannabinoid 1 (CB1) receptor in discrete neuronal populations (Monory et al., 2007) by THC results in a well-characterized trait spectrum that includes a reduction in spontaneous locomotion, hypothermia, and antinociception. Genetic or pharmacological deletion of the CB1 receptor abrogates these effects (Rinaldi-Carmona et al., 1994; Compton et al., 1996; Ledent et al., 1999; Zimmer et al., 1999; Huestis et al., 2001). Thus, quantification of mobility, hypothermia, and antinociception following THC exposure is a robust method to measure response to THC mediated specifically through CB1 receptor signaling pathways. In our previous work, we demonstrated that B6 mice exhibited a greater reduction in THC-induced spontaneous locomotion relative to D2 mice and that D2 mice exhibited a greater reduction in THC-induced antinociception relative to B6 mice. We detected minimal strain differences in THC-induced hypothermia, however, females were much more sensitive to this effect than males. Moreover, our previous findings demonstrated that differences between B6 and D2 in THC responses are caused by pre-existing genetic variation in effectors of endogenous cannabinoid receptor signaling.

Remarkably, genetic variation at the CB1 receptor locus (Cnr1) is not the cause of variation in THC responses between B6 and D2. In other previous work, we quantified striatal protein levels of the CB1 receptor and found them to be significantly higher in B6 relative to D2 mice (Parks et al., 2019). There are no high impact segregating variants within the Cnr1 locus of B6 and D2 strains that could account for the observed variation in CB1 receptor levels and there was also no evidence of genetic modulation of Cnr1 transcript expression levels among BXD strains (Parks et al., 2019). Thus, the causal genes and variants in the CB1 signaling pathway that contribute to genetic differences in response to THC among B6, D2, and the BXDs remain to be discovered.

Identification of genetic variation in CB1 receptor signaling pathways has important behavioral and physiological implications that extend beyond cannabis and cannabinoids. This system plays an important role in modulating synaptic plasticity, behavior, and reward circuitry. Moreover, disruptions of endocannabinoid system function and signaling contribute broadly to learning and addiction processes for both cannabinoid and non-cannabinoid drugs (Parsons and Hurd, 2015). Recent studies quantified initial response to several drugs of abuse (i.e., cocaine, alcohol, morphine) or intravenous cocaine self-administration among a large set of BXDs (Philip et al., 2010; Dickson et al., 2016). BXD traits associated with these studies are available in an open data repository (GN)[1] and are of sufficient strain depth to facilitate well-powered correlation analysis between initial responses to THC and initial response or self-administration of other drugs of abuse. The existence of strong correlations between cannabinoid and non-cannabinoid drug response traits is evidence of shared genetic and/or biological co-regulation, potentially due to genetic variation in CB1 receptor signaling pathways.

In this study, we leverage the BXD family to address two main questions. Our first question is whether B6 and D2 differences in THC responses are heritable and, if so, segregate amongst the BXD family. Our second question is whether initial responses to THC among the BXD population are related to responses to other drugs of abuse. Answers in the affirmative would greatly support the use of the BXD family to identify DNA variants and genes in cannabinoid signaling pathways and shared causal networks mediating responses to THC and other abused substances.

Materials and Methods

Experimental Subjects

The BXD family (Figure 1A) is maintained in a large breeding colony at The University of Tennessee Health Science Center. BXD mice were made available for this study as part of a Mouse Strain and Pilot Projects program. Mice from up to 20 BXD strains were provided after weaning (∼21 days) as they became available in the colony. Mice were transported in same sex cages housed with littermates across campus from the breeding colony to the housing and testing facility in an air conditioned and covered vehicle by the UTHSC Lab Animal Care Unit several weeks prior to any handling. At least 1 week prior to testing, animals were housed individually and handled daily. Handling consisted of lifting each individual in either cupped hands (avoiding any lifting by the tail) or a plastic lid from a pipette tip box. Food (Envigo 7921 standard feed) and water were provided ad libitum and animals were maintained on a 12 h:12 h light:dark cycle. Shepard’s cob was provided as standard cage bedding. All testing was performed during the light cycle from 0700 to 1600 h. All procedures were approved by The University of Tennessee Health Science Center Institutional Animal Care and Use Committee.

The experimental design of the study called for inclusion of three males and three females from each BXD strain from one to two litters and a target age at testing of 70–150 days-of-age (adult). A summary of actual ages, sex, strains, and litter numbers obtained for the study can be found in Supplementary Table 1. We were able to obtain 2–5 males and females for all strains except BXD1, BXD11, and BXD15, for which, no female data was obtained and a single male was tested for BXD11. Average age for BXDs included in the study was 115 days with a standard deviation of ± 16 days. Age ranges from 57 to 155 days with an average range between 99 and 127 days. Strain and sex groupings consist of a mixture of target age ranges (Supplementary Figure 1), with the exception of female BXD77 females (all aged 152 days) and BXD75 males (aged 153 and 155 days).

Three strategies were incorporated into the study to control for batch effects due to the somewhat arbitrary availability of BXD subjects. First, strain and sex were counterbalanced across batches when possible. Based on availability from the BXD breeding colony, mice were assigned to testing cohorts consisting of a maximum of 32 mice tested each week. Within each cohort, mice were assigned to one of up to four batches. Each batch consisted of up to eight mice that were phenotyped together within the same test session. A summary of cohort and batch assignments for each sex and strain can be found in Supplementary Table 1. Second, instead of using a matched control group, we quantified treatment response relative to each subject’s own baseline (see description of difference scores as a method to calculate Initial Responses to THC below). Previously we found quantification of difference scores to be a robust and more direct measurement of drug response that is more resistant to batch effects that may differentially impact control and experimental groups (Mulligan et al., 2018; Parks et al., 2020). Finally, for inbred panels, strain mean responses are used for QTL mapping as opposed to mapping trait data at the level of the individual (see section “Materials and Methods” for “Genetic Mapping” below).

THC and Vehicle Formulation

Δ9-Tetrahydrocannabinol was formulated in an ethanol:cremophor:saline (5:5:90) vehicle followed by filter sterilization. The resulting formulation was stored in the dark at 4°C in a septum sealed vial. Vehicle (VEH) was prepared in the same manner. THC and VEH were administered by intraperitoneal (i.p.) injection at a dose of 10 mg/kg such that a 30 g mouse received a 100 μl injection volume. Each formulation was used for the entire testing schedule of each cohort (see details below). Previously we demonstrated that each THC formulation contains 95.2% of the initial amount of THC when measured 8 days after preparation (Parks et al., 2020).

Testing Schedule

Treatment and trait measurements are described in Figures 1B,C and in Parks et al. (2020). BXD mice of each sex and strain were included (see Supplementary Table 1 for exact numbers). On the first day (Day −1), no injections were given. On the second day (Day 0, baseline), all mice received an injection of VEH (100 μl per 30 g). On Day 1, mice received an injection of 10 mg/kg THC. On each day, body temperature, tail withdrawal latency in response to a thermal stimulus, and activity (time mobile) in the open field (OF) were measured at multiple time points in the same animal. Rectal body temperature was measured using a ThermoWorks digital thermometer with a mouse rectal probe adaptor at time 0 and 60 min post-injection. Tail withdrawal latency was measured 60 min post-injection by gently restraining each mouse in a 50 mL conical tube. The tail was then submerged ∼2 cm into a 52°C water bath. The latency to remove the tail was recorded to the nearest second. In most cases, tail flick latency was measured by two independent experimenters, and the resulting scores were averaged. Time spent mobile or immobile (periods of no movement lasting for 3 or more seconds) was measured 30 min post-injection in a 40 cm × 40 cm × 40 cm OF over a 10 min interval using video recording and ANY-maze (Stoelting) tracking software.

Statistical Analysis of THC Treatment Effects

BXD response data from day 0 (VEH treatment) to day 1 (THC treatment) were used in a three-way between-subjects (treatment × strain × sex) analysis of variance (ANOVA). The purpose of this analysis was twofold: we wanted to determine if THC treatment had a large and significant effect on responses and quantify any interactions between treatment and sex that would justify additional sex-specific analyses. ANOVA was performed using base functions in R (anova.lm function). Effect sizes (partial omega-squared or ωp2) for each parameter in the ANOVA were calculated in R using the omega_squared function in the effectsize package. Response during habituation on day −1 was not analyzed and individual BXD data points were only removed for technical reasons (e.g., equipment failure during recording).

Statistical Analysis of Initial Responses to THC

Because THC treatment effects were significant and of moderate to large effect for all traits, initial responses to THC were calculated as the difference between day 1 and day 0 (baseline) for each trait (Figure 1B). Difference scores provide a more direct method to demonstrate initial THC responses and are more robust to the effects of batch (note that counterbalancing was not possible in the study design). The low number of replicates among BXD strains precluded outlier analysis (see Supplementary Table 1). BXD male and female responses were combined for the antinociception and mobility initial sensitivity traits and separated for the hypothermia initial sensitivity trait based on detection of sex interaction effects between treatment and strain for hypothermia in the three-way ANOVA to evaluate THC treatment effects. Strain averaged initial sensitivity BXD traits for mobility and antinociception (Record IDs 21481 and 21483, respectively) and strain and sex averaged initial sensitivity for the BXD hypothermia trait (Record IDs 21485 and 21487) have been uploaded to the GeneNetwork (GN2)[2] web service (RRID:SCR_002388) for the Mouse species, the BXD Family group, the Traits and Cofactors type, and the BXD Published Phenotypes dataset.

To quantify and visualize the effect of strain on initial response to THC for each trait (y), one-way between-subjects ANOVA in the form of y ∼ strain was performed using base functions in R (anova.lm function). Effect sizes (omega-squared or ω2) for each ANOVA were calculated in R using the omega_squared function in the effectsize package. Heritability (h2), an estimate of the proportion of the trait variation explained by genetic factors, can be calculated based on the variance components of a one-way ANOVA (y ∼ Strain). We used ω2 as an estimate of h2.

Initial sensitivity to THC in the parental B6 and D2 strains has been profiled extensively using identical phenotyping and analysis methods (Parks et al., 2020). For comparison to BXDs, summarized parental data is included in Figure 2. Parental data includes THC initial responses quantified in 17 male and nine female B6 mice and 17 male and six female D2 mice. The age range of parental mice overlaps that of the current BXD study.

FIGURE 2

Variation in Initial Response to THC Among BXDs. A significant effect of THC treatment relative to baseline (p < 0.001) was observed for all traits. Initial response to THC is shown as the difference between baseline and initial THC treatment. Strains are shown on the x-axis and initial response (difference between baseline response on day 0 and response to THC on day 1) is shown on the y-axis. Negative values indicate a decrease in response on day 1 compared to day 0 and positive values indicate an increase in response. Parental strains are indicated in blue (B6) and red (D2) and BXD strains are shown in black. (A) There is a significant effect of strain (p < 0.001) on initial response to THC for time spent mobile in the open field 30 min post-injection. Male and female responses were combined as there were no significant interaction effects involving sex. Responses to THC varied ∼90-fold among strains for the mobility trait. (B) The antinociceptive effects of THC were quantified using tail withdrawal latency in response to a thermal stimulus at 60 min post-injection. Male and female responses were combined as there were no significant interaction effects involving sex. Response to THC varied nearly 8-fold across strains, however, no significant effects of strain were observed for the antinociception trait. (C,D) There is a significant main effect of strain (p < 0.05) on initial hypothermic response 60 min post-injection of THC for both females (C) and males (D). Hypothermic response to THC in females varies 13-fold across strains compared to 24-fold variation in response in males. Summarized B6 and D2 responses shown for comparison (Parks et al., 2020).

Genetic Mapping

Genetic mapping is performed using strain means and not individual level data. Power in genetic mapping studies is primarily determined by the number of unique individual genomes. Recombinant inbred panels allow for replication of genomes within strain which can further increase power to detect QTLs. Here we applied two different methods to compute QTL probability given strain genotypes and initial responses to THC using GN (Mulligan et al., 2017). The first method applied was traditional simple regression (Haley-Knott or HK) and the second was genome-wide efficient mixed model association [GEMMA; (Zhou and Stephens, 2012)] with the “leave one chromosome out” (LOCO) option [rationale for using both methods described in Mulligan et al. (2018)]. Prior to performing genome scans, the strain-averaged distribution for each trait was checked for normality using the probability plot function in GN. All traits were approximately normally distributed and used for QTL mapping. Note that only strain-averaged BXD data is informative and used for mapping while parental strain data (i.e., B6 and D2) are not included. For both HK and GEMMA, a dense panel of 7,321 markers was used for mapping. QTL were considered noteworthy if: (1) a genome-wide suggestive level (p < 0.63, equivalent to a 63% probability of a false positive or one false positive per genome scan) was reached following permutation (1000 tests) using the HK method, and (2) the peak marker association [−log(p) > 3] detected using GEMMA overlapped the QTL mapped by HK.

Genome-wide efficient mixed model association permutation was performed for each trait as described previously (Mulligan et al., 2018). For trait 21481, the 95th percentile (significant) threshold and 67th percentile (suggestive) threshold were 4.8 and 3.8 LOD, respectively. For trait 21483, the 95th percentile (significant) threshold and 67th percentile (suggestive) threshold were 4.5 and 3.6 LOD, respectively. For trait 21485, the 95th percentile (significant) threshold and 67th percentile (suggestive) threshold were 4.7 and 3.7 LOD, respectively. For trait 21487, the 95th percentile (significant) threshold and 67th percentile (suggestive) threshold were 4.7 and 3.7 LOD, respectively.

Quantitative trait loci scans on small subsets of BXDs, as in our study, are expected only to identify QTLs of large effect. For example, a mapping population comprised of 20 BXD strains with four replicates per strain is only powered to detect a QTL explaining 50% of the trait variance, whereas a population size of ∼80 strains with four replicates per strain is well powered to detect a QTL explaining 20% of the trait variation at 80% power (Andreux et al., 2012; Ashbrook et al., 2021).

Candidate Gene Search

A 1.5-LOD drop from the peak marker (HK mapping) was used to define an approximate ∼95% confidence interval (CI) for a QTL of interest (Alonso-Blanco et al., 2006). For each QTL CI, a complete list of mouse reference genes was generated using the UCSC Genome Browser (RRID:SCR_005780) Table Browser Tool (Group: Genes and Gene Predictions; Track: NCBI RefSeq; Genome: GRCm38/mm10).

Genes in each QTL CI were then prioritized based on several criteria. First, all variants distinguishing B6 and D2 within each QTL CI were identified using Sanger’s Mouse Genome Project Mouse_SNPViewer/rel-15050 (RRID:SCR_011784). Specifically, SNPs, InDels, and structural variants (SVs) distinguishing B6 and D2 and their predicted consequences were retrieved and compared (by gene symbol) to reference genes located within the boundaries of each QTL CI. Next, legacy mRNA datasets available in GN for BXD brain tissue were queried to identify genes located in each QTL CI whose expression was modulated by variants located within or near the location of the gene itself – a cis expression QTL (eQTL). Brain region data sets (cortex, hippocampus, hypothalamus, nucleus accumbens, and striatum) were selected based on known functional regulation by cannabinoid receptor signaling, moderate to high expression of CB1 receptor in each region, and functional involvement of each region in initial response traits (e.g., motor, hypothermia, and nociception). Datasets included: Hippocampus Consortium M430v2 (Jun06) (GN Accession: GN110 and GEO Series: GSE84767); INIA Hypothalamus Affy MoGene 1.0 ST (Nov10) (GN Accession: GN281 and GEO Series: GSE36674); HQF BXD Neocortex ILM6v1.1 (Feb08) RankInv (GN Accession: GN284); VCU BXD NA Et vs. Sal M430 2.0 (Oct07) (GN Accession: GN156); VCU BXD PFC Sal M430 2.0 (Dec06) RMA (GN Accession: GN135; GEO Series: GSE28515); HQF Striatum Affy Mouse Exon 1.0 ST Gene Level (Dec09) RMA (GN Accession: GN399). Cis eQTL data are summarized in Supplementary Table 2. Probes that overlap polymorphic SNPs and have cis eQTLs with higher expression associated with the B6 (B) allele (ILM6v1.1 and M430 data sets only) are flagged in this table as expression measurements that could be biased in the direction of B alleles (higher expression relative to D2, or D, alleles) due to technical probe hybridization artifacts. Finally, a search for relevant biological function (i.e., endocannabinoid/cannabinoid receptor signaling) was conducted using the GeneCup (genecup.org; formerly known as RatsPub)[3] web service (RRID:SCR_021281) to search through PubMed (RRID:SCR_004846) literature abstracts for associations between cannabinoid-related terms and gene symbols for each QTL CI. A secondary search was performed for high priority candidate genes using GeneCup to place these genes in a broader context of addiction and psychiatric diseases.

BXD Exploratory Phenome Analysis

Two recent studies quantified initial responses to multiple drugs of abuse (cocaine, alcohol, and morphine) or intravenous cocaine self-administration among BXDs (Philip et al., 2010; Dickson et al., 2016). Both studies also included large numbers of BXDs, including at least 12 of the strains used to generate our THC initial response data. Associated drug response traits from these studies have been deposited into the GN BXD Published Phenotypes dataset as strain-averaged trait data. We selected BXD drug response and cocaine self-administration traits from both of these studies for comparison with THC initial response traits generated by our study. Trait data was obtained with the Get Any search function in GN using the following input: 2681503 19958391 (PubMed IDs). This resulted in 790 traits. The BXD legacy traits mined from GN were retained for correlation analysis if there were at least 12 matched BXD strains shared with our THC initial response traits. This resulted in retention of 757/790 traits. Traits were further filtered if descriptions included the terms “Morphine response,” “morphine withdrawal,” “Cocaine response,” “Intravenous cocaine,” OR “Ethanol response” and DID NOT include the terms “Saline control,” “inactive,” OR “unconditioned.” A final total of 302 traits (excluding the four initial THC response traits) were included for correlation analysis. Traits were assigned to general categories based on phenotype and whether the trait was measured in females, males, or both. Category descriptions along with GN record IDs, trait descriptions, and values for drug response and cocaine self-administration can be found in Supplementary Table 3. The Pearson correlation coefficient and corresponding p-values were calculated for each of the four THC initial response traits and all 302 BXD legacy drug response traits. Significant correlations were associated with a p < 0.05. Note that all correlations were performed on BXD strain averaged trait data and do not include parental B6 and D2 strains.

Results

Acute THC Treatment Has a Profound and Significant Effect on Mobility, Hypothermia, and Antinociception Initial Response Traits Among the BXD Family

We found that treatment (i.e., baseline values on day 0 compared to day 1 values following THC treatment in the same individual) had significant (all p-values < 0.001) and moderate to large main effects on mobility, hypothermia, and antinociception in our BXD cohort. Relative to the VEH treatment on day 0, THC treatment on day 1 greatly and significantly reduced mobility (F1,142 = 417.50, p < 0.001, ωp2 = 0.66, 90% CI [0.59, 0.71]). THC treatment caused a large and significant decrease in body temperature (F1,144 = 398.50, p < 0.001, ωp2 = 0.65, 90% CI [0.57, 0.70]). Acute THC treatment resulted in a moderate and significant increase in tail withdrawal latency (F1,141 = 88.75, p < 0.001, ωp2 = 0.29, 90% CI [0.19, 0.38]).

Relative to treatment effects, main effects of strain or sex were generally smaller. Strain effects contributed significantly to response variation for mobility (F19,142 = 6.13, p < 0.001, ωp2 = 0.31, 90% CI [0.13, 0.35]), hypothermia (F19,144 = 3.42, p < 0.001, ωp2 = 0.17, 90% CI [0.00, 0.18]), and antinociception (F19,141 = 2.38, p < 0.01, ωp2 = 0.11, 90% CI [0.00, 0.08]). Significant main effects of sex were observed for hypothermia (F1,144 = 5.43, p < 0.05, ωp2 = 0.02, 90% CI [0.00, 0.07]) and antinociception (F1,141 = 5.40, p < 0.05, ωp2 = 0.02, 90% CI [0.00, 0.07]).

We observed a moderate and significant treatment-by-strain interaction effect for mobility (F19,142 = 4.90, p < 0.001, ωp2 = 0.26, 90% CI [0.07, 0.29]) and a smaller significant treatment-by-strain interaction effect for hypothermia (F19,144 = 3.04, p < 0.001, ωp2 = 0.15, 90% CI [0.00, 0.05]). Smaller effect trends for sex-by-treatment (F1,144 = 2.62, p = 0.11, ωp2 = 0.007, 90% CI [0.00, 0.05]), strain-by-sex (F16,144 = 1.34, p = 0.18, ωp2 = 0.02, 90% CI [0.00, 0.00]), and treatment-by-strain-by-sex (F16,144 = 1.71, p = 0.05, ωp2 = 0.05, 90% CI [0.00, 0.00]) interactions were observed for the hypothermia trait.

Significant Strain Effects Account for Variation in Initial Responses to Acute THC Treatment Among BXDs

To demonstrate and visualize the direct effect of acute THC treatment in each BXD strain, we calculated the difference between day 1 and baseline day 0 as the initial THC response for each trait (Figure 1B). Based on the presence or absence of treatment-by-sex interaction effects, male and female initial THC responses were combined by strain for the mobility and antinociception traits and separated by sex for the hypothermia trait. Initial responses to THC measured in 20 BXD strains are summarized in Figure 2 (distribution within males and females of each strain can be found in Supplementary Figure 2).

For mobility, strain had a large and significant effect on initial response to THC (F19,87 = 5.87, p < 0.001, ω2 = 0.46, 90% CI [0.25, 0.52]). Average mobility across BXD strains was reduced by 238.40 ± 28.08 s (all values reported as mean ± standard error, Figure 2A). Strains varied an astonishing 90-fold in their initial response to the immobilizing effects of THC. Similar to the sensitive B6 parental strain (−437.30 ± 19.49 s), BXD16 (−439.70 ± 26.03 s), and BXD87 (−468.13 ± 24.69 s) demonstrated large reductions in mobility after acute THC. In contrast, BXD71 (−6.47 ± 6.47 s) and BXD83 (−5.06 ± 7.50 s) demonstrated lower sensitivity to THC, even when compared to the less sensitive D2 parental strain (−164.48 ± 21.14 s).

Acute THC resulted in an average 1.20 ± 0.11 s increase in tail withdrawal latency across all strains measured (Figure 2B). Initial analgesic/antinociceptive response to THC varied nearly 8-fold across strains with an increase in tail withdrawal latency from 2.60 s in the most sensitive strain (BXD14) to only 0.33 s in the least sensitive strain (BXD87). However, strain did not have a significant effect on initial analgesic response to THC (F19,86 = 1.31, p = 0.20, ω2 = 0.05, 90% CI [0.00, 0.00]).

Strain had a moderate and significant effect on initial hypothermic response to THC in females (F16,35 = 2.42, p < 0.05, ω2 = 0.30, 90% CI [0.00, 0.30]). In males the effect was also significant, but larger (F19,37 = 4.07, p < 0.001, ω2 = 0.51, 90% CI [0.03, 0.56]). Average body temperature across BXDs was reduced by 3.36 ± 0.33°C in females and by 2.42 ± 0.27°C in males. BXD62 (−6.07 ± 1.35°C) and BXD101 (−5.73 ± 0.60°C) females and BXD16 (−4.86 ± 0.35°C) males demonstrated greater relative sensitivity to the hypothermic effects of THC (Figure 2C). In contrast, BXD83 (−0.45 ± 0.25°C) females and BXD11 (0.20°C), BXD48a (−0.45 ± 0.30°C), and BXD74 (−0.57 ± 0.58°C) males demonstrated reduced sensitivity to the hypothermic effects of THC (Figure 2D). BXD females and males varied in initial response to the hypothermic effects of THC by 13-fold and 24-fold, respectively, and the range of responses in the BXDs exceeded that of the parental B6 and D2 strains.

Variation in Initial Responses to THC Are Heritable Among BXDs

Previously, we identified putative heritable strain differences in acute response to a single THC injection (10 mg/kg) in the parents of the BXD family – B6 and D2 (Parks et al., 2020). Here we test the hypothesis that polymorphic variants between these strains segregate among the BXD progeny and cause variation in THC initial responses. If responses to THC are variable and heritable among BXDs, then QTL mapping in the BXD population has the potential to identify causal variants in effectors of CB1 receptor signaling.

In general, THC response traits demonstrated high heritability (Table 1). Tail withdrawal latency was an exception. Mobility and hypothermia traits demonstrated higher heritabilities, indicative of a strong contribution of genetic variants to phenotypic variation.

TABLE 1

Estimates of heritability among BXD strains.

THC IS Trait	h2
Time Mobile (F + M)	0.46
Hypothermia (F)	0.30
Hypothermia (M)	0.51
Antinociception (F + M)	0.05

Initial Response to the Effects of THC on Mobility and Hypothermia in Males Are Correlated and Co-regulated by a Locus on Chr 11

Examination of trait correlation structure can provide evidence of underlying co-regulation of trait expression due to biological factors. Based on the presence or absence of treatment-by-sex interaction effects, male and female initial THC responses were combined by strain for the mobility and antinociception traits and separated by sex for the hypothermia trait resulting in a total of four traits for correlation analysis (Figure 3A). We observed a significant (p < 0.01) and positive correlation between time mobile (males and females combined) and hypothermia in males, but not in females, following the first exposure to THC (Figures 3A,D). As expected, based on the presence of a strong treatment-by-sex interaction effect of THC, hypothermic responses to THC in males and females were modestly and positively correlated, but these correlations were not significant. Modest and positive (but not significant) correlations were also observed between the hypothermic response in females and time mobile following THC treatment. In contrast, THC-induced antinociception was uncorrelated with both the mobility trait, and female and male hypothermia traits.

FIGURE 3

THC Initial Response Trait Correlations and Genetic Co-Regulation. (A) Trait correlation network constructed after calculating Pearson’s correlation on ranked trait values to generate Spearman’s rank order correlation coefficient (Rho or P). The hypothermia trait was split by sex due to a condition-by-sex interaction effect on body temperature change in response to THC. Changes in time mobile in the open field and body temperature in males following a single THC (10 mg/kg) exposure are significantly (p < 0.01) correlated. Weak (not significant) correlations are also observed between mobility and hypothermia in females following a single THC exposure and between changes in male and female hypothermia in response to THC. (B) Genetic regulation of THC initial response traits is observed, albeit at a suggestive level with no traits passing the threshold for genome-wide significance (p < 0.05). Interval maps for each trait shown as a different line color corresponding to the correlation network in (A). The strength of association (LOD) on the y-axis is plotted for each trait across the genome (megabase position on each chromosome or Chr) on the x-axis. Change in mobility and male body temperature are regulated by the same suggestive locus on Chr 11 (zoomed region in boxed area). (C) Use of an alternative linear mixed-model QTL mapping method (GEMMA, see section “Materials and Methods”) resulted in replication of genetic co-regulation of initial motor and hypothermia THC response traits from the same Chr 11 interval. Strength of association or –log(P) values (y-axis) shown for each marker (blue dots) within the Chr 11 QTL interval for mobility and temperature (males) initial THC response traits. Mapping with linear mixed-models can account for population family structure or kinship within the BXD population that is not addressed using traditional interval mapping. SNP density plotted in orange on the x-axis. (D) Scatterplot describing the relationship between average change in time mobile in the open field in response to THC (X-axis) and average change in male body temperature in response to THC (Y-axis) for each BXD strain based on Spearman’s rank order correlation.

Heritability estimates were high for several initial THC response traits (i.e., change in time mobile in the open field in response to THC; Table 1), indicating strong regulation of trait expression by genetic variants. To explore the potential for future genetic mapping studies of these traits in a larger BXD panel, we performed QTL scans for all traits. As expected, no genome-wide significant QTLs (LOD > ∼3.5) were detected for any trait (Figure 3B) using traditional interval HK mapping methods. However, suggestive QTLs (LOD > 2.5) were detected for all THC initial response traits (Table 2). Moreover, the highly correlated mobility and male hypothermia traits were co-regulated by a QTL on Chr 11 with a 1.5 LOD drop confidence interval between 80 and 90 Mb (Figure 3B).

TABLE 2

Summary of suggestive QTLs for THC initial response traits.

Trait	LOD	Chr	1.5 LOD CI (Mb)	GEMMA−log(p) (Peak Chr:Mb)
Time Mobile	2.8#	11	80,036,498–90,077,597	3.2 (Chr 11:89,612,498)
Hypothermia (M)	3.2# 2.5#	11; 19	80,036,498–89,612,498 27,531,578–36,395,869	3.9# (Chr 11:82,187,018) 3.7# (Chr 19:32,887,054)
Hypothermia (F)	2.6#	13	107,382,731–115,929,261	3.7# (Chr 13: 110,449,36)
Antinociception	2.9# 2.5#	1 4	168,144,357–172,235,364 102,851,020–109,169,766	3.9# (Chr 1: 170,949,421) 3.3 (Chr 4: 107,548,653) 4.3# (Chr 5: 53,611,846)

Shared genetic co-regulation of mobility and male hypothermia traits on Chr 11 was also replicated using a different QTL mapping method (GEMMA, see section “Materials and Methods”) that leverages linear mixed-models and accounts for population family structure (Figure 3C). Suggestive QTLs for initial response to the hypothermic (Chr 19 in males and Chr 13 in females) and analgesic (Chrs 1 and 4) effects of THC also replicated across both QTL mapping methods (Table 2). However, the mobility (Chr 11) and antinociception QTL (Chr 4) did not meet the suggestive threshold using GEMMA. Additional suggestive QTLs for antinociception (Chr 5) were detected by GEMMA only (Table 2).

Effectors of Endocannabinoid Signaling Identified as Positional Candidates for THC Initial Response QTLs

As a first step toward identification of variants and genes modulating initial responses to THC, we prioritized several suggestive QTLs for exploration of putative candidate genes. From Table 2, we selected the Chr 11 locus modulating both mobility and male hypothermia in response to THC and loci on Chrs 19 (male hypothermia), 13 (female hypothermia), Chr 1 (antinociception), and Chr 4 (antinociception). All QTLs were replicated using both HK and GEMMA mapping methods. Following a standard workflow (see section “Materials and Methods”) we identified a number of genes near or overlapping SNPs, InDels, or SVs, some of which are predicted to have an impact on transcript or protein integrity (e.g., missense, stop gain or loss, splice region variant, or frameshift) (Supplementary Table 4). Based on evidence in support of genetic regulation of gene expression by cis eQTLs in naive BXD brain tissue (Supplementary Table 2) and literature associations with relevant cannabinoid terms detected by GeneCup (see section “Materials and Methods”), we nominated several positional candidates for each prioritized QTL. Top positional candidates include Ndufs2 (Chr 1 antinociception QTL); Scp2 (Chr 4 antinociception QTL); Rps6kb1/p70s6K (Chr 11 mobility and male hypothermia QTL); Pde4d (Chr 13 female hypothermia QTL); and Pten (Chr 19 male hypothermia QTL) (Table 3 and Supplementary Table 4). These genes are likely to play a role in modulating initial response to THC based on location within modulatory QTL, putative functional sequence variants between B and D haplotypes that segregate among BXD progeny, evidence of genetic control of expression in brain tissue, and previously reported involvement in endocannabinoid/cannabinoid receptor signaling pathways. Below we summarize the evidence in support of these top positional candidates.

TABLE 3

Summary of genes with variants and high impact variants within QTL CIs.

QTL	Genes	SNP or InDel	SV	High Impact	cis eQTL	Priority Genes
Chr 11: Mobility and Hypothermia (M)	212	130	34	60	33	Rps6kb1/p70s6kC, AatfC, AcacaHCP. Akap1YCN, Ankfn1Y, Bcas3Y, Car4S, Ccl1, Ccl11, Ccl3, Ccl5, CltcC, Cuedc1HC, DgkeHYC, Dhrs11Y, Dhx40H, Dynll2Y, Gdpd1/Gde4, Ggnbp2C, Heatr6HYCN, Mks1Y, Mrm1C, Mrps23H, Msi2HC, Myo19Y, Ppm1dHN, Scpep1HY, Slfn8C, Trim25CN, Usp32HN, Wfdc18C
Chr 19: Hypothermia (M)	105	53	6	3	5	PtenN, 2700046G09RikHYCN, Atad1HN, Papss2HCN,
Chr 13: Hypothermia (F)	107	46	8	8	12	Pde4dCS, Ercc8C, Map3k1C, Plk2N, Plpp1P, Rab3cHC, Zswim6H
Chr 1: Antinociception	91	60	11	29	23	Ndufs2HYCNP, Adamts4HCN, Alyref2CN, Arhgap30H, Cd84CN, Fcer1gHYC, Fcgr3YCS, Klhdc9Y, Nit1HYCNP, Pfdn2HYNP, PpoxHCNPS, SdhcHCNPS, Ufc1HCNP, Usp21C
Chr 4: Antinociception	107	68	15	14	16	Scp2CP, Acot11HYN, Cc2d1bY, Kti12N, Mroh7YC, Ndc1HC, PodnC, Rab3b, Tceanc2YCS, Txndc12C, Zcchc11HN, Zfyve9HCN, Zyg11bHNP

NADH Dehydrogenase [Ubiquinone] Iron-Sulfur Protein 2 (Ndufs2) Is a Positional Candidate for Modulation of the Analgesic Response to THC

Ndufs2 encodes a core subunit of the mitochondrial Complex 1. The gene also contains both missense and putative splice region variants and demonstrates evidence of genetic regulation of expression in the form of cis eQTLs detected in naive BXD hippocampus, hypothalamus, neocortex, nucleus accumbens, and prefrontal cortex. Moreover, Ndufs2 has been implicated in mediating response to cannabinoids. Activation of mitochondrial CB1 receptors resulted in a decrease in PKA-dependent phosphorylation of oxidative phosphorylation proteins, including NDUFS2, and a concomitant decrease in brain mitochondrial function that was associated with cannabinoid-induced synaptic depression and amnesia (Hebert-Chatelain et al., 2016).

Sterol Carrier Protein 2 (Scp2) Is a Positional Candidate for Modulation of the Analgesic Response to THC

The SCP2 protein has been proposed to regulate brain endocannabinoid levels and, thus, endocannabinoid system function (Liedhegner et al., 2014; Martin et al., 2019). The Scp2 gene contains a putative splice region variant and demonstrates genetic modulation of expression (cis eQTL) in naive BXD neocortex and prefrontal cortex.

Ribosomal Protein S6 Kinase (Rps6kb1) Is a Positional Candidate for Modulation of the Effects of THC on Both Mobility and Hypothermia in Males

The gene Rps6kb1 encodes the 70-kDa ribosome protein S6 kinase (p70S6K) and contains a predicted splice region variant. Expression of Rps6kb1 is modulated by a cis eQTL in neocortex. Activation of the G-protein coupled CB1 receptor by THC triggers cascades of intracellular signaling, including activation of phosphoinositide-3 kinase (PI3K)/Akt/glycogen synthase kinase 3 (GSK-3) and subsequent activation of the serine/threonine kinase mammalian target of rapamycin (mTOR). P70S6K is a downstream target of mTOR whose activation is associated with protein synthesis and regulation of different cellular states (e.g., survival, growth, or autophagy). The CB1/mTOR/p70S6K signaling pathway in hippocampal GABAergic interneurons was found to mediate THC-induced long-term memory deficits (Puighermanal et al., 2009, 2013).

Phosphodiesterase 4D, cAMP Specific (Pde4d) Is a Positional Candidate for Modulation of the Hypothermic Response to THC in Females

Phosphodiesterase 4 hydrolyzes cAMP and plays a role in mediating cAMP signaling pathways and neuroadaptations favoring drug reinforcement, tolerance, and dependence (Cherry and Davis, 1999; Perez-Torres et al., 2000; Muschamp and Carlezon, 2013; Nestler, 2016). Phosphodiesterase inhibitors have shown some efficacy in reducing drug seeking behavior or intake of psychostimulants, alcohol, and opioids in preclinical models [reviewed in Olsen and Liu (2016)]. There are four isoforms of Phosphodiesterase 4, including the positional candidate PDE4D, which expresses 11 alternative splice variants and exhibits widespread expression in brain (Olsen and Liu, 2016). Expression of Pde4d is modulated by cis eQTLs in naïve BXD neocortex and striatum and the gene locus contains a missense variant as well as a putative splice region variant.

Phosphate and Tensin Homolog (Pten) Is a Positional Candidate for the Modulation of the Hypothermic Response to THC in Males

Highly expressed in brain, Pten encodes a dual protein and lipid phosphatase enzyme that has been proposed to interact with neurotransmitter receptors, including NMDA receptors (Ning et al., 2004) and serotonin 5-HT2C receptors (Ji et al., 2006), and regulate their activity. Disruption of the interaction between PTEN and serotonin 5-HT2C receptors in the ventral tegmental area (VTA) has been shown to inhibit the firing rate of dopaminergic VTA neurons projecting to the nucleus accumbens (Ji et al., 2006). Moreover, disruption of this interaction also blocked the ability of THC to enhance the firing rate of dopaminergic VTA neurons and subsequent THC-induced conditioned place preference (Ji et al., 2006) suggesting that interactions between PTEN and serotonin receptors may mediate some of the behavioral effects of cannabinoids and other drugs of abuse (Maillet et al., 2008). The expression of Pten in naive BXD nucleus accumbens is modulated by a cis eQTL and there are 59 variants located within or near the gene locus, although the impact of these variants on gene function or expression is not clear.

Positional Candidates Identified as Putative Modulators of Initial Responses to THC Have Been Previously Associated With Drugs of Abuse

For each high priority positional candidate or QTL gene (QTG) we explored broader associations with over 300 addiction and psychiatric disease key words using a newly developed resource: GeneCup (see section “Materials and Methods”). This search finds sentences in PubMed that contain both the gene of interest and keywords relevant to drugs of abuse, addiction, and psychiatric disease. Our overall goal was to place provisional endocannabinoid signaling QTGs identified in our study into a broader biological context. Given the broad modulatory effects of the endocannabinoid system on brain function, we hypothesized that THC initial response QTGs would also be associated with diverse signaling in response to other drugs of abuse and psychiatric disorders.

Pten (Chr 19 male Hypothermia QTG), Rps6kb1/p70S6K (Chr 11 Mobility and male Hypothermia QTG), and Pde4d (Chr 13 female Hypothermia QTG) were the most well connected QTGs in terms of associations with drugs of abuse and psychiatric disease (Supplementary Figure 3). All three of these QTGs are associated with terms related to stress, anxiety, depression, and schizophrenia. Rps6kb1/p70S6K and Pten were associated with terms related to alcohol, amphetamine, cocaine, nicotine, opioid, phychedelics, and sensitization. Pde4d was associated with terms related to nicotine and psychedelics. Pten, Rps6kb1/p70S6K, and Pde4d were also associated with neuroplasticity, neurotransmission, signaling and transcription terms in support of their known roles in these biological processes. Moreover, these three QTGs have also been implicated in brain regions relevant to drug abuse and psychiatric disease as evidenced by their literature co-citation with relevant terms (e.g., accumbens, amygdala, cortex, habenula, hippocampus, or striatum).

Although less well-connected, the QTGs Ndufs2 (Chr 1 antinociception QTL) and Scp2 (Chr 4 antinociception QTL) were also associated with neuroplasticity, signaling, transcription, stress, and hippocampus terms (Supplementary Figure 3). In addition, Ndufs2 was associated with schizophrenia terms and Scp2 was associated with alcohol-related terms.

Initial Responses to THC Are Correlated With Initial Responses to Other Drugs of Abuse and Cocaine Self-Administration

From our background list of 302 BXD legacy traits, we found 43 significant (p < 0.05) correlations (Supplementary Table 5) between 40 drug response or self-administration traits and one or more of our THC initial response traits. Several interesting correlation patterns were evident between initial responses to THC and BXD legacy cocaine, ethanol, and morphine response traits (Table 4). One of the most striking patterns was an enrichment of significant and positive correlations between THC-induced antinociception and initial locomotor responses to morphine (50 mg/kg i.p.) from 90 to 150 min post-injection. Other patterns emerged as well: (1) significant and positive correlations between THC-induced hypothermia (females) and naloxone-induced (30 mg/kg i.p.) withdrawal from morphine (50 mg/kg i.p.), (2) significant and positive correlations between cocaine self-administration and THC-induced immobility, (3) sex-specific and significant correlations between THC-induced hypothermia and preference for cocaine, and (4) significant correlations between both THC-induced immobility and hypothermia (females) and initial motor response and tolerance to alcohol/ethanol.

TABLE 4

Summary of significant (p < 0.05) correlations between initial response to THC and BXD legacy drug response traits.

The role of sex on trait covariation was not explicitly quantified in this exploratory analysis. However, strong correlations were detected between each THC trait and BXD legacy trait measurements both when sex was combined or considered separately. Finally, correlation patterns did not appear to be driven by genetic co-regulation of cannabinoid and drug response traits from shared loci (i.e., peak QTLs for correlated traits were non-overlapping). Instead, these correlation patterns indicate shared biological mechanisms and more complex genetic modulation driving response to THC and other drugs of abuse.

Discussion

In this study, we demonstrate that strain differences in initial responses to THC between B6, D2, and their BXD progeny are heritable (Figure 2 and Table 1). Thus, genetic mapping of THC response traits in the BXD family will lead to the identification of genomic loci, including variants and genes, that modulate response to cannabinoids. Candidate QTGs will almost certainly be direct or indirect effectors of endocannabinoid signaling given that the THC response traits measured in our study are entirely downstream of CB1 receptor signaling. Following this line of reasoning, we identified multiple QTLs and positional candidate endocannabinoid signaling QTGs (Pten, Rps6kb1/p70S6K, Pde4d, Ndufs2, Scp2) that may control variation in initial response to THC (Tables 2, 3). Many of these candidates mediate responses to other drugs of abuse (Supplementary Figure 3). Moreover, we found significant correlations between initial responses to THC and behavioral responses to cocaine, alcohol, and morphine (Table 4). Taken together, we provide strong evidence that gene variants in endocannabinoid signaling pathway genes are responsible for individual variation in initial THC responses and are likely to cause variation in the behavioral responses to other drugs of abuse.

Our genetic screen for natural variation in traits modulated by effectors of endocannabinoid signaling was designed to identify specific genes and pathways involved in the response to THC. Targeted genetic deletion has revealed insight into the integral role of CB1 receptor in this response. However, little is known about the role of natural variation in response to cannabinoids. A handful of genes (e.g., CADM2, NCAM1, NRG1, CSMD1, and CHRNA2) have been associated with cannabis use or dependence in human genome-wide association studies (Han et al., 2012; Sherva et al., 2016; Stringer et al., 2016; Pasman et al., 2018; Demontis et al., 2019). The underlying mechanisms whereby these genes mediate use or dependence to cannabinoids or interact with the endogenous cannabinoid system has yet to be elucidated, with one notable exception. Heterozygous Nrg1 (Neuregulin) mutant mice, in which one copy of the gene contains a deletion of the transmembrane region, are more sensitive to the initial locomotor suppressant and behavioral effects of THC (Boucher et al., 2007) and develop more rapid tolerance to the motor suppressant effects (Boucher et al., 2011) relative to control mice. In naïve heterozygous Nrg1 mutant mice relative to controls, these differences in cannabinoid responses were preceded by modest increases of CB1 receptor levels in substantia nigra and significant decreases in both thalamic NMDA receptor levels and striatal dopamine D2 receptor levels (Newell et al., 2013). The association of NRG1 variants with cannabis dependence in humans and independent preclinical evidence for a role of Nrg1 in initial response to cannabinoids is striking. Identification of additional genes that modulate endocannabinoid signaling and cannabinoid response is non-human genetic populations is clearly warranted and could help clarify associations determined in humans.

To date, no human studies have addressed the impact of genetic variation on acute behavioral responses to cannabis or derived cannabinoids. Greater understanding of the genes and pathways mediating responses to cannabinoids is important for at least two reasons. First, identification of variants in endocannabinoid signaling pathways will help us better understand biological responses to cannabinoids, including risk of cannabinoid dependence, and provides the opportunity to decouple unwanted side effects with therapeutic properties of cannabinoid drugs (e.g., motor impairment versus analgesic effects of THC). Second, endocannabinoid pathways represent common points of convergence that are highly relevant for understanding the molecular response to many drugs of abuse and molecular mechanisms underlying addiction processes. To this end, we have identified strong correlations between initial responses to THC and initial response or self-administration of cocaine, morphine, and/or alcohol among BXD strains that could be indicative of shared genetic and biological co-regulation.

The endocannabinoid system (e.g., lipid ligands, biosynthetic and catabolic enzymes, G-protein coupled receptors and effectors of signaling) plays an important role in modulating synaptic plasticity, behavior, and reward circuitry. Cannabis, cannabinoids and non-cannabinoid drugs of abuse have been shown to disrupt endocannabinoid system function and signaling, and these changes may contribute to addiction processes (Parsons and Hurd, 2015). For example, exposure to alcohol and opioids has been shown to increase the level of endocannabinoid system ligands (AEA or 2-AG) in rodent brain reward regions (Caille et al., 2007; Alvarez-Jaimes et al., 2009a,b; Ceccarini et al., 2013). Moreover, activation of CB1 receptors in some rodent models has been shown to enhance the rewarding effects of alcohol and opioids through both dopamine dependent and independent pathways (Ellgren et al., 2007; Serrano and Parsons, 2011; Panagis et al., 2014). Specifically, adolescent THC treatment enhanced adult heroin self-administration in Long-Evan male rats (Ellgren et al., 2007) and loss of CB1 receptors in mice or receptor antagonism in rats and mice reduces or attenuates alcohol preference and consumption while receptor activation enhances both behaviors (Serrano and Parsons, 2011). THC and other CB1 receptor agonists are often found to mediate reinstatement of drug-seeking behavior for alcohol or heroin in rats (Panagis et al., 2014). Both cocaine (Fourgeaud et al., 2004; Liu et al., 2005; Grueter et al., 2006; Pan et al., 2008) and alcohol (Clarke and Adermark, 2010; Adermark et al., 2011) exposure have been shown to alter CB1 receptor-mediated plasticity in rodent brain reward regions (e.g., striatum, ventral tegmental area, nucleus accumbens, and/or bed nucleus of the stria terminalis). CB1 receptor-mediated plasticity and signaling is also important for mediating response to opioids (Guegan et al., 2016; Zhao et al., 2017). Thus, the endocannabinoid system converges at many points with other neurotransmitter systems to mediate drug response and alter behavior following drug exposure.

We hypothesized that correlations between initial response to THC and responses to alcohol, cocaine, or opioids in the BXD population resulted from the convergence of genes and pathways that mediate both effects. Indeed, we observed striking correlations between THC and opioid responses (enhanced sensitivity to the analgesic effects of THC and enhanced sensitivity to the motor stimulant effects of morphine), THC and cocaine responses (enhanced sensitivity to the motor suppressant effects of THC and higher active lever presses for cocaine), THC and ethanol responses (enhanced sensitivity to the motor suppressant and hypothermic effects of THC and enhanced sensitivity to the motor stimulant effects of ethanol), and THC and withdrawal from alcohol and opioids (enhanced sensitivity to the motor suppressant effects of THC and higher levels of withdrawal from alcohol in contrast to lower levels of precipitated withdrawal from opioids). The results of these exploratory correlational analyses hint at underlying shared genetic and biological regulation. Of interest, no single shared locus could account for the covariance in drug response traits, indicating complex regulation by multiple loci, genes, and variants. Importantly, our exploratory analysis identified individuals of the BXD population with extreme responses to multiple drugs of abuse. These individuals can be leveraged to directly test predictions about drug response and to identify mediators of behavioral responses to both cannabinoid and non-cannabinoid drugs using both genetic (QTL mapping) and pharmacological (through targeting of specific pathways) methods.

Our study marks the first and largest attempt to quantify genetic factors mediating initial responses to THC in males and females of a genetic population. While we have made an important contribution to the field, there are some limitations of our study that will need to be addressed by future experiments. First, THC traits were profiled in a subset of 20 BXD strains which limits our ability to detect QTLs of small effect size. All QTLs identified in our study are suggestive after accounting for the effects of genome-wide testing. It is likely that these loci will be replicated in an independent or larger cohort of BXD strains. We emphasize that this replication will be required to confirm QTLs and QTGs prioritized by our study. Other aspects of our study design limit translatability somewhat and should be addressed by future studies. These include the use of a single major component of cannabis (THC), use of a single dose (10 mg/kg), and quantification of acute physiological responses to THC (as opposed to chronic exposure, dependence, and/or withdrawal). Finally, we were agnostic to genetic variation in cannabinoid metabolic pathways in this study.

Drug metabolism is an important biological process that impacts response to drugs of abuse and susceptibility to addiction. Based on trait correlation structure and genetic mapping, it is unlikely that strain differences in THC metabolism contributed to trait variation in all three of our initial response traits (e.g., mobility, hypothermia, and antinociception). It will be important in the future to establish whether sex or strain differences in THC metabolism contribute to variation in response and whether there are functional variants that impact THC metabolism in rodents. Several functional human variants (alleles) have been identified for cannabinoid metabolism genes [e.g., CYP3A4 and CYP2C9, reviewed in Hryhorowicz et al. (2018)]. The impact of these variants on response to cannabinoids and addiction is not clear, although a small study reported a trend toward increased sensitivity to THC in individuals homozygous for the CYP2C9∗3 allele associated with slower clearance of THC (Sachse-Seeboth et al., 2009).

Ultimately, our study demonstrates the feasibility of genetic mapping in the BXDs to identify underlying biological and genetic factors contributing to response to THC and other drugs of abuse. We have taken important first steps toward the identification of loci and genes that modulate initial responses to THC. We have also provided evidence of covariation among drug response traits, although the underlying causal factors remain elusive. Finally, we have identified BXD individuals with extreme responses to acute THC exposure. Identification of BXD individuals with greater initial sensitivity to THC will facilitate more detailed genetic analysis of intake, reward, withdrawal, tolerance, and the behavioral impact of cannabinoids on motivation, cognition, and health. The genetic architecture of behavioral and pharmacological responses to cannabis and cannabinoids is difficult to reconstruct using currently available human data sets. A better understanding of trait correlation structure and regulation by underlying genetic variation in diverse rodent models is expected to help bridge these gaps and direct the focus of future human studies.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by The University of Tennessee Health Science Center Institutional Animal Care and Use Committee.

Author Contributions

BJ, BM, and MM designed the study. BM, HC, and RW provided databases, tools, reagents, and animals. CP and MM performed the experiments and wrote the manuscript. CP, CR, PP, and MM analyzed the data. BJ, BM, RW, HC, and MM performed the major editing. The study was carried out in the laboratory of MM with intellectual guidance and support from BJ, BM, and RW. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer MD declared a past co-authorship with the authors HC, RW, and MM to the handling editor.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.