Opiate responses are controlled by interactions of Oprm1 and Fgf12 loci in rodents: Correspondence to human GWAS findings

We mapped high-precision time-series data (15 min bins for 3 hours) generated for ~ 700 adult BXD mice across 105 morphine- and naloxone-related traits using new sequence-derived marker maps and a linear-mixed model. We confirm a previously mapped sex-independent effect of initial locomotor responses to morphine (50 mg/kg ip) that maps precisely to Oprm1 on chromosome (Chr) 10, with the linkage score reaching −log10P of ~12.4 (with a high B allele) at 75 min and exhausted by 160 min. We detected a new modulator of opiate locomotor activation in both sexes on Chr 16, with a peak linkage that climbs from 105 through to 180 min after injection. This locus includes one compelling candidate—fibroblast growth factor 12 (Fgf12). We also detected a strong, but transient epistatic interaction between these two loci. Single nuclei transcriptomic analyses in rats demonstrates that expression of Oprm1 and Fgf12 mRNA covary in one specific subtype of Drd1 medium spiny neurons. Our Bayesian network analysis identified that a cascade of MAP kinases—Mapk8ip2, Map3k11, and Map3k12—are part of the Oprm1–Fgf12 network. This is the first demonstration of a time-dependent epistatic interaction modulating drug response in mammals with interesting mechanistic implications. Analysis of OPRM1 and FGF12 gene networks in human GWAS data highlights enrichment of signals associated with substance use disorder.


Introduction
Genetic factors play significant roles in the risk for opioid use disorder (OUD) (Crist et al., 2019).
Classical twin studies have estimated that approximately half of the risk for OUD is genetic (Kendler et al., 2003;Tsuang et al., 1998). However, it has been challenging to discover genetic variants that influence the risk for OUD in humans or in animal models (Crist et al., 2019;Mistry Results QTL Mapping of Morphine-Induced Locomotor Responses. A total of 64 members of the BXD family were used in balanced male and female cohorts of roughly 6 cases per sex and strain, as described in Philip et al. (2010). Locomotion was recorded for 180 min after an acute morphine injection, with data for each 15 min combined into one bin. Data for the 10th time bin (135-150 min) were lost and were not included in our analysis. Variation in locomotor responses showed a skewed distribution (Supplementary Figure 1). We therefore used the quantile normalized form of these data (Supplementary Figure 2) in our analysis.
The results of QTL mapping for each time bin was shown in Figure 1 (Genome-wide significance level was 3log 10 P value of 3.77 for the BXD family). We found, during the 03135 min intervals, a significant QTL on Chr 10 at~7.6 Megabases (Mb) in males and at~8.6 Mb in females ( Figure 2a-b), with maximal 3log 10 P of 10.06 for males and 9.60 for females. This QTL was also reported in the original analysis by Philip et al. (2010). For the first 135 minutes, this peak remained as the most significant QTL. However, at the 1503165 min bin of locomotion data, we detected a second significant QTL, which was located on Chr 16 (~25330 Mb, Figure   3a-b). The highest 3log 10 P was 4.28 for males and 10.56 for females. This is a novel association not identified by Philip et al (2010).
Males and females showed similar patterns of QTL results at both Chr 10 and Chr 16 locus.
The time bins where the associations were found and locations of the loci were almost identical. This was confirmed by strong positive correlation in morphine-induced locomotion response between sexes. Examples at the 45-60 min and 165-180 min time bins were shown in Figure   1c-d.
Candidate Gene Identification for Morphine-Induced Locomotor Response. MOR is expressed in multiple brain regions (Sivalingam et al., 2020) and is involved in opiate-induced reward and locomotor response (Contet et al., 2004). Oprm1, which encodes the opioid receptor mu 1, is located in the Chr 10 locus (Figure 2c-d) and is the most likely candidate gene.
Genenetwork contains many brain region-specific gene expression data sets. We found several datasets where the Chr 10 locomotion locus also regulated the expression of Oprm1, i.e. a cis-eQTL. For example, Trait 1451709_at from NAc is shown in Figure 2e, and Trait 1441791_at from hippocampus is shown in Figure 2f. This evidence further supports Oprm1 as the candidate gene for the QTL on Chr 10. Examining the haplotypes, we found the D allele has higher expression in both the NAc (Figure 2e) and hippocampus (Figure 2f). We next examined the QTL located on Chr 16 ( Figure 1). This is a novel QTL, with Fgf12 as the possible candidate gene (Figure 3c-d, 4). 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. We found cis-eQTL for Fgf12 in striatum (Figure 3 e,g, peak 3log 10 P of 4.43, Trait ID: 4886596) and ventral tegmental area (VTA) (Figure 3 f,h, peak 3logP of 4.55, Trait ID: 1426432_at). We further scanned the phenotypes associated with the Chr 16 locus (i.e., PheWAS) using systems-genetics.org (Li et al., 2018), and found that Fgf12 significantly affects responses to various drugs, including morphine and naloxone ( Figure 5). Together, these data provided strong support for Fgf12 to be the candidate gene for the Chr 16 locus. Both candidate genes are summarized in Table 1.

Potential regulatory mechanisms of gene expression. We queried the variant browser in
GeneNetwork and found that all SNPs located in Oprm1 and Fgf12 regions in the BXD family were non-coding with a possible impact on transcript level as opposed to protein function. We further examined structural variants located in the vicinity of the Oprm1 and Fgf12 genes identified using linked-read genomic sequence data in GeneNetwork (Watson and Ashbrook, 2020). We found two large deletions within the 5' UTR of Oprm1 at about 6.2 and 6.7 Mb on Chr 10. Although neither of these variants overlap with known cis-regulatory elements, there is a  Supplementary Table 3. These data provided likely mechanisms linking genetic variants to expression changes.
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 (0390min) and near Fgf12 at later time points (+90 min). We tested for a possible epistatic interaction between these distinct loci across the entire time series. Figure 7 highlights the epistatic interaction for all four possible two-haplotype combinations on Chrs 10 and 16 (BB, BD, DB, DD) in both sexes. This interaction effect is strongest in the 45-60 min period4precisely at the midpoint between the early-additive Oprm1 effect and the late-additive Fgf12 effect (Fig 8c-g, e.g., male trait BXD_11331). Those BXDs that inherit the B haplotype of Oprm1 (defined as rs29339674, 6.7 Mb) when combined with the D haplotype of the Fgf12 locus (defined as rs4169220 at 31.6 Mb)4have unexpectedly high levels of locomotion compared to all of the other two-locus combinations. The LOD score of the epistatic component is in the range of 2.5-4.8. (MF = trait BXD_11845, LOD=4.2; F trait = BXD_11588, LOD=4.8; M trait = BXD_11331, LOD=2.5) and is significant by permutation testing at a genome-wide threshold of 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 LOD score 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) has an epistatic effect with a LOD of 3.0; an additive effect near Oprm1 with a LOD 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 LOD of 11.4.
By 75390 minutes (trait BXD_18847) the interaction effect has dropped to a LOD of 2.2 and the model has a cumulative LOD of 9.40. In marked contrast, over the two intervals from 90 to 120 minutes we do not detect significant epistasis. At best, the interaction LOD is only in the range of 1.532.0. In the 120-135 and 150-165 minute 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 LOD scores of 2.1 and 4.2, respectively.
Finally, in the last interval4165-180 minutes (trait BXD_11585)4we 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 LOD that is 4.9.
Cell-type specific co-expression analysis of Oprm1 and Fgf12 in the rat nucleus accumbens 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 core obtained from two female heterogeneous stock (HS) rats (Carrette et al., 2021). After filtering out low-quality nuclei and potential doublets, a total number of 4,495 high-quality nuclei transcriptomes remained with a median number of 3,363 transcripts (unique molecular identifiers) and 1,861 genes detected per nucleus. We then performed SCT transformation , dimensional reduction, and clustering using Seurat (Stuart et al., 2019) and identified 17 cell-type clusters

Constructing a Causal Bayesian Network.
We hypothesize that genetic variants in the Chr 10 and Chr 16 QTL regions are driving the Oprm1 and Fgf12 differential expression in the BXD family and influencing morphine-induced locomotor response. We used literature to guide us in identifying potential signaling molecules that can moderate the interaction of these two genes (Bhushan et al., 2017;Buchsbaum et al., 2002) to be included in our causal network model. For example, the FGF12 gene (FHF1 protein) (Sochacka et al., 2020) and MAPK8IP2 (Islet Brain-2 or IB2) are linked to sodium channels (SCN2A=Nav1.2 and SCN8A=Nav1.6) as part of a scaffold. The scaffold protein that MAPK8IP2 encodes is predicted to play a role in modulating the c-Jun amino-terminal kinase signaling pathway. It has been demonstrated that this protein interacts with and controls the activity of the kinases MAPK8/JNK1 and MAP2K7/MKK7 (O'Leary et al., 2016).
Using this prior literature and the available gene expression data from BXD mice, we tested Scn2a, Scn8a, Map3k11, Map3k12 and Mapk8ip2 as candidates using a BN framework available in GeneNetwork ( Figure 9). Expression levels of these genes in relevant brain regions of the BXD mice were retrieved from Genenetwork, together with genotype and behavior data.
These data were then transmitted directly to the Bayesian Network Web Server (BNWS) (Ziebarth and Cui, 2017). In the BNWS, tiers are used to organize the nodes (i.e. loci, gene expression, phenotypes, etc.) in a network into different layers. Each tier represents a different level of abstraction or complexity, with nodes at higher tiers representing more general or abstract concepts, and nodes at lower tiers representing more specific or concrete concepts.
We used 4 tiers for chromosomal location (green), candidate gene (blue) , signaling molecules (purple), and phenotype (red). The phenotype and genotype data we provided supported a causal model (

QTL Mapping of Naloxone induced withdrawal responses.
After the morphine locomotion tests were complete, 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). We found 11 genome-wide significant QTLs among 5 naloxone induced withdrawal behavioral responses ( Figure 6). The genomic location for these QTLs are listed in Table 2. We also identified 61 suggestive (3log 10 P at or above 2.34) QTLs in males and 65 in females (Supplementary Table 2 We also evaluated the biological function of potential candidate genes using a literature mining tool Genecup (Gunturkun et al., 2022). We found multiple candidate genes in the QTL on Chr 14 for naloxone induced jumps 0315 min after injection, which include Slc7a7 and Slc7a8.
Both genes belong to the solute carrier family 7. Solute carriers are transmembrane transporters that help drugs and other substances to cross biological membranes (Hu et al., 2020). These candidate genes are listed in Table 1.

Integrating murine data with human GWAS results.
There is support for our results in human GWAS of OUD. OPRM1 variant rs179971 has been associated with OUD in the first well-powered meta-analysis (Zhou et al., 2020) and confirmed in a recent study (Hatoum et al., 2022a). A human GWAS of opioid cessation reported a signal for FGF signaling pathway, where the effect of FGF12 was significant (p=0.0015, odd ratio=1.24) (Cox et al., 2020). In gene-based analyses by Deak et al. (2022) for OUD, there was a significant signal at FGF12 (p = 0.006).
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 genes that overlap with FGF12 and OPRM1 in the GTEx sample (GTEx Consortium, 2013) (GTExv8). We then examined whether there is gene-based enrichment for this network via MAGMA gene-set enrichment analysis (de Leeuw et al., 2015) in the GWAS by Zhou et al.(2020). Our initial analysis was performed for accumbens, putamen, caudate and adipose (as a negative control) for the FGF12 lists. This gave us only 8 tests to correct for (OUD and addiction GWAS in the four tissues). For OUD specifically, there was no significant enrichment for the original 500 gene lists, likely due to the large size of the lists and the lower power we have in SUDs GWAS. In the GWAS by Deak et al. (2022), where FGF12 was nominally significant, there was nominal significance for putamen and caudate but not for accumbens or adipose tissue, and neither significant result passed correction for these 8 tests (p-values at~0.01).
Upon realizing that this data-set had very low power and our gene list cut-off was too large, we decided to specify our list. The network of genes was defined by examining genes correlated with OPRM1 and FGF12 found in the specific tissue (see methods). Consistent with SNP and gene-based results, there was nominal enrichment for the addiction risk factor (Hatoum et al., 2022b) (p=0.038), but no enrichment for the network in the OUD GWAS (p=0.493).

Discussion
We analyzed time-dependent behavioral responses to morphine and naloxone collected from the BXD family of mice over a decade ago (Philip et al., 2010) using WGS-based genetic markers and linear mixed models. We confirmed an association between locomotor response and a region on Chr 10 that overlaps Oprm1. We discovered a novel association on a Chr 16 locus that overlaps Fgf12 in both sexes. These two loci had a significant but transient epistatic interaction between 45390 min after morphine injection. Transcriptomic data from NAc in rats showed 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.
The proximal Chr10 locus associated with early phase locomotor response to morphine contains the Oprm1 genes. 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).
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. For example, when MOR was deleted from D1 neurons, mice showed hyperlocomotion in response to the opioids. But when MOR was deleted from A2a neurons, mice showed an increased movement in response to opioids (Severino et al., 2020). These data indicated that Oprm1 is a strong candidate gene for the Chr 10 locus associated with morphine induced locomotion response.
The locomotion data showed BXDs with the D allele at the Oprm1 locus were less sensitive to a 50 mg/kg dose of Morphine (Figure 7). They demonstrate a blunted motor activity response similar to what was seen in low doses of morphine (5 to 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 showed 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 are in agreement with those observed in the parental strains of BXD mice, where locomotor response to morphine from C57BL6 was more pronounced and lasted longer compared to those from DBA2 mice (Murphy et al., 2001). This was accompanied by greater morphine-induced dopamine release in the NAc in C57BL6 than in DBA2 mice (Murphy et al., 2001). Therefore, it is likely that differences in morphine-induced dopamine release is involved in the differential locomotor response to morphine in the BXD mice.
The association between the Chr 16 locus (Figure 1 a,b) and morphine induced locomotion 1503180 minutes after injection is a novel finding. Over the past decade, with the expanded BXD lines, refinements in genetic markers and improved mapping methods have notably improved the possibility of detecting significant association between genotypes and phenotypes; QTLs previously considered non-significant for traits might now have genome-wide significance in the larger database (Watson and Ashbrook, 2020).
While there are several other 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 3 e,f). Fgf12 is a member of the FGF family. These growth factors regulate many cellular functions, including proliferation, survival, migration, and differentiation (Yun et al., 2010). Studies have implicated members of the FGF family, such as Ffg2 (Even-Chen and Barak, 2019a), FGF receptor 1 (Even-Chen and Barak, 2019b), and FGF receptor 2 (Blackwood et al., 2019) in substance abuse. More recently, it was demonstrated that administration of an FGF21 analog in non-human primates decreased alcohol consumption by 50% via an amygdalo-striatal circuit (Flippo et al., 2022). Unlike these secreted FGFs, Fgf12 is an intracellular protein that serves as cofactors 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 been shown to play a role 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 prenatal strains of the BXD mice (Murphy et al., 2001). Further, mice lacking Drd1 were not responsive to acute cocaine (Karlsson et al., 2008;Xu et al., 1994)  been the subject of intense research. A recent snRNA-seq study of the NAc in rat brains also detected this specific subtype, which was labeled Grm8-MSN subtype (Savell et al., 2020).
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.
We further explored the signaling network associated with Oprm1 and Fgf12. Kinases play key roles in the behavior and expression of both candidate genes (Bhushan et al., 2017;Buchsbaum et al., 2002). More specifically, FGF12 was shown to interact with the MAP kinase scaffolding protein, IB2 (MAPK8IP2) (Schoorlemmer and Goldfarb, 2001). In addition, FGF12 also interacts with voltage gated sodium channel NAV1.5 (Ornitz and Itoh, 2015). We included gene expression data of several of these MAP kinases and sodium channels in our BN. The Beyesian network searching algorithm excluded several of the kinases and all the sodium channels. The resulting network (Figure 9) indicated that Mapk8ip2, Mak3k11 (Mlk3) and Map3k12 (Dlk) form a highly interactive network mediating the interaction between Oprm1 and 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, warranting the need for future targeted studies to understand relationships between these genes and their associated phenotypes. 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, we also found significant QTL for the differences in locomotion before and after naloxone injection  Figure 6). 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 (-log 10 P=5.34 in males and -log 10 P=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.
Lastly, 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). All of these specifics across OUD stages also come with continually evolving definitions, and, therefore, also contributes to uncertainty and potential loss of significant results.
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 mammals4a 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 data sets 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.

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 (839 weeks) reared at the Oak Ridge National Laboratory in 200732008. 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 h and were exported as the sum of 15 min bins. Then, each mice 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 minutes. Naloxone's effects on locomotor activity levels, and signs of withdrawal were recorded for 15 minutes, including the number of jumps, fecal boli, urine puddles, wet dog shakes, instances of abdominal contraction, salivation, ptosis and abnormal posture were recorded for 15 minutes). The number of animals used per strain was 7.8 (mean ± SD) for females and 6.6 for males. Testing occurred at about 8 to 9 weeks after birth.

QTL Mapping in GeneNetwork
We reanalyzed 105 phenotypes from the Philip et al. (2010) study available in genenetwork.org.
The trait IDs and their most significant loci are shown in Supplementary Table 1. 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 algorithm with the LOCO option, which corrects family structure, to perform genetic mapping. 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 3log 10 P 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. Phenome-wide association studies (PheWAS) examine the relationship between a SNP and their association to collections of phenotypes. We performed a PhWAS on Oprm1 and Fgf12 using systems-genetics.org (Li et al., 2018).
The network for a gene of interest was defined by the top 500 genes it was correlated with in behavioral relevant brain regions in the GTEx database (GTEx Consortium et al., 2017) (e.g., substantia nigra, spinal cord, putamen, pituitary gland,hypothalamus, hippocampus, frontal cortex, cerebellum, cerebellar hemisphere, cerebellar cortex, and NAc, etc.). Control networks were generated using additional tissues that are not likely involved in addiction (e.g., whole blood, subcutaneous tissue, transformed fibroblast cells, etc.). These networks were calculated using transcriptome data archived in GeneNetwork. Functional enrichment of these networks were 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 2 female HS rats were obtained from the oxycodone tissue repository at UCSD (Carrette et al., 2021). The brain tissues were collected from rats sacrificed after 4 weeks of prolonged abstinence from oxycodone intravenous self-administration. The rats used in this study displayed a low addiction index-like behavior, which is calculated using several behavioral measures, as detailed in Carrette et al, (Carrette et al., 2021). Brain tissue was extracted and snap-frozen (at −30°C). Cryosections of~500 microns (Bregma 2.28-0.72 mm) were used to dissect the NAc core punches on a −20°C frozen stage. Punches from 3 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 6,000 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 by Qubit dsDNA HS Assay Kit (Thermo-Fischer 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
For this study, we only analyzed data generated by the RNA library type (GEX).
Raw sequencing data were converted to FASTQ format using bcl2fastq (Illumina). The FASTQ data were firstly aligned to the Rattus Norvegicus mRatBN7.2/rn7 genome and then aggregated into one sample file using CellRanger arc 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 16000, 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 as 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  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 correlation was computed in R 4.2.1 using cor.test and p<0.05 was considered significant.

Human Translational Data
First, 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 (GTEx Consortium et al., 2017) and GeneNetwork and then extracted out the shared genes with a liberal cut-off 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 found to be 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.

Causal Model by Bayesian Network
Lastly, we created a causal model of how these genes collaborate as underlying mechanisms behind the locomotion behavioral phenotype. To do this, we used the BNWS for Biological Network Modeling (https://bnw.genenetwork.org/sourcecodes/home.php) (Ziebarth et al., 2013).
To generate the model we first clicked on the "Create a Model" tab. On the "Create a Model" page, we uploaded our BXD morphine raw data file.
Once the data file was uploaded, we specified our model parameters. This included the number of nodes needed in the network, and a few other relevant parameters. We separated our nodes into 4 tiers. to label chromosome location, candidate gene, signaling molecules, and phenotype each as a separate tier. These are used to organize the nodes into different layers.
Each tier represents a different layer of complexity. Once specified, we clicked the "Create Model" button to start building our BN.

Data availability
snRNA-seq data has been deposited into NCBI GEO (accession number: GSE214388)