1. Genetics and Genomics
  2. Neuroscience

Acute opioid responses are modulated by dynamic interactions of Oprm1 and Fgf12

  1. Paige M Lemen
  2. Yanning Zuo
  3. Alexander S Hatoum
  4. Price E Dickson
  5. Guy Mittleman
  6. Arpana Agrawal
  7. Benjamin C Reiner
  8. Wade Berrettini
  9. David George Ashbrook
  10. Mustafa Hakan Gunturkun
  11. Xusheng Wang
  12. Megan K Mulligan
  13. Caleb J Brown
  14. Eric J Nestler
  15. Francesca Telese
  16. Robert W Williams  Is a corresponding author
  17. Hao Chen  Is a corresponding author
  1. University of Tennessee Health Science Center, United States
  2. University of California, San Diego, United States
  3. Washington University School of Medicine, United States
  4. Marshall University, United States
  5. Ball State University, United States
  6. University of Pennsylvania, United States
  7. Nash Family Department of Neuroscience and Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, United States
https://doi.org/10.7554/eLife.108845.3
Share this article
Cite this article
  1. Paige M Lemen
  2. Yanning Zuo
  3. Alexander S Hatoum
  4. Price E Dickson
  5. Guy Mittleman
  6. Arpana Agrawal
  7. Benjamin C Reiner
  8. Wade Berrettini
  9. David George Ashbrook
  10. Mustafa Hakan Gunturkun
  11. Xusheng Wang
  12. Megan K Mulligan
  13. Caleb J Brown
  14. Eric J Nestler
  15. Francesca Telese
  16. Robert W Williams
  17. Hao Chen
(2026)
Acute opioid responses are modulated by dynamic interactions of Oprm1 and Fgf12
eLife 14:RP108845.
https://doi.org/10.7554/eLife.108845.3
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

This study integrates large-scale behavioral, genetic, and molecular analyses in animal models to investigate morphine response. Utilizing high-quality, time-series Quantitative Trait Loci (QTL) mapping, the work provides compelling evidential support for novel, time-dependent genetic interactions (epistasis). A fundamental result of this rigorous analysis is the discovery of a novel Oprm1-Fgf12-MAPK signaling pathway, which offers new insights into the mechanisms of opioid sensitivity.

https://doi.org/10.7554/eLife.108845.3.sa0

Abstract

We generated time-series data for 105 morphine- and naloxone-related traits across ~700 BXD mice (64 diverse strains for both sexes) for 3 hr after a single morphine injection. Variations in responses were mapped using genome sequencing-based genotypes. The locomotor responses to morphine mapped to the µ opioid receptor gene (Oprm1) on chromosome (Chr) 10 with a peak linkage of 12.4 (–logp). The B allele inherited from C57BL/6J was associated with up to 60% higher activity. This effect climaxed at 75 min but was exhausted by 160 min. A second major modulator of locomotion emerged after approximately 100 min. This locus was located on Chr 16 with peak linkage of 10.6 in females and included one compelling candidate, fibroblast growth factor 12 (Fgf12). A strong and transient epistatic interaction existed between the Oprm1 and Fgf12 loci during a short time window (45–75 min). In heterogeneous stock rats, we demonstrated that Oprm1 and Fgf12 were co-expressed in one subtype of Drd1+ medium spiny neuron. A Bayesian network analysis supported an Oprm1-to-Fgf12 network that involves a MAP kinase cascade that modulates FGF12 phosphorylation and locomotor activation. OPRM1 and FGF12 networks in human genome-wide association study (GWAS) data highlight enrichment of signals associated with substance use disorder. This study represents the first demonstration of a time-dependent epistatic interaction modulating drug response in mammals and the first linkage of Fgf12 to opioid-induced behavior.

Introduction

Understanding the molecular genetics of opioids is crucial for treatment of opioid use disorder (OUD). Great progress has been made recently in genome-wide association studies (GWAS) of OUD. Genetic variants of several genes have reached genome-wide significance, such as KCNC1, KCNC2 (Gelernter et al., 2014), CNIH3 (Nelson et al., 2016), RGMA (Cheng et al., 2018), KDM4A (Sanchez-Roige et al., 2021), and OPRM1 (Zhou et al., 2020). Recently, Deak et al., 2022, have identified an additional set of 19 independent risk loci for OUD. Most of these variants explain a relatively small fraction of the heritability of OUD and require validation.

The most replicated gene among these studies is the opioid receptor gene OPRM1 (Deak et al., 2022; Gaddis et al., 2022; Zhou et al., 2020), which encodes the μ opioid receptor (MOR). Morphine, the prototypic MOR agonist, inhibits adenylyl cyclase and decreases the transmission of nociceptive information (Yam et al., 2018). Activation of MOR also increases dopaminergic neuronal transmission to the nucleus accumbens (NAc) by inhibiting GABAergic interneurons in the ventral tegmental area (VTA) (De Vries and Shippenberg, 2002). Co-expression of both MORs and dopamine D1 receptors (DRD1)—but not D2 receptors (DRD2)— is required for the initial locomotor response to morphine (Severino et al., 2020). MOR activation in striatal neurons is also a sufficient signal of opioid reward (Cui et al., 2014). However, the downstream molecular networks associated with OPRM1 activation and signaling remain incompletely understood.

The lack of large and well-phenotyped populations is a major impediment to further discovery of specific genes and mechanisms that modulate different phases of OUD (Berrettini, 2017). Compared to human GWAS, genetic mapping using model organisms requires a much smaller sample size. Phenotypes can be measured more accurately and studied mechanistically with better control over environmental variables. Large families of recombinant inbred strains (RI) are especially powerful for identifying genetic drivers responsible for trait variation, such as responses to opioids. These families are made by crossing two or more inbred strains, followed by intercrossing and inbreeding progeny for more than 20 generations (Williams et al., 2001). The BXD RI family of mice was created by crossing C57BL/6J (B6) with DBA/2J (D2) (Ashbrook et al., 2021; Taylor et al., 1973). This family has been used to study the effects of alcohol and many psychoactive drugs (Belknap et al., 1993; Crabbe, 1998; Parks et al., 2021; Belknap et al., 1995), along with numerous other molecular and complex traits.

The original data were described in Philip et al., 2010, who identified the Oprm1 locus. Here, we significantly extended these analyses by using whole-genome sequence-based genetic maps (Ashbrook et al., 2021) and the GEMMA linear mixed model (LMM) (Zhou and Stephens, 2012) for mapping, which improves mapping power and controls for population structure better than previously used methods. In addition to Oprm1, our analysis identified a new morphine response locus on chromosome 16 (Chr 16). We exploited complementary rodent datasets—both whole-brain proteomics and single-nucleus RNA-seq (snRNA-seq)—and identified a new candidate gene—fibroblast growth factor 12 (Fgf12). We found a strong and transient epistatic interaction between the Oprm1 and Fgf12 loci in a short time window after morphine injection. The snRNA-seq data highlights a single subtype of Drd1-positive cell type in which Oprm1 and Fgf12 are co-expressed. We then created a probabilistic causal model formalizing a molecular scheme in which gene variants modulate key signaling molecules to control locomotion after morphine injection. Lastly, we integrated rodent data with comparable human gene expression to provide a better translational context for understanding the molecular mechanisms and cellular cascades that may underlie initial variation in human OUD responses.

Results

QTL mapping of morphine-induced locomotor responses

A total of 64 fully inbred strains from the BXD family (Ashbrook et al., 2021) were used. Each strain was represented by an average of six males and six females (Philip et al., 2010). Locomotion was recorded for 180 min after an acute morphine injection, with movement data binned for each 15 min period. All data were generated by co-authors (PED, GM). Data for the 10th time bin (135–150 min) were lost. Variation in locomotor responses revealed a skewed distribution (Figure 1—figure supplement 1). We therefore quantile normalized data (Figure 1—figure supplement 2) for most analyses, although we note that linkage statistics are robust with respect to this transformation.

Quantitative trait locus (QTL) linkage was computed using a subset of ~7000 informative sequenced-based markers selected via linkage disequilibrium pruning and using the latest version of GEMMA v0.98.5 (Prins et al., 2021, https://github.com/genetics-statistics/GEMMA) as implemented in https://genenetwork.org/. The genome-wide –logp significance threshold (p<0.05), determined by 1000 permutations, is approximately 3.77 for the BXD family. A list of genome-wide significant loci for both morphine-induced locomotion and naloxone-induced withdrawal was summarized in Table 1, with more details provided in Supplementary file 1. We analyzed both sexes separately and jointly to identify potential sex-specific genetic modifiers of morphine response because significant sex by strain interaction was detected in our prior analysis (Philip et al., 2010).

Table 1
Genome-wide significant loci morphine-induced locomotion or naloxone-induced withdrawal.
Locus position (Mb, GRC38)
LocusSex†Peak –logp‡ChrProxPeakDistal1.5 –logp intervalAdd. effect§GN short phenotype descriptionTime period (min)GN BXD trait
Mor1aF4.42172.377.678.46.1–0.42Morphine-induced locomotion, distance traveled0–1511583
Nalx1aMF3.961124.0128.5141.817.9–0.43Naloxone-induced withdrawal, salivationNA11875
Nalx1bMF4.141163.3167.2168.24.8–0.44Naloxone-induced withdrawal, salivationNA11875
Mor1bF4.321173.4173.5176.22.8–0.38Morphine-induced locomotion, distance traveled150–16511584
Mor4aM4.06419.721.733.113.4–0.39Morphine-induced locomotion, distance traveled45–6011331
Mor5aF5.0159.89.813.43.50.52Morphine-induced locomotion, distance traveled105–12011581
Mor5bF3.86576.182.290.414.20.41Morphine-induced locomotion, distance traveled30–4511587
Nalx6aMF/F4.696112.1116.9122.810.7–0.42Naloxone-induced withdrawal, horizontal activity0–1511871
Mor7aM/F5.28796.096.897.51.5–0.44Morphine-induced locomotion, distance traveled0–1511326
Mor9aF3.90983.183.986.13.00.45Morphine-induced locomotion, distance traveled120–13511582
Mor10a*M/F10.53105.65.69.64.0–0.59Morphine-induced locomotion, distance traveled30–4511330
Nalx12aMF/F4.911216.627.628.511.90.46Naloxone-induced withdrawal, horizontal activity0–1511871
Mor12aF4.031269.670.874.75.2–0.42Morphine-induced locomotion, distance traveled0–1511583
Mor12bF5.711298.699.6102.53.9–0.44Morphine-induced locomotion, distance traveled30–4511587
Nalx14aM4.001449.152.161.212.1–0.41Naloxone-induced withdrawal, number of jumpsNA11336
Nalx14bMF4.3614118.4118.5119.41.0–0.25Naloxone-induced withdrawal, horizontal activity0–1511871
Mor16a*F/M10.561627.027.529.92.8–0.61Morphine-induced locomotion, distance traveled165–18011585
Nalx16bMF4.021657.559.762.24.60.42Naloxone-induced withdrawal, salivationNA11875
  1. *

    Significant for morphine and naloxone treatments; best candidate for Mora10a is Oprm1; best for Mor16a is Fgf12.

  2. †

    MF = male and female joint mapping. M/F=significant in both sexes, but –logp value applies to first entry.

  3. ‡

    –logp values given for quantile normalized traits but verified on other scales. Values are genome-wide p≤0.05 at ≥3.80.

  4. §

    Positive additive effects mean higher trait values linked to D genotypes. Negative values linked to high B values.

  5. Strong candidate gene for Mor7a is Tenm4.

Throughout the first 2 hr after injection of morphine (0–135 min bins), we detected a single and highly significant locus on Chr 10 (Figure 1) with peak linkage that extends from about 5.7–8.9 Mb in males and from 8.9 to 9.6 Mb in females, and a maximal –logp score of 10 and higher for males and 9.6 and higher for females (Figure 2) using quantile normalized locomotor data. This strong locus was reported by Philip et al., 2010. But we also now detect a strong second locus for locomotor activation on Chr 16 between 25 and 30 Mb, that only emerges 90 min after injection and reaches its temporal peak more than 2 hr after injection (Figure 3). Linkage is significant in both sexes but is several orders of magnitude stronger in females than males—linkages of 10.6 and 4.28, respectively.

Figure 1 with 2 supplements see all
Download asset Open asset
Time series of quantitative trait loci (QTLs) for morphine-induced locomotor response.

Locomotion data were quantile-normalized and mapped against whole genome sequencing (WGS)-based genotype using GEMMA in https://genenetwork.org/. QTLs for male (a) and female (b) are stacked up by increasing 15 min time intervals. The dotted lines represent the genome-wide significance level of ~3.77. The color shaded areas indicate consistent associations across time bins for the Chr 10 locus (red) and the Chr 16 locus (blue). Strong correlation in morphine-induced locomotor response between male and female BXD strains is shown in Pearson’s correlation scatter plots during the 45–60 min (c) and 165–180 min (d) time frames.

Figure 2
Download asset Open asset
Quantitative trait loci (QTLs) of morphine-induced locomotion between 45 and 60 min on Chr 10.

(a) QTL for males with a peak –logp of 10.06 (n=63 strains). (b) QTL for females with a peak –logp of 9.60 (n=64 strains). (c) Zoomed-in view of QTL in males. (d) Zoomed-in view of QTL in females. (e) Cis-eQTL for Oprm1 in the nucleus accumbens (NAc) of the BXDs, with a peak –logp of 10.01 (n=34 strains). (f) Cis-eQTL for Oprm1 in the hippocampus of the BXDs, with a peak –logp of 9.7 (n=67 strains) on Chr 10 at 5.6 Mb. (g) Oprm1 neighborhood in BXD family with SNP densities. (h) Haplotype map of the eQTL region. The ‘B’ of BXD is the mother, and ‘D’ is the father. GEMMA with LOCO was used for all association mapping.

Figure 3
Download asset Open asset
Quantitative trait loci (QTLs) of morphine-induced locomotion between 165 and 180 min on Chr 16.

(a) QTL in males with a peak –logp of 4.28 (n=63 strains). (b) QTL for females with a peak –logp of 10.56 (n=64 strains). (c) Zoomed-in view of the QTL in males. (d) Zoomed-in view of the QTL in females. (e) Cis-eQTL in the striatum of the BXDs, with a peak –logp of 4.43. (f) Cis-eQTL in the ventral tegmental area (VTA) of the BXDs, with a peak –logp of 4.06. (g) Histogram of the normalized expression of Fgf12 in striatum and zoomed-in view of the QTL region. (h) Histogram of the normalized expression of Fgf12 in VTA and zoomed-in view of the QTL region. GEMMA with LOCO was used for all association mapping.

Males and females display roughly similar patterns of linkage on both Chr 10 and Chr 16. This is consistent with the generally strong positive sex correlations in morphine-induced locomotion responses across all 15 min intervals. Correlation coefficients peak at about 0.84–0.88 from 15 to 75 min after injection—the phase during which morphine induces the steepest increase in locomotion—from about 29 to 71 m per 15 min bin (Figure 1c and d). This is also the interval that corresponds to the peak linkage in both sexes near Oprm1. Thereafter, correlations drop but are still generally above 0.70 through to 120 min. During the last hour, correlations are lower—ranging from 0.57 to 0.68. This drop in correlation corresponds to the period during which effects of the Fgf12 locus are strongest and, conversely, effects of the Oprm1 locus are weakest. The distance traveled in the final interval (165–180 min) is down to a level slightly lower than that in the first interval—about 25 m—and the final correlation between sexes is 0.68. This is the interval during which the Fgf12 locus has a remarkably high linkage of 9.7 and 5.3 in females and males, respectively. In contrast, at this late stage, Oprm1 has linkage below 2.5 in both sexes.

Candidate gene identification for morphine-induced locomotor response

The μ1 opioid receptor protein (MOR or OPRM1) is expressed in multiple brain regions and is involved in opioid-induced reward and locomotor response (Contet et al., 2004). In mice, the very large Oprm1 gene is on Chr 10 between 6.76 and 7.04 Mb, just proximal to the linkage peak (Figure 2c and d, GRCm38 assembly). In the BXD family, there are over 40 known SNPs, indels, and larger variants segregating in this gene locus, but to the best of our knowledge, none alter protein sequence. Given its location just proximal to the QTL peak, its large size, and the high level of genetic variation, Oprm1 is an uncontroversial candidate for differences in acute morphine locomotor activation (Table 2). But there is some additional support—expression of Oprm1 mRNA has been quantified in many brain regions across the BXD family, and variants in and around this gene clearly modulate its expression most strongly in NAc and hippocampus (cis-eQTLs with –logp values >7.0, see examples in Figure 2e and f). Variation in the expression levels of Oprm1 from the B and D haplotypes could contribute causally to the Chr 10 locus. Of functional importance, the D allele has roughly 50% higher expression in both NAc (Figure 2e) and hippocampus (Figure 2f), an interesting observation given that BXD strains that inherit the B allele tend to have a much stronger initial locomotor response.

Table 2
List of candidate genes.
Candidate geneChrStart (Mb)Stop (Mb)Male trait IDFemale trait IDPhenotype
Tenm479697.5BXD_11326BXD_11572Locomotion 0–15 min after morphine injection
Oprm1105.69.6BXD_11331BXD_11588Locomotion 45–60 min after morphine injection
Oprm1105.69.6BXD_11339BXD_11596Differences in locomotion before vs. after injection
Slc7a7/slc7a81449.161.2BXD_11336BXD_11850Number of jumps after naloxone-induced withdrawal
Fgf121627.029.9BXD_11328BXD_11585Locomotion 165–180 min after morphine injection
Fgf121627.029.9BXD_11339BXD_11596Differences in locomotion before vs. after injection

  1. Potential candidate genes in QTL regions that were further investigated in this study are provided. For details on any given trait, please visit the corresponding GN website by embedding the trait ID in the URL such as https://genenetwork.org/show_trait?trait_id=11331&dataset=BXDPublish (for trait BXD_11331, morphine response (50 mg/kg ip), locomotion (open field) from 45 to 60 min after injection in an activity chamber for males [cm]).

We restricted our analysis of candidate genes for the novel Chr 16 locus to a 1.5 –logp confidence interval (corresponding to a 95% confidence interval) extending from 27 to 29 Mb (Figure 1). This small region contains only five protein-coding genes and one gene model (Gm10823) in the following proximal-distal order: Ostn, Uts2d, Ccdc50, Gm10823, Fgf12, and Mb21d2 (Figure 4). Fgf12, previously known as Fhf1, is the largest protein-coding gene in this interval (~600 kb), with a promoter located close to Mb21d2 (Figures 3c, d, 4). This is also the strongest biological candidate gene (Table 2) and is already known to interact with the C-terminal region of three sodium channel proteins—SCN9A (Wildburger et al., 2015), SCN8A (Liu et al., 2003), and SCN2A (Wildburger et al., 2015). Fgf12 encodes a protein that is a member of the fibroblast growth factor (FGF) family. FGFs are involved in response to alcohol (Even-Chen et al., 2017) and morphine (Flores et al., 2010) in key brain reward regions.

Figure 4
Download asset Open asset
Chr 16 locus at 165–180 min after morphine injection.

(a) Zoomed-in view of the Chr 16 locus, from 26 to 29 Mb, showing individual SNPs (blue dots) and genes (purple horizontal lines) in the region. (b) A screenshot of GeneNetwork showing SNPs that inhabit this region. Almost all B vs. D SNPs (orange hash along x-axis) are restricted to two regions.

The majority of the Fgf12 gene is situated in a region that is almost identical-by-descent between the B and D parental haplotypes. This is highlighted in Figure 4 as the very low SNP density. There are 12 known non-coding variants between B and D haplotypes and no known coding variants; however, there are very significant cis-acting expression QTLs associated with Fgf12 expression differences in the striatum (Figure 3e and g, –logp of 4.4) and the VTA (Figure 3f and h, –logp of 4.1). The B allele is associated with expression that is as much as 60% higher than that of the D allele. We identified no coding variants for the other genes in this interval, and Fgf12 showed the most robust cis-eQTL evidence.

An epistatic interaction between Oprm1 and Fgf12 loci modulates locomotor response to morphine

We detect large additive effects for the locomotor responses to morphine near Oprm1 at early time points (0–90 min) and near to Fgf12 at late time points (+90 min). We tested for a possible epistatic interaction between these distinct loci and discovered a strong but transient interaction in both sexes (Figure 5). The epistasis is strongest in the 45–60 min period—precisely at the midpoint between the early and additive Oprm1 effect and the late and additive Fgf12 effect (Figure 5c–g, e.g. male trait BXD_11331). Those BXDs that inherit the B haplotype of Oprm1 (defined as rs29339674, 6.7 Mb), as well as the D haplotype of the Fgf12 locus (defined as rs4169220 at 31.6 Mb), have unexpectedly high levels of locomotion compared to all of the other two-locus combinations (Figure 5, Figure 5—figure supplement 1). The –logp value of the epistatic component is in the range of 2.5–4.8 (MF = trait BXD_11845, –logp = 4.2; F trait = BXD_11588, –logp = 4.8; M trait = BXD_11331, –logp = 2.5) and is significant at p<0.05. The full model that includes this interaction term, as well as both additive effects at Oprm1 and at Fgf12, has a remarkably high –logp value of 11.9–13.2 (MF = 13.2, F=11.9, M=12.0). This epistatic interaction is preserved as the time window moves forward. For example, maps of the 60–75 min interval (trait BXD_18846) have an epistatic effect with a –logp of 3.0; an additive effect near Oprm1 with a –logp of 8.2; and an additive effect at Fgf12 of merely 0.14. The full model in this interval still has a very high cumulative –logp of 11.4. By 75–90 min (trait BXD_18847), the interaction effect has dropped to a –logp of 2.2, and the model has a cumulative –logp of 9.40. In marked contrast, over the two intervals from 90 to 120 min, we do not detect significant epistasis. At best, the interaction –logp is only in the range of 1.5–2.0. In the 120–135 and 150–165 min intervals, the originally powerful Oprm1 additive effect has faded to 2.74 and 1.15, respectively. There is no detectable interaction with the Fgf12 locus. However, the purely additive effect at Fgf12 is now much stronger and is detectable for the first time as an independent additive effect with –logp values of 2.1 and 4.2, respectively. Finally, in the last interval—165–180 min (trait BXD_11585)—we detect no linkage to the Oprm1 locus, or any epistasis, but we again pick up a highly significant additive effect at the Fgf12 locus with a peak –logp, i.e., 4.9.

Figure 5 with 1 supplement see all
Download asset Open asset
Epistatic interaction among genotype combinations.

(a) Males and (b) females with different combinations of the ‘B’ (i.e. B/B) or ‘D’ (i.e. D/D) genotypes yield different distances traveled after morphine injection over 120 min. Distance traveled (m) between the two genotypes for each loci is shown for each time point in (c–f) males and (g–j) females. Error bars represent standard errors.

Cell-type-specific co-expression of Oprm1 and Fgf12 in the rat NAc core

To further examine whether Oprm1 and Fgf12 were co-expressed in the same cells of the NAc, we performed snRNA-seq using a droplet-based approach (10x Genomics). Nuclei were isolated from microdissected NAc cores obtained from two female heterogeneous stock (HS) rats (Carrette et al., 2021). The NAc was chosen due to its central role in opioid reward and the observed strain differences in morphine-induced locomotion (Velásquez et al., 2019). After filtering out low-quality nuclei and potential doublets, a total of 4495 high-quality nuclei transcriptomes remained with a median number of 3363 transcripts (unique molecular identifiers) and 1861 genes detected per nucleus. We then performed SCT transformation (Hafemeister and Satija, 2019), dimensional reduction, and clustering using Seurat (Stuart et al., 2019) and identified 17 cell-type clusters (Figure 6a).

Figure 6 with 2 supplements see all
Download asset Open asset
Oprm1 and Fgf12 are positively correlated in rat nucleus accumbens (NAc) suggested by single-nucleus RNA-seq (snRNA-seq).

(a) UMAP visualization of cell clusters from 4495 nuclei from rat NAc core. (b) Dot plot showing the expression level of cell-type marker genes in cell clusters. The shade of dots denotes normalized and scaled average expression, and the size of dots denotes the percentage of cells expressing the gene in each cell cluster. (c, d) Violin plots indicating the normalized and scaled expression level of Oprm1 (c) and Fgf12 (d) across cell clusters. (e–g) Scatter plots showing the correlation relationship between Oprm1 and Fgf12 in all cells (e), D1-MSN-3 (f), and D1-MSN-2 (g). Nuclei counts (n), Pearson’s correlation coefficients (r), and p-values (p) are labeled on each panel. UMAP, uniform manifold approximation and projection; D1-MSN, Drd1-expressing medium spiny neuron; D2-MSN, Drd2-expressing medium spiny neuron; GABA, GABAergic inhibitory neuron; Glut, glutamatergic excitatory neuron; Ol, oligodendrocyte; OPC, oligodendrocyte precursor cell; Ast, astrocyte; Mg, microglial cells.

We annotated the cell clusters based on the expression level of established marker genes of the major cell types: neurons (Snap25, Syt1), excitatory neurons (Slc17a7), inhibitory neurons (Gad1), dopaminergic neurons (Ppp1r1b, Foxp2, Blc11b), oligodendrocytes (Mbp), oligodendrocyte precursor cells (Pdgfra), astrocytes (Aqp4), and microglia (Arhgap15). In addition, we identified seven subtypes of dopaminergic-receptor neurons. We used expression of Drd1 and Tac1 to define D1-type medium spiny neurons (D1-MSNs). The expression of Drd2 and Penk was used to define D2-type medium spiny neurons (D2-MSNs) (Figure 6b, n=1768 MSNs total, see Figure 6—figure supplement 1f). MSNs represented the majority (91.8%) of neuronal cells, as expected by previous knowledge of the NAc cellular composition.

Oprm1 was uniquely expressed in the D1-MSN-3 subtype, while Fgf12 was expressed across all cell populations (Figure 6c and d). We then performed Pearson’s correlation analysis to test the idea that the epistatic interaction between Oprm1 and Fgf12 is mediated by a specific cell type. While a weak but significant positive correlation (r=0.08, p=1.8e-8) between the expression of Oprm1 and Fgf12 (Figure 6e) was detected when all cells were included, a much stronger correlation (r=0.35, p=6.1e-8, Figure 6f) was found in D1-MSN-3 cells. In contrast, D1-MSN-2 had a significant weak negative correlation (r = –0.12, p=0.02, Figure 6g). In conclusion, D1-MSN-3 is the only cell cluster expressing a high level of Oprm1 in rat NAc, in which the expression levels of Oprm1 and Fgf12 are significantly and positively correlated.

Constructing and testing Bayesian networks

We hypothesize that genetic variants in the Chr 10 and Chr 16 QTL regions are driving differential expression of Oprm1 and Fgf12 in the BXD family and thereby modulating morphine-induced locomotor responses. We used the literature to guide our selection of mediators acting between Oprm1 and Fgf12 (Bhushan et al., 2017; Buchsbaum et al., 2002). For example, FGF12 (Sochacka et al., 2020) and MAPK8IP2 are linked to a few sodium channel components—SCN2A (Nav1.2) and SCN8A (Nav1.6) (Schoorlemmer and Goldfarb, 2001; Seiffert et al., 2022)—as part of a membrane-associated scaffold concentrated at the axon hillock. The MAPK8IP2 scaffold protein is thought to modulate JUN amino-terminal kinase signaling and also interacts with and controls the activity of MAPK8/JNK1 and MAP2K7/MKK7 (O’Leary et al., 2016).

We tested Scn2a, Scn8a, Map3k11, Map3k12, and Mapk8ip2 as candidate mediators of the Oprm1 and Fgf12 loci and variation in locomotor activity. Expression levels of these transcripts in striatum, VTA, and NAc of BXD strains were passed from GeneNetwork into the Bayesian Network (BN) Webserver (Ziebarth and Cui, 2017) to define the most probable network structure. We constrained the nodes into four tiers: Tier 1 contains the two loci as instrumental variables; tier 2 contains expression estimates of the two prime candidate genes, Oprm1 and Fgf12; tier 3 contains the sodium channels and MAP kinases that potentially mediate both additive and epistatic effects on locomotion; and tier 4 contains the locomotor outcomes expressed in four time-series bins.

The most strongly supported model (Figure 7) recapitulated the cis-eQTL results by indicating that the Chr 10 and Chr 16 loci control the expression of Oprm1 and Fgf12 mRNA. This model also supports the causal role of Oprm1 and its downstream MAP kinases on locomotor response during the first 120 min after injection. In contrast, the last 45 min are dominated by a direct Fgf12 effect. Further, Oprm1 and Fgf12 directly modulate Map3k11, which in turn modulates Mapk8ip2 and Map3k12 mRNA expression and regulates locomotion from 0 to 75 min. Mapk8ip2 is embedded in this network with multiple connections.

Figure 7 with 1 supplement see all
Download asset Open asset
Modeling the mechanistic interactions between genetic variants and phenotypes using Bayesian network.

This network illustrates the relationship between genetic variants in the Chr 10 and Chr 16 quantitative trait locus (QTL) regions and the differential expression of candidate genes Oprm1 and Fgf12 in the BXD family. Additional gene expression data of MAP kinases were retrieved from https://genenetwork.org/. Morphine-induced locomotion responses at different time bins are also included in the network. This causal hypothesis was developed using the Bayesian network framework available in GeneNetwork, where the arrow is indicative of the direction of the causal relationship. Additional genes included in the input of the network construction but had no connection to the final network were excluded in the illustration.

This BN aligns well with whole-brain proteomics data (Figure 7—figure supplement 1). In agreement with the mRNA data, Oprm1 and Fgf12 are expressed with opposite polarities, in which the B allele is associated with high expression of Oprm1 and low expression of Fgf12 and vice versa for the D allele. The three MAP kinases are also expressed in higher levels in animals that inherit B alleles than D alleles.

QTL mapping of naloxone-induced withdrawal responses

Once the locomotion test was completed (180 min after the morphine injection), naloxone was injected (30 mg/kg in isotonic saline at a volume of 10 ml/kg) to induce a morphine withdrawal response (Philip et al., 2010). Locomotion and other behavioral measurements were taken for both males and females. QTL analysis of these traits were not included in Philip et al., 2010. There are 16 genome-wide significant QTLs among 5 naloxone-induced withdrawal behavioral responses (Figure 8, Supplementary file 1).

Figure 8 with 1 supplement see all
Download asset Open asset
Quantitative trait loci (QTLs) for naloxone-induced morphine withdrawal responses in male and female BXD mice.

Behavior data were quantile-normalized and mapped against whole-genome sequencing (WGS)-based genotypes using GEMMA in https://genenetwork.org/. (a) Number of jumps 15 min after naloxone injection in males. (b) Number of jumps 15 min after naloxone injection in females. (c) Change in locomotion measured by the last 15 min of morphine locomotion response minus the first 15 min after naloxone injection in males, with a peak on Chr 10 and a peak on Chr 16. (d) Change in locomotion measured by the last 15 min of morphine locomotion response minus the first 15 min after naloxone injection in females, with a peak on Chr 16. (e) Horizontal activity (distance traveled) 0–15 min after naloxone injection for males. (f) Horizontal activity (distance traveled) 0–15 min after naloxone injection for females, with a peak on Chr 6 and a peak on Chr 12. (g) Number of beam breaks in an open field 0–15 min after naloxone injection in females, with a peak on Chr 6 and a peak on Chr 12. (h) Salivation level in females. Dashed lines: genome-wide significance threshold. Dotted lines: threshold for suggestive significance.

Figure 8c, d displays the differences in locomotion before and after naloxone injection, with a QTL peak on Chrs 10 and 16 that overlap the morphine locomotor QTLs in Figures 2 and 3 and that are likely to correspond to variants in or near to Oprm1 and Fgf12. Our results show that the genotype effect on locomotion was similar before and after naloxone injection (Figure 8—figure supplement 1). We evaluated the biological function of naloxone candidate genes using GeneCup (Gunturkun et al., 2022). Candidates for the locus on Chr 14 for naloxone-induced jumps include Slc7a7 and Slc7a8 (Table 2). Both are transmembrane amino acid transporter proteins (Hu et al., 2020) that dimerize with SLC3A2 (Hurkmans et al., 2022; Torrents et al., 1998; Torrents et al., 1999). Of relevance, SLC7A8 (LAT2) is an L-DOPA transporter (Crocco et al., 2024).

Integrating murine data with human GWAS results

There is strong support for OPRM1 in human OUD GWAS. OPRM1 variant rs1799971 has been associated with OUD in several studies (e.g. p=8 × 10–10 in Zhou et al., 2020, p=2 × 10–8 in Hatoum et al., 2022a, and p=4.92 × 10–9, OR = 1.046, in Deak et al., 2022). A human GWAS of opioid cessation reported linkage to the FGF signaling pathway in which the effect of FGF12 was nominally significant (p=0.0015, odds ratio = 1.24) (Cox et al., 2020). In gene-based analyses of OUD by Deak et al., 2022, there was a nominally significant signal for FGF12 (p=0.006), with the top FGF12 SNP rs1553460 also reaching nominal significance (OR = 1.015, p=0.021).

Considering the epistatic interaction between OPRM1 and FGF12 and the complexity of human genetic signals (many genes of small effects), it is possible that these two genes function in a network that contains many other genes. We thus developed a list of 500 genes that overlap with FGF12 and OPRM1 in the GTEx sample (GTEx Consortium, 2013) (GTEx v8). We then examined whether there is enrichment of members of this network using MAGMA (de Leeuw et al., 2015) in the GWAS by Zhou et al., 2020. Our initial analysis was performed for NAc, putamen, caudate, and adipose tissue for the FGF12 lists (the latter as a negative control).

For OUD, no significant enrichment was observed for the original 500-gene lists, likely due to the lists’ large size and the limited statistical power of the SUD GWAS. In the GWAS conducted by Deak et al., 2022, where FGF12 was nominally significant, we found nominally significant enrichment in two of the brain regions, putamen and caudate. Given the low power of this dataset and the broad initial gene list, we refined our analysis to focus on genes specifically correlated with OPRM1 and FGF12 in relevant tissues (see Methods). No enrichment was observed for the gene network in the OUD GWAS (p=0.493). However, this refined approach revealed nominal enrichment (p=0.038) for the addiction risk factor (Hatoum et al., 2022b).

Discussion

We analyzed time-dependent behavioral responses to morphine and naloxone collected from the BXD family of mice (Philip et al., 2010) using WGS-based genetic markers and LMMs. We discovered a novel association on a Chr 16 locus that overlaps Fgf12 in both sexes, in addition to confirming an association between locomotor response and a region on Chr 10 that overlaps Oprm1. Further, these two loci had a significant but transient epistatic interaction between 45 and 90 min after morphine injection. According to our transcriptomic data from NAc in rats, Oprm1 and Fgf12 are colocalized in a specific subtype of D1-MSN, and their expression levels are positively correlated. Analysis of both OPRM1 and FGF12 in human GWAS data demonstrated an enrichment of signals associated with SUD phenotypes and a modest corroboration of variants in the FGF12 locus on Chr 3q28.

Additive effect of the Oprm1 locus

The proximal Chr10 locus associated with early phase locomotor response to morphine contains the Oprm1 gene. This locus has been associated with the antinociceptive effect of morphine (Bergeson et al., 2001). Oprm1 encodes the primary receptor responsible for morphine-induced locomotor response (Severino et al., 2020; Smith et al., 2009) and is also responsible for its addiction liability (Ballester et al., 2022; Zhang et al., 2020). Notably, the Oprm1 mRNA cis-eQTL in the NAc and hippocampus precisely overlaps with the Mor10a QTL interval (Figure 2e and f), providing strong secondary evidence for Oprm1 as the causal gene. Morphine-induced locomotion was affected by Oprm1 variants (Popova et al., 2019). For example, in the Oprm1 A112G knock-in mouse, morphine-induced hyperactivity was blunted (Mague et al., 2009). Further, the effect of MOR on locomotion depends on the neuronal cell types involved. For example, when MOR was deleted from D1 neurons, mice had hypolocomotion in response to opioids. By contrast, when MOR was deleted from A2a (adenosine A2a receptor-expressing) neurons, mice displayed increased movement in response to opioids (Severino et al., 2020). These data indicate that Oprm1 is a strong candidate gene for the Chr 10 locus associated with morphine-induced locomotion response.

BXDs with the D allele at the Oprm1 locus were less sensitive to a 50 mg/kg dose of morphine (Figure 6). The locomotion data demonstrate a blunted motor activity response similar to what was seen in low doses of morphine (5–10 mg/kg) (Mori et al., 2004), while BXDs with the B allele demonstrate higher sensitivity at the same dose; strains with the B allele had an initial biphasic increase in motor activity followed by a drop. This suggests that the B allele of Oprm1 is associated with greater sensitivity or locomotion activity. These data match the differences between the parents of the BXD family. The locomotor response to morphine in C57BL/6J is more pronounced and lasts longer compared to that of DBA/2J (Murphy et al., 2001). This dimorphism may be linked to differences in morphine-induced dopamine release in the NAc in C57BL/6J than in DBA/2J mice (Murphy et al., 2001).

Additive effect of the Fgf12 locus

The association between the Chr 16 locus (Figure 1a and b) and morphine-induced locomotion 150–180 min after injection is a novel finding, as is the detection of a time-delimited epistatic interaction of the Chr 16 locus with the Oprm1 locus. The use of a high-density WGS-based marker set and the LMM (GEMMA) allowed us to detect this and other novel loci not previously detected. Of all positional candidate genes in the Chr 16 locus (Figure 4), Fgf12 was the most biologically relevant and was supported by strong cis-eQTL in striatum and VTA (Figure 3e and f). Studies have implicated members of the FGF family, such as Fgf2 (Even-Chen and Barak, 2019a), FGF receptor 1 (Even-Chen and Barak, 2019b), and FGF receptor 2 (Blackwood et al., 2019) in substance use. It has recently been demonstrated that administration of an Fgf21 analog in non-human primates decreased alcohol consumption by 50% via an amygdala-striatal circuit (Flippo et al., 2022). Unlike these secreted FGFs, Fgf12 is an intracellular protein that serves as a cofactor for voltage-gated sodium channels and other molecules and is involved in intracellular signaling events (Schoorlemmer and Goldfarb, 2001). While the full spectrum of its function remains unknown, Fgf12 has demonstrated a role in locomotion. For example, mice with a null mutation of Fgf12 have ataxia (Goldfarb et al., 2007). Elevated stress responses have been associated with a reduction in Fgf12 expression in the prefrontal cortex in rats (McCreary et al., 2016). FGF12 expression is also elevated in the anterior cingulate cortex of patients with major depressive disorder (Evans et al., 2004). Variants in human FGF12 have been linked to cognitive decline (Mohammadnejad et al., 2020), major depression (Uher et al., 2013), schizophrenia (Stertz et al., 2021), and sleep quality (Byrne et al., 2013). Recent evidence also supports the druggability of other intracellular FGFs, such as FGF14 (Dvorak et al., 2025; Singh et al., 2025) and FGF13 (Singh et al., 2025), through their interactions with sodium channels. Our findings suggest that FGF12 could represent a novel therapeutic target for OUD by modulating its interaction with sodium channels.

The epistasis of morphine locomotor activation

We identified a highly significant epistatic interaction between the Oprm1 and Fgf12 loci 45–90 min after morphine injection (Figure 6), which erodes in the following 30 min, providing strong pharmacological constraints on mechanisms. Those strains that inherit both the B/B genotype at Oprm1 and the D/D genotype at Fgf12 have significantly higher activity levels compared to those strains that inherit the other three genotype combinations. Our mapping extends in intervals of 15 min over a 3 hr period. We can divide this interval into four phases (Figure 1). The first phase after the morphine injection is influenced only by DNA variants in the Oprm1 locus. The second phase—from 45 to 90 min—is characterized by the strong but transient epistatic interaction of Oprm1 and Fgf12 loci. The third phase is a kind of hiatus that extends from 90 to 120 min in which we do not detect any epistatic interactions. Only Oprm1 is active during this phase. In the final phase—from 120 to 180 min—the Oprm1 effect is exhausted and now the Fgf12 locus acts independently for the first time with a purely additive effect on activity. The hiatus phase may be explained mechanistically by effects of other loci that we are not yet able to detect reliably. Note, for example, the transient female-specific locus on Chr 5 that fills part of the hiatus. This complex time course of transitioning in genetic control on morphine-induced locomotion is likely caused by molecular cascades triggered by the activation of the OPRM1 receptor and its signaling molecules and downstream molecular network linked to FGF12. To some extent, we have tried to bridge the gap between the genetics of morphine responses to possible molecular cascades using Bayesian causal modeling as in Figure 8.

Complex traits, such as OUD, are controlled by many genes that interact with each other and their environment (Crist et al., 2019). While improved statistical approaches (D’Silva et al., 2022; Wei et al., 2014) and novel machine learning methods (Chicco and Faultless, 2021) are starting to enable the study of epistasis in human data, most human genetic studies either did not have the power to detect epistasis or ignored its effect (Carlborg and Haley, 2004; Uffelmann et al., 2021). In contrast, many epistatic interactions have been identified using model organisms (Mackay, 2014), mostly due to the ability to obtain data from individuals with controlled genotypes and measuring phenotypes in well-controlled environments (Campbell et al., 2018). Identifying gene-gene interactions in model organisms can provide candidates to be evaluated in the human population for its potential in improving the prediction of individual disease risk and the application of personalized therapies.

Cell-type-specific gene expression in NAc

While gene interaction is nearly always measured statistically without regard to biological mechanism, we explored the utility of snRNA-seq in identifying cellular and molecular mechanisms of this interaction by examining gene expression in NAc, a brain region relevant to the biological effect of morphine. We focused on NAc because differences in morphine-induced dopamine release in NAc correspond with the locomotor response in the parental strains of the BXD mice (Murphy et al., 2001). Further, mice lacking Drd1 are not responsive to acute cocaine- (Karlsson et al., 2008; Xu et al., 1994) or morphine- (Wang et al., 2015) induced locomotor response. While our analysis identified 17 distinctive cell types (Figure 7a), Oprm1 was detected almost exclusively in the D1-MSN-3 subtype that has high expression of Foxp2 and low expression of Ppp1r1b (i.e. Darpp-32; Figure 7b and c). The total number of genes (n_feature RNA) and other quality control (QC) parameters in D1-MSN-3 are comparable with other neuronal cell types. In contrast, Fgf12 was expressed in most of the neuronal and glial cell types (Figure 7d). The expression levels of Oprm1 and Fgf12 were positively correlated (r=0.36, p=6e-8) in D1-MSN-3 but not in D1-MSN-2 cells, indicating a potential cell-type-specific mechanism for the epistatic interaction.

MSNs in the NAc have been the subject of intense research. A recent snRNA-seq study of the NAc in rat brains also detected this specific subtype expressing high levels of Oprm1, which is selectively labeled by the expression of Chst9 (Andraka et al., 2023; Savell et al., 2020). These data suggest that D1-MSN-3 is a highly likely location for the interaction between Orpm1 and Fgf12. Notably, this cell type is conserved also in primate and human NAc snRNA-seq data (Andraka et al., 2023), suggesting that it may play a critical role in opioid responses across species.

To further confirm our results and allow for analysis of correlations, we queried data from the Ratlas (https://day-lab.shinyapps.io/ratlas/), which reveals that Fgf12 is indeed widely expressed in all cell types in other datasets, and that Oprm1 was not detected in either Drd1 or Drd2 MSN neurons. This is most likely due to the lack of detection of the MSN subclass that is equivalent to the D1-MSN-3 subtype that we identified (Figure 6, Figure 6—figure supplement 1). However, cultured cell data do express Oprm1 in both Drd1 and Drd2 cells (Figure 6—figure supplement 2). While out of the scope of this paper, future work will be needed to assess the impact of the Oprm1-Fgf12 interaction on the activity of the D1-MSN-3 population in mediating opioid responses.

Bayesian causal network modeling of OPRM1 and FGF12

We included gene expression data of several MAP kinases and sodium channels with known relationships with Oprm1 and Fgf12 as the input for our BN. The resulting network (Figure 8) indicated that Mapk8ip2, Mapk3k11 (Mlk3), and Map3k12 (Dlk) form a highly interactive network mediating the interaction between Oprm1 and Fgf12 in regulating their effects on morphine-induced locomotor response. In addition to validating the results on genetic influence on Oprm1 and Fgf12 expression and morphine-induced locomotor response, the network identified Map3k11 as a potential nexus that links Oprm1 and Fgf12. Map3k11 further activates other MAP kinases such as Mapk8ip2 and Map3k12 to modulate locomotion before 120 min. Thus, Map3k11 is a potential mediator of the epistasis between Oprm1 and Fgf12. Our BN did not contain a direct relationship between Fgf12 and Mapk8ip2, the encoded proteins for which directly bind to each other (Bhushan et al., 2017). It is possible that these two proteins do not regulate each other’s mRNA levels, which was used in our BN. While our modeling suggests that the epistatic interaction between Oprm1 and Fgf12 is mediated through Map3k11, further research, such as those using Mendelian randomization studies and other forms of mediation analysis, could provide further understanding of the interactions between MAP kinases, FGF12, and OPRM1.

Genetic modulation of naloxone effects

Naloxone is a potent opioid antagonist that reverses opioid overdoses (Theriot et al., 2022) and is often used in treatment settings or as a primary overdose prevention method. Pharmacogenomics has primarily focused on individualized treatment for pain, with less priority on OUD. However, understanding the role genetics plays in naloxone’s mechanism could initiate potentially more personalized and effective treatments for OUD. We identified multiple significant associations between behavioral measures of the effect of naloxone on Chrs 3, 5, 7, 10, 11, 12, 13, 14, and 18. These loci are likely to contain novel genes that play a role in OUD. For example, Slc7a7 and Slc7a8 are candidate genes for the number of jumps after naloxone injection. Both of these genes are a part of the solute carrier family, which play a large role in the absorption of drugs and xenobiotics (Lin et al., 2015), and are expressed in the brain. Similar to the morphine locomotion time-series phenotypes, there are significant QTL for the differences in locomotion before and after naloxone injection (Figure 8c and d) on Chrs 10 and 16. Mice with the B allele at the Fgf12 locus still had elevated locomotion during 165–180 min compared to those with the D allele. Injection of naloxone eliminated the remaining effect of morphine and brought locomotion in both genotypes close to baseline level (Figure 8—figure supplement 1). The difference in locomotion before vs. after naloxone injection thus amplified the net effect of morphine at this late phase and further confirmed that this effect was mediated via the MOR. The highly significant association (–logp = 5.34 in males and –logp = 9.57 in females) between the Chr 16 locus and this phenotype provided a strong confirmation for its role in late phase response to morphine. Analysis of these genes in the homologous human regions in larger populations are needed to translationally confirm and extend the validity of these results.

Integration with human data

While our attempt to integrate murine data with human GWAS results demonstrated that the pipeline works, the lack of enrichment for the network in OUD GWAS could likely be due to the small sample size in the original human data for that specific SUD. This could also likely be due to the lack of data specification, such as the differences between examining physical dependence vs. the development of OUD, initial use vs. problematic use, and withdrawal vs. relapse (Sanchez-Roige et al., 2021).

In summary, our study confirmed that the Oprm1 locus strongly modulates the early phase locomotor response to morphine and identified a novel association between Fgf12 locus and the late phase response to morphine. To the best of our knowledge, the epistatic interaction between the Oprm1 and Fgf12 loci during the middle phase of locomotor response is the first demonstration of a transient time-dependent epistatic interaction modulating drug response in mammals—a finding with interesting mechanistic implications. Our snRNA-seq analysis suggests that the D1-MSN-3 subpopulation of dopaminergic-receptor neurons might be critical in the epistatic interaction between Oprm1 and Fgf12. Our study further shows that the interaction between the Fgf12 and Oprm1 genes is mediated by Mapk8ip2, Map3k11, and Map3k12, and that the activity of these proteins is modulated by the presence of Fgf12 variants. Finally, this work demonstrates how high-quality Findabile, Accessible, Interoperabile, and Reusable (FAIR+) phenotypes can be used with updated datasets to yield striking results, and how joint mouse and human neurogenomic and GWAS results can be merged at gene and network levels for bidirectional validation of SUD variants and molecular networks.

Methods

Animal and behavior data collection

Morphine-induced locomotor responses and naloxone-induced withdrawal were recorded for 64 BXD RI strains. All animals were young adults (8–9 weeks) reared at the Oak Ridge National Laboratory in 2007–2008. Detailed methods were reported in the original publication by Philip et al., 2010. Briefly, the testing protocols included giving each mouse a single i.p. injection of morphine sulfate (50 mg/kg in isotonic saline at a volume of 10 ml/kg), followed by immediately placing the mice into an activity chamber (43.2 cm L×43.2 cm W × 30.4 cm H, ENV-515, Med Associates, St Albans, VT, USA). Each chamber contained two sets of 16 photocells placed at 2.5 and 5 cm above the chamber floor. Activity was measured as photocell beam breaks and converted into horizontal distance traveled (cm). This behavior and rearing were recorded for 3 hr and were exported as the sum of 15 min bins. Then, each mouse in these cases received an injection of naloxone (30 mg/kg in isotonic saline at a volume of 10 ml/kg, i.p.) and were immediately returned to the chambers for an additional 15 min. Naloxone’s effects on locomotor activity levels and signs of withdrawal were recorded for 15 min, including the number of jumps, fecal boli, urine puddles, wet dog shakes, instances of abdominal contraction, salivation, ptosis, and abnormal posture. The number of animals used per strain was 7.8 (mean ± SD) for females and 6.6 for males. Testing occurred at about 8–9 weeks after birth.

QTL mapping in GeneNetwork

We reanalyzed 105 phenotypes from the Philip et al., 2010, study available in https://genenetwork.org/, a database containing phenotypes and genotypes, and also serves as an analysis engine for QTL mapping, genetic correlations, and phenome-wide association studies (Mulligan et al., 2017; Sloan et al., 2016; Watson and Ashbrook, 2020). The trait IDs and their most significant loci are shown in Supplementary file 2. The data for mapping QTLs consist of polymorphic genetic markers and quantitative trait values for strain means. Statistical heritability was also evaluated to estimate the degree of variation among the morphine and naloxone traits due to genetic variation. The original analysis (Philip et al., 2010) was performed with the expanded BXD strain to evaluate complementary behavioral phenotyping using Haley-Knott regression in GeneNetwork. In our reanalysis, we first quantile-normalized the trait data. We then used the newly implemented GEMMA (Zhou and Stephens, 2012) with the LOCO option, which corrects family structure, to perform genetic mapping using approximately 7000 markers obtained from whole-genome sequencing. Candidate genes were selected from the confidence interval of a one-LOD drop-off from the peak statistical significance as determined previously (Philip et al., 2010). The criterion for genome-wide significance was a –logp value of 3.77 for the BXDs. Loci labeled by genetic markers could act independently at all time points or possibly in an epistatic interaction at a few time points. Therefore, we also examined epistatic interactions between two loci by performing a pair-scan implemented in GeneNetwork (v1). This scan generated a matrix map output. Functional enrichment of these networks was obtained using WebGestalt (Zhang et al., 2005). The biological relevance of candidate genes in the context of SUD was queried using GeneCup (Gunturkun et al., 2022). We conducted most early exploratory analysis using the web version of GeneNetwork. For final analysis, we used the API interface of GeneNetwork (Mulligan et al., 2017) to conduct standardized analysis and generate uniform figures for all phenotypes of interest.

Brain samples for snRNA-seq

Brain samples from two female HS rats were obtained from the oxycodone tissue repository at UCSD (Carrette et al., 2021). These rats were selected based on their low addiction-linked behavior phenotypes from a larger cohort characterized in the oxycodone self-administration paradigm. The brain tissues were collected from rats euthanized after 4 weeks of prolonged abstinence from oxycodone intravenous self-administration (Carrette et al., 2021). Brain tissue was extracted and snap-frozen (at −30°C). Cryosections of ~500 μm (Bregma 2.28–0.72 mm) were used to dissect the NAc core punches on a −20°C frozen stage. Punches from three sections were combined for each rat.

snRNA-seq library preparations

snRNA-seq library was performed using the Chromium Next GEM Single Cell Multiome Reagent Kit A (catalog number 1000282) following Chromium Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Kits User Guide (10x Genomics). Approximately 10,000 nuclei were loaded per reaction, targeting recovery of 6000 nuclei after encapsulation. After the transposition reaction, nuclei were encapsulated and barcoded. Next-generation sequencing libraries were constructed following the User Guide. Final library concentration was assessed using Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific), and post-library QC was performed using Tapestation High Sensitivity D1000 (Agilent) to ensure that fragment sizes were distributed as expected. Final libraries were sequenced using the NovaSeq6000 (Illumina).

Bioinformatic analysis of snRNA-seq data

Raw sequencing data were converted to FASTQ format using bcl2fastq (Illumina). The FASTQ data were first aligned to the Rattus norvegicus mRatBN7.2/rn7 genome and then aggregated into one sample file using CellRanger version 2.0.0 using default parameters. The output files were analyzed using Seurat version 4.1.1 (Stuart et al., 2019). Nuclei with gene numbers between 500 and 6000, RNA counts between 1000 and 16,000, the percentage of mitochondrial gene reads lower than 2%, and the percentage of small 40S or large 60S ribosomal (Rps and Rpl) gene reads lower than 1% were considered high-quality cells and kept for further analyses. To eliminate potential doublets, we used DoubletFinder 2.0.3 and removed 5% nuclei identified as doublets (McGinnis et al., 2019).

High-quality singlet nuclei were then normalized and scaled using SCT transformation (Hafemeister and Satija, 2019) with the percentage of mitochondrial genes, RNA counts, and sample ID as covariates. The dimensional reduction was performed using PCA, and the first 30 PCs were used for KNN graph construction and clustering using the Louvain algorithm. Uniform manifold approximation and projection (UMAP) was used for visualization of the clusters.

The dot plots and violin plots were generated using Seurat (Stuart et al., 2019). Pearson’s correlation was computed in R 4.2.1 using cor.test, and p<0.05 was considered significant.

Human translational data

We reviewed OUD and other relevant published GWAS data in previous literature with variants in OPRM1 and FGF12. We also extended our review to looking at potential roles of FGF12 in humans specifically. Then, all the files for FGF12 and OPRM1 correlations for each individual brain region, as well as some controls, were downloaded using GTEx v8 (Battle et al., 2017) and GeneNetwork and then extracted out the shared genes with a liberal cutoff of 500 genes. We ran a genome-wide enrichment analysis using MAGMA, a powerful tool that identifies genes and gene-sets associated with a specific disease for analysis (de Leeuw et al., 2015). However, upon running our genome-wide enrichment analysis for OUD, we realized our initial search was too broad and much larger than in previous literature (see Results section). To reduce the amount of noise, we identified genes that overlap between FGF12 and OPRM1 in a given tissue and identified genes correlated specifically to each OPRM1 and FGF12 in several tissues. To do this, we developed a list of genes (and Ensembl IDs) in each brain region to examine the number of occurrences each gene is seen across this data. We narrowed down our search to the top genes with high counts and compared them to the same number of genes with a count of only one. The count is the number of times the gene is associated with both FGF12 and OPRM1 in all the regions/tissue datasets in GeneNetwork. Association is defined as among the top 500 genes ranked by correlation with FGF12 or OPRM1.

BN modeling

The BN server located at https://bnw.genenetwork.org/sourcecodes/home.php was used. Raw data files containing genotypes for the Oprm1 and Fgf12 loci, morphine-induced locomotion, and expression levels of relevant genes (sodium channels, MAP kinases) were transferred from GeneNetwork. We separated our nodes into four tiers to label chromosome location, candidate gene, signaling molecules, and phenotype (morphine locomotion), each as a separate tier. These are used to organize the nodes into different layers. Each tier represents a different layer of complexity. We permitted causal connections to flow from the first tier to the second tier, as well as from the second tier to the third tier. We also allowed interaction between nodes located within the second or the third tiers.

Data availability

All behavioral data are available from https://genenetwork.org/ (see Supplementary file 2 for trait IDs). All snRNA-seq data have been deposited into NCBI GEO (accession number: GSE214388).

The following data sets were generated
    1. Zuo Y
    2. Telese F
    3. Chen H
    (2023) NCBI Gene Expression Omnibus
    ID GSE214388. Opiate responses are controlled by interactions of Oprm1 and Fgf12 loci in the murine BXD family: Correspondence to human GWAS findings.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE214388

References

    1. Crabbe JC
    (1998)
    Provisional mapping of quantitative trait loci for chronic ethanol withdrawal severity in BXD recombinant inbred mice
    The Journal of Pharmacology and Experimental Therapeutics 286:263–271.
  54. Book
    1. Mulligan MK
    2. Mozhui K
    3. Prins P
    4. Williams RW
    (2017) GeneNetwork: a toolbox for systems genetics in
    In: Schughart K, Williams RW, editors. Systems Genetics: Methods and Protocols. New York, NY: Springer. pp. 75–120.
    https://doi.org/10.1007/978-1-4939-6427-7
  74. Book
    1. Theriot J
    2. Sabir S
    3. Azadfard M
    (2022)
    Opioid AntagonistsStatPearls
    Treasure Island (FL): StatPearls Publishing.

Article and author information

Author details

  1. Paige M Lemen

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1774-0385

  2. Yanning Zuo

    University of California, San Diego, San Diego, United States
    Contribution
    Data curation, Formal analysis, Investigation, Writing – original draft
    Competing interests
    No competing interests declared

  3. Alexander S Hatoum

    Washington University School of Medicine, St. Louis, United States
    Contribution
    Data curation, Formal analysis, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8002-7267

  4. Price E Dickson

    Marshall University, Huntington, United States
    Contribution
    Data curation, Investigation
    Competing interests
    No competing interests declared

  5. Guy Mittleman

    Ball State University, Muncie, United States
    Contribution
    Data curation, Investigation
    Competing interests
    No competing interests declared

  6. Arpana Agrawal

    Washington University School of Medicine, St. Louis, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  7. Benjamin C Reiner

    University of Pennsylvania, Philadelphia, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  8. Wade Berrettini

    University of Pennsylvania, Philadelphia, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared

  9. David George Ashbrook

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Data curation, Formal analysis, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7397-8910

  10. Mustafa Hakan Gunturkun

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared

  11. Xusheng Wang

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Formal analysis, Writing – original draft
    Competing interests
    No competing interests declared

  12. Megan K Mulligan

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  13. Caleb J Brown

    Nash Family Department of Neuroscience and Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  14. Eric J Nestler

    Nash Family Department of Neuroscience and Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, United States
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7905-2000

  15. Francesca Telese

    University of California, San Diego, San Diego, United States
    Contribution
    Data curation, Formal analysis, Supervision, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared

  16. Robert W Williams

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    labwilliams@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8924-4447

  17. Hao Chen

    University of Tennessee Health Science Center, Memphis, United States
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    hchen3@uthsc.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2680-6921

Funding

National Institute of Development Administration (U01DA053672)

  • Robert W Williams
  • Hao Chen

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Ethics

This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (#23-0418) of the University of Tennessee Health Science Center.

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.108845. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2025, Lemen et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 653
    views
  • 24
    downloads
  • 2
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Citations by DOI

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Paige M Lemen
  2. Yanning Zuo
  3. Alexander S Hatoum
  4. Price E Dickson
  5. Guy Mittleman
  6. Arpana Agrawal
  7. Benjamin C Reiner
  8. Wade Berrettini
  9. David George Ashbrook
  10. Mustafa Hakan Gunturkun
  11. Xusheng Wang
  12. Megan K Mulligan
  13. Caleb J Brown
  14. Eric J Nestler
  15. Francesca Telese
  16. Robert W Williams
  17. Hao Chen
(2026)
Acute opioid responses are modulated by dynamic interactions of Oprm1 and Fgf12
eLife 14:RP108845.
https://doi.org/10.7554/eLife.108845.3

Share this article

https://doi.org/10.7554/eLife.108845

Categories and tags

Research organisms

Further reading

    1. Genetics and Genomics
    2. Neuroscience

    Opioid Use Disorder: Identifying new targets in the fight against opioids

    Fernanda Laezza
    Insight

    Experiments reveal that a time-dependent epistatic interaction influences how mice respond to opioids, and that intracellular fibroblast growth factors also influence opioid sensitivity.