Behavioral and transcriptomic markers of susceptibility to escalate fentanyl intake.

Behavioral and transcriptomic markers of susceptibility to escalate fentanyl intake - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Search Log in Primary site navigation Search Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Transl Psychiatry . 2026 Apr 22;16:299. doi: 10.1038/s41398-026-04034-1 Behavioral and transcriptomic markers of susceptibility to escalate fentanyl intake Jack Keady 1, # , Richard Charnigo 2, # , Jakob D Shaykin 1, 3, # , Emily R Prantzalos 1 , Mengfan Xia 4 , Emily Denehy 3 , Cody Bumgardner 5, 6 , Justin Miller 5, 6, 7, 8 , Pavel Ortinski 4 , Michael T Bardo 3 , Jill R Turner 1, ✉ Author information Article notes Copyright and License information 1 College of Pharmacy, University of Kentucky, Lexington, KY 40506 USA 2 Department of Biostatistics, University of Kentucky, Lexington, KY 40506 USA 3 Department of Psychology, University of Kentucky, Lexington, KY 40506 USA 4 Department of Neuroscience, University of Kentucky, Lexington, KY 40506 USA 5 Department of Pathology and Laboratory Medicine, University of Kentucky, Lexington, KY 40506 USA 6 Division of Biomedical Informatics, Department of Internal Medicine, University of Kentucky, Lexington, KY 40506 USA 7 Sanders-Brown Center on Aging, University of Kentucky, Lexington, KY 40506 USA 8 Department of Microbiology, Immunology, and Molecular Genetics, University of Kentucky, Lexington, KY 40506 USA ✉ Corresponding author. # Contributed equally. Received 2025 May 13; Revised 2026 Mar 12; Accepted 2026 Mar 30; Collection date 2026. © This is a U.S. Government work and not under copyright protection in the US; foreign copyright protection may apply 2026 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ . PMC Copyright notice PMCID: PMC13234145  PMID: 42020354 Previous version available: This article is based on a previously available preprint posted on bioRxiv on January 29, 2025: " Behavioral and genetic markers of susceptibility to escalate fentanyl intake ". Abstract The “loss of control” over drug consumption, known as escalation of intake in opioid use disorder (OUD), is well-established in preclinical rodent models. However, little is known about how antecedent behavioral characteristics, such as valuation of hedonic reinforcers prior to drug use, influence fentanyl intake trajectories. Moreover, it is unclear if distinct escalation phenotypes are driven by transcriptomic markers predictive of OUD susceptibility. Male and female Sprague-Dawley rats (n = 72) were trained in a sucrose reinforcement task using a progressive ratio schedule. Individual differences in responsivity to sucrose were hypothesized to predict escalation of fentanyl intake. Rats underwent daily 1-h acquisition sessions for i.v. fentanyl self-administration (2.5 µg/kg; FR1) for 7 days, followed by 21 6-h escalation sessions, then tissue from prefrontal cortex was collected for RNA sequencing and qPCR. Latent growth curve and group-based trajectory modeling were used, respectively, to evaluate the association between sucrose reinforcement and fentanyl self-administration and to identify whether distinct escalation phenotypes can be linked to gene expression patterns. Sucrose breakpoints did not predict fentanyl acquisition nor change during escalation but did predict fentanyl intake on the first day of extended access. Permutation analyses found no associations between behavior and single gene expression, either overall or within our ascertained phenotypes. However, weighted genome correlation network analysis (WGCNA) and gene set enrichment analysis (GSEA) determined several gene modules linked to escalated fentanyl intake, including genes coding for voltage-gated potassium channels, calcium channels, and excitatory synaptic signaling. Transcription factor analyses identified EZH2 and JARID2 as potential transcriptional regulators associated with escalated fentanyl intake. Further, these modules were enriched for genome-wide association study (GWAS) term categories relating to substance use disorders. Escalation of opioid intake largely differs from motivation for natural rewards like sucrose. Further, the gene networks associated with fentanyl escalation suggest that engagement of select molecular pathways were associated with “addiction prone” behavioral endophenotypes, potentially representing druggable targets for OUD. Our extended in silico analysis found that gene networks associated with the “addiction prone” high escalating rats were enriched for genes with known risk alleles for substance use disorders and identified transcription factors that may regulate these networks, highlighting the importance of integrating findings from translational preclinical models, supporting patient-centered treatment options for OUD. Subject terms: Pharmacogenomics, Predictive markers, Addiction, Molecular neuroscience Introduction Escalation of intake is a prominent feature of opioid use in humans. Diagnostic criteria of opioid use disorder (OUD) incorporate measures of increased drug consumption, such as changes in frequency or dose, as a factor directly impacting severity of opioid use across time. While increased opioid use may be perceived as tolerance to a unit dose of the drug, clinical data indicate that the small changes in pharmacokinetic parameters for orally administered opioids after chronic use [ 1 , 2 ] are insufficient to explain the 10–20 fold increase of drug intake common with escalation [ 3 ]. In the meantime, escalated intake of opioids has been linked to a significantly elevated risk of developing both opioid and non-opioid substance use disorders, as well as other adverse outcomes, in a recent retrospective study of over 50,000 chronic pain patients in the Veterans Healthcare Administration system [ 4 ]. These considerations suggest that escalation of opioid use represents a phenotypic signature that can predict substance use in the future and may potentially capture other traits characteristic of risk. Escalation of drug intake can be observed in rodent models under conditions of extended access to opioids in self-administration (SA) sessions. Extended access SA sessions last anywhere from 3–23 h, similar to the time window of opioid use in people with OUD [ 5 ], whereas short (1–2 h) SA sessions serve as a useful assay for controlled, but limited intake [ 6 , 7 ]. Moreover, extended access SA sessions are advantageous in their ability to produce increased heroin- and fentanyl-induced reinstatement and heightened signs of withdrawal compared to rodents with short access [ 8 – 10 ]. Many factors may influence escalated use, including intermittency of access, hedonic value changes, withdrawal symptoms, habit formation, and environmental (e.g., context or social) events [ 11 ]. While contributing to significant heterogeneity in behavioral patterns associated with fentanyl use, these observations also highlight the possibility that behavioral characteristics present prior to initiation of fentanyl use may impact the emergence of OUD. Several preclinical studies have assessed antecedent sensitivity to non-drug reward to predict later drug reinforcement across several drugs of abuse, but have reported mixed findings [ 12 – 16 ]. Thus, it is our goal to ascertain an antecedent model with high predictive validity to identify heterogeneous subpopulations that may be at risk for high opioid intake, prior to initiation of fentanyl use. The prefrontal cortex (PFC) is an important site of both neuronal activity underlying addictive behavior and changes in gene expression induced by opioid exposure [ 17 – 19 ]. Widespread transcriptional activation in the PFC serves as a regulatory response following chronic opioid exposure [ 20 ]. Moreover, while acute opioid exposure has acute pharmacological actions throughout the central nervous system (CNS), counter-adaptational transcriptomic signatures persist in the PFC [ 21 ]. Therefore, analyses of PFC transcriptomic changes provide a lens through which individual fentanyl escalation phenotypes may be captured, with the goal of uncovering the neural mechanisms underlying the “loss of control” phenomenon. In the present study, we aimed to: 1) examine whether individual differences in motivation for non-drug reward influence the likelihood and extent of subsequent fentanyl use, 2) characterize individual differences in fentanyl escalation, and 3) determine the transcriptomic signatures associated with escalation-prone behavioral phenotypes. Methods Subjects Male and female Sprague-Dawley rats (PND 50–55; n = 72; n = 36 male; n = 36 female) were individually housed in a colony room with ad libitum water access and were maintained on a 12-h/12-h light/dark cycle with lights on at 0700 h. Prior to sucrose or fentanyl self-administration sessions, rats were moved from the colony room to the testing room and weighed, allowing for time to acclimate to the testing room. All experimental procedures were followed in accordance with the NIH Guide for the Care and Use of Laboratory Animals (8 th edition, 2011) and approved by the University of Kentucky Institutional Animal Care and Use Committee. Sample sizes were based on prior fentanyl self-administration and transcriptomic studies from our laboratory and others using comparable designs; no a priori statistical power calculation was performed. Animals were not randomized to experimental groups, as behavioral phenotypes were identified post hoc using data-driven latent growth curve and group-based trajectory modeling. Sucrose self-administration Self-administration experiments were conducted in ventilated, sound attenuating operant chambers, equipped with a house light, two nose-poke holes, a pellet dispenser and a food receptacle (Med-Associates Inc., East Fairfield, VT, USA). During the experiment (Fig. 1A ), rats had ad libitum access to water and were given 20 g of normal chow per day after each operant session to control for satiation given the caloric value of sucrose and maintain rats at ~85–95% of free-feeding weight [ 22 – 25 ]. All subjects were trained 1 hr daily on a fixed-ratio one (FR1) schedule of reinforcement between 0800-1000, in which each nose poke into a reinforced active hole delivered a single 45 mg sucrose pellet. Pokes into the inactive hole had no programmed consequences. Each sucrose pellet was followed by a 20 s timeout period during which house light went off, serving as a cue to signal the timeout, and nose pokes had no scheduled consequences. After 10 days of FR1 nose-poke training, the rats progressed to the FR3 schedule (three nose pokes for one sucrose pellet) for 3 days and then to FR10 (ten nose pokes for one sucrose pellet) for 8 days. Following FR training, the rats were placed on a progressive ratio (PR) schedule of reinforcement for 7 days, during which successive sucrose pellets could be earned according to an increasing number of reinforced nose-pokes based on the formula: [5e (pellet # * 0.2) ] – 5 [ 26 ]. This progressive ratio formula, originally developed by Richardson and Roberts [ 26 ], is widely used to assess motivation for reward, where e denotes Euler’s number, the base of the natural logarithm. The session ended when rats failed to reach the next nose-poke criterion within 1 hr. The final reinforced nose-poke criterion achieved was recorded as the “breakpoint” value for that session. On the last day of sucrose PR, the rats were given ad libitum chow. Fig. 1. Motivation for Sucrose as Predictor of Fentanyl Escalation. Open in a new tab A Experimental timeline. B , C Sucrose self-administration. B PND (62–82) Number of sucrose pellets that male (blue circles) and female (red circles) rats received over the course of FR 1 (day 1–10), FR 3 (day 11–13), and FR 10 (day 14–21) sessions. C Progressive ratio breakpoints for male and female rats across the 7 days of PR training. D , E Fentanyl self-administration. Rats underwent D 1 h acquisition and E 6 h escalation sessions for fentanyl (2.5 µg/kg/inf). F , G Latent growth curve model considering hypothetical rats that behave at exactly the average breakpoint (black line) over the last three days of sucrose PR, at one standard deviation above the average (red line), and at one standard deviation below the average (blue line). The relationship between sucrose PR and fentanyl infusions is visually appreciated as the magnitude of difference between the 3 growth curves [** indicates p < 0.01 as determined by Bonferroni’s multiple comparison tests, ## indicates p < 0.01 for the main effect of sex, ### indicates main effect of time. Data in Fig. 1B–E represent mean ± SEM]. Catheterization and fentanyl self-administration On PND 90, rats were anesthetized with a ketamine/xylazine/acepromazine (75, 7.5, 0.75 mg/kg) mixture at a 0.15 ml/100 g body weight (i.p.) and were implanted with chronic indwelling jugular catheters with a cannula secured to the skull with dental acrylic anchored by four screws (Fig. 1A ). Rats were given 5–7 days to recover after surgery and catheters were flushed daily with a solution containing gentamycin, heparin, and saline. Catheter patency was monitored throughout the study; if a rat showed an unexpected decrease in fentanyl intake for two or more consecutive sessions, patency was assessed via a post-session intravenous bolus of morphine (15 mg/kg in 0.2 ml), followed by a saline flush, with the expected rapid cataleptic response used as confirmation. To avoid potential carryover of learned motor patterns from sucrose responding, distinct operant responses were used for each reinforcer: nose pokes for sucrose and lever presses for fentanyl. Rats were then trained to self-administer intravenous (i.v.) fentanyl (2.5 µg/kg/infusion) using a two-lever operant procedure under FR 1 schedule of reinforcement in a two-lever operant conditioning chamber (ENV-001; MED Associates, St. Albans, VT) using methods described previously by Malone et al. [ 9 ]. Increasing FR values during sucrose self-administration were used to facilitate transition to a progressive ratio schedule for assessing motivation, whereas fentanyl self-administration remained at FR1 to allow evaluation of intake patterns over time. Briefly, adult male and female Sprague-Dawley rats (PND 97) were trained to self-administer fentanyl (2.5 µg/kg/infusion) in a 2-lever operant conditioning chamber using an autoshaping procedure. During the first 35 min of a 1-hr session, one lever (active lever) was extended into the chamber, followed 15 sec later by retracting the lever and delivering an infusion accompanied by 20-sec illumination of 2 cue lights, one located above each lever; pressing the lever during the 15-sec extension led to immediate delivery of the infusion and cue illumination. The intention of this procedure was to train rats to approach the active lever, while controlling for the number of drug exposures across rats. Two hr following the autoshaping session, rats were placed back into the chamber with both levers continuously extended for 1-hr and allowed to earn a fentanyl infusion and 20-sec cue light illumination with each active lever press (FR1); the 20-sec cue light illumination signaled a timeout period during which additional presses on the active lever were not reinforced. This protocol occurred across 7 consecutive days, at which point autoshaping sessions were omitted and the FR1 session was extended to 6 hr across 21 days. Session duration was increased to 6 h to model extended-access conditions that promote escalation of opioid intake and mimic clinically relevant patterns of excessive use. Of the initial cohort, 63 rats completed the acquisition phase and 58 rats remained for escalation analysis following exclusions for incompletion, catheter patency loss, or surgical failure. Behavioral data analysis General remarks Data analysis, visualization, and statistical calculations were carried out using SAS (Version 9.4), R (Version 3.4.1 or later), GraphPad Prism (version 9.4 or later) and MATLAB (R2022A); the PROC TRAJ add-on to SAS (trajv9-64) [ 27 , 28 ] was employed, as was the scatterplot3D package for R (version 0.3–44) [ 29 ]. Formal hypothesis tests used alpha = 0.05. Further details about methods appear in the supplemental file 1 . Quantifying sucrose preference Breakpoint values were recorded for each rat throughout the 7-day sucrose PR protocol. The average breakpoint over the final three days was calculated for each rat, as responding had stabilized and did not differ between sexes during this period. These three-day averages were converted to z - scores by subtraction of the overall mean (47.5) and then division by the overall standard deviation (29.8). Latent growth curve modeling of fentanyl acquisition and escalation We used PROC MIXED in SAS (Version 9.4) to fit a latent growth curve model to describe the expected number of fentanyl infusions as a quadratic function of time, with coefficients of the quadratic function (hereafter, “coefficients”) permitted to depend linearly on sucrose PR protocol z-scores . Fentanyl acquisition and escalation phases were modeled separately. Partial data due to a chewed or disconnected infusion line were included for all 7 days (acquisition phase) or all 21 days (escalation phase). The latent growth curve model included random effects for rat-specific adjustments to the coefficients with an unstructured covariance matrix. Group-based trajectory modeling (escalation phase) We used the PROC TRAJ add-on to SAS [ 27 , 28 ] with the MULTGROUPS option to perform group-based trajectory modeling simultaneously for three outcomes measured daily during the 21-day escalation phase: active lever presses resulting in fentanyl infusions (reinforced presses), nonreinforced active lever presses during signaled timeouts (futile presses), and non-reinforced inactive lever presses (inactive presses). These additional measures were included to capture behavioral processes such as locomotor activity and frustrative responding that may accompany escalation and assist in distinguishing subgroups; however, clusters were primarily differentiated by reinforced presses. Partial data were included from rats not having valid information for all 21 days. Futile presses and inactive presses were sufficiently non-normal in their distributions that we log transformed them (after the addition of 1, to avoid taking a logarithm of 0) before applying PROC TRAJ under censored normal model specifications (with upper limits of 500 for reinforced presses and 10 for the log transformed outcomes); however, for ease of interpretation, figures in the main manuscript and supplemental file 1 convey the original, non-transformed values. We found that a model with 4 groups exhibited a much better log likelihood than a model with 3 or fewer groups. Group-based trajectory modeling (acquisition phase) The group-based trajectory modeling methodology was mostly similar to that for the escalation phase, with the following exceptions. First, because non-normality was less pronounced in the acquisition phase data, the futile press and inactive press outcomes were not log transformed. Second, the censored normal model specifications were altered (with upper limits of 30 for reinforced presses, 70 for futile presses, and 40 for inactive presses). Third, we adopted a 5 group model. For further details on the classification of the five groups, please refer to the “Group-Based Trajectory Modeling” section in the Results. Cross-classifications into trajectory groups We used PROC FREQ in Version 9.4 of SAS with the EXACT option to test for association between trajectory group assignments derived from the escalation and the acquisition experiments; this test was based on the 56 rats assigned by trajectory modeling to non-anomalous (N > 2) groups in both experiments. Tissue collection, RNA-sequencing, secondary analysis Tissue collection and RNA extraction Eighteen h after the last self-administration session, the timepoint in line with when daily sessions were run, rats were euthanized and brain tissue was collected. Total RNA was extracted from the PFC using TRIzol (Invitrogen) and an RNeazy Mini Kit (QIAGEN) according to the manufacturers’ instructions. RNA samples were quantified using the NanoDropOne (Thermo Scientific) and were confirmed to be of good quality (RNA 260/280 > 2.0). RNA-sequencing and alignment Tissue from rats classified into fentanyl escalation Groups 1 (M = 5, F = 3), 2 (M = 4, F = 4), and 3 (M = 4, F = 4) were selected for RNA-sequencing. No rats from fentanyl escalation Group 4 were used, as there were only N = 2, and these rats were considered anomalous. These 24 samples were sent to GENEWIZ from Azenta Life Sciences (South Plainfield, NJ) for next-generation sequencing (standard RNA-seq service). The remaining samples were used as biological replicates for validation via RT-PCR: N = 4 - 7 per escalation group (Group 1: N = 4 [M = 1, F = 3]; Group 2: N = 7 [M = 6, F = 1]; Group 3: N = 5 [M = 1, F = 4]). Sample quality was verified (RNA integrity numbers >9) and RNA libraries were prepped with PolyA selection before sequencing on an Illumina NovaSeq 6000 S4 flow cell, generating 2 × 150 base pair (bp) paired-end reads. RNA - sequencing data was preprocessed and aligned to the reference genome using Galaxy (24.0.1.dev0), a web-based cloud computing platform for bioinformatic analysis, as previously described [ 30 – 33 ]. Briefly, fastq files from the NovaSeq 6000 were checked for quality using FastQC (Galaxy Version 0.74+galaxy0). Trim Galore (Galaxy Version 0.6.7+galaxy0) was used to remove adaptor content, low quality reads ( < 20 PHRED score), and read length ( < 20 bp) from fastq files, between 0.3% to 0.4% of reads were removed from each sample for not meeting quality requirements. Trim Galore output was aligned to the rat genome (rn6, downloaded from iGenome, Illumina) using STAR (Galaxy Version 2.7.11a+galaxy0). Average input read counts for the rat PFC samples were 29.39 million, with an average of 80.85% uniquely mapped reads. Gene counts were determined using featureCounts (Galaxy Version 2.0.3+galaxy2). Differential expression Differential RNA expression was determined using pyDESeq2 (0.4.9) at a false discovery rate (FDR) < 0.05, filtering out genes with less than 10 total reads per row, with fentanyl escalation groups treated as the between-group comparisons [ 34 ]. Sex was included as a variable in differential expression, but not treated as a between-group factor due to lack of sex differences observed behaviorally. Information regarding permutation-based testing for single gene differential expression can be found in Suppl. File 1 . Weighted gene co-expression network analysis We used the Python package PyWGCNA (2.0.3) to conduct Weighted Gene Co-Expression Network Analysis (WGCNA) on the sequencing results with minimum module size = 30, MEDissThres = 0.2, and a soft power 17 (Suppl. Fig. 1 )[ 35 ] on the corrected counts from pyDESeq2’s median of ratio. Functional enrichment of modules identified via WGCNA were performed via the GSEAPY (1.1.3) enrichr function [ 36 , 37 ], using the 2023 Gene Ontologies (Biological Process, Cellular Component, and Molecular Function)[ 38 , 39 ] as well as the 2023 GWAS catalog [ 40 ] with all gene names within a given module as input. Transcription factor analysis Binding Analysis for Regulation of Transcription (BART2 v = 2.1) was used to predict potential transcription factors regulating each module, using published mouse and human chromatin immunoprecipitation sequencing libraries [ 41 ]. Quantitative RT-PCR Quantitative reverse transcriptase PCR was performed as described previously [ 42 – 44 ] on PFC samples from escalation Groups 1,2, and 3. Briefly, RNA was isolated using the RNeasy Mini kit (Qiagen) and complementary DNA (cDNA) was synthesized from 500 ng of isolates RNA. qPCR reactions were assembled using Thermo Scientific Maxima SYBR Green master mix, cDNA, and primers diluted down to final concentration of 5 nM (Integrated DNA Technologies, Coralville, IA) designed via Primer3Plus. Levels of mRNA were determined using the 2 -ΔΔCT method [ 45 ], and genes of interest were normalized to the housekeeping gene TATA-Binding Protein ( Tbp ) after assessing the genes of interest at values in relation to Tbp and Glyceraldehyde 3-phosphate dehydrogenase ( Gapdh ). Gene expression values were normalized to fentanyl escalating Group 1, and primer sequences are shown in Table 1 . Genes chosen for RT-PCR validation were selected by generating a priority score by combining results from the differential expression analysis and WGCNA with the following formula: priority score = − l o g 10 adjusted P * | log 2 fold change | * log 10 ( connectivity ) 1 Where adjusted P is the adjusted P value and | log 2 fold change | is the absolute log2(FC) for a given gene between two escalation phenotypes, and connectivity is the sum of connection strength of a given gene with all other genes within the same module generated by the WGCNA package [ 35 , 46 ]. Priority scores were generated for genes in modules that correlate to our observed fentanyl escalation phenotypes. This formula was chosen to generate priority scores as it allows us to prioritize genes for RT-PCR validation that have differential expression, or trends towards differential expression, between observed phenotypes and are highly connected within a module correlated to our observed phenotypes. Therefore, the priority score reflects genes with potential role in mediating shifts in the transcriptomic profiles of the escalation phenotype. Table 1. Sequence of primers used in qPCR. Gene Forward Primer Reverse Primer Grin2a GATGACAACCACCTCAGCATT TTCTTCACATTCATCCCTTCG Grin2b TGCAAGGGGTTCTGTATTGAC AGTGATCCCACTGCCATGTAG Slc1a2 TGGTACAAGCCTGTTTCCAAC GCCGAAAGCAATAAAGAATCC Tnr TCCATGGGACAGTAACCTCTG TAGCTGGATAGAGCCTGGTGA Erbb4 TGGAGGAAAGCCCTATGATG CTGGGGGACCAAATATTCCT Kcna2 AGCTAGGAGAAGAGGCAATGG GTTTCCAGACAGAAGCTGACG Kcnma1 GTAACGTGGACACCCTTGAGA TGATCGTTGCCAGGAATTAAC Kcnh7 GTGGTAAGTGATCCTGCCAAA AGGAGACGTGCAGTCTTCAAA Kcnq3 TACCTGGTGGAGAAGGATGTG TCCTTCCCAGGTTTTAGGTGT Tbp GAGCTCTGGGATTGTACCACA CACCATTTTCCCAGAACTGAA Gapdh AGTCTACTGGCGTCTTCACCA TGCATTGCTGACAATCTTGAG Open in a new tab Grin2a Glutamate Ionotropic Receptor NMDA Type Subunit 2A; Grin2b Glutamate Ionotropic Receptor NMDA Type Subunit 2B ; Slc1a2 Solute Carrier Family 1 Member 2; Tnr Tenascin R; Erbb4 , Erb-B2 Receptor Tyrosine Kinase 4; Kcna2 , Potassium Voltage-Gated Channel Subfamily A Member 2; Kcnma1 , Potassium Calcium-Activated Channel Subfamily M Alpha 1; Kcnh7 , Potassium Voltage-Gated Channel Subfamily H Member 7; Kcnq3 , Potassium Voltage-Gated Channel Subfamily Q Member 3; Tbp TATA-Box Binding Protein; Gapdh , Glyceraldehyde 3-phosphate dehydrogenase. Differences in gene expression between escalation groups were tested using a one-way ANOVA in GraphPad Prism 10.0 software package (GraphPad Software, CA). RT-PCR data were normalized to escalation Group 1 (low-trajectory intake). To examine relationship between qPCR results and escalation phenotypes, a combined z-score for the biological replicates including reinforcers, inactive responses, and futile responding for each rat was calculated over the last 7 self-administration sessions using the following formula: Z total = X 1 − μ { 1 , 1 } σ { 1 , 1 } − X 2 − μ 2 , 1 σ 2 , 1 + X 3 − μ { 3 , 1 } σ { 3 , 1 } 3 2 Where X 1 is reinforcers, X 2 is futile presses, and X 3 is inactive presses, all averaged for each individual rat across the last 7 self-administration sessions. Then let μ { k , 1 } and σ k , 1 represent mean and standard deviation of the kth outcome (k = 1,2,3) for Group 1 biological replicates. SPSS Statistics software package 29 (IBM, NY) was then used to perform linear regressions treating qPCR data as the independent variables [9 genes of interest] and Z total as the dependent variable. The stepwise method with F value of 2 for entry and 1 for removal was used in these models. We also performed a principal component analysis via scikit-learn (1.5.1) on the biological replicates, with X 1 , X 2 , and X 3 . SPSS Statistics software package 29 (IBM, NY) was then used to perform linear regressions treating qPCR data as the independent variables [9 genes of interest] and the first principal component (PC1) as the dependent variable. The stepwise method with F value of 2 for entry and 1 for removal was used in these models. Additional supplementary methods can be found in Suppl. File 1 . Results Sucrose and fentanyl self-administration We completed behavioral profiling of sucrose self-administration on a PR schedule of reinforcement in 72 rats (36 males, 36 females; Suppl. Fig. 2 ) that were tested in 6 squads of 12 rats over an 8-month period. Figure 1B displays reinforcers earned across 21 days of sucrose self-administration, prior to the PR task. A mixed-effects ANOVA indicated a significant sex x time interaction (F (20, 1398) = 4.83, p < 0.001). To probe the interaction, a Bonferroni multiple comparisons test determined that male rats earned significantly more sucrose pellets on session 11 compared to females (p = 0.002). Figure 1C displays sucrose breakpoints as a function of time. A mixed-effects ANOVA indicated a significant sex x time interaction (F (6, 420) = 7.091, p < 0.001), and a Bonferroni multiple comparisons test indicated a significant difference in PR breakpoints for sucrose between males and females at session 2 (p = 0.001). The data indicated a decrease of sucrose breakpoints from the first to the second and third (for males) day of PR training followed by a more gradual decrease for males only over the remaining days. Average breakpoint values across the last 3 days of PR training for each rat were used for associations with fentanyl seeking in subsequent analyses as responding was stable and no differences were observed between sexes. After seven days of sucrose PR, the rats went on to fentanyl self-administration training, progressing from 1-h sessions for 7 days (“acquisition phase”, Fig. 1D ) to 6-h sessions for 21 days (“escalation phase”, Fig. 1E ). A total of 63 rats (n = 32 male; n = 31 female) completed all acquisition sessions and 58 rats (n = 28 male; n = 30 female) were included in escalation analyses; rats excluded from these analyses were due to catheter patency or surgical failure. Figure 1D displays fentanyl self-administration during the acquisition phase. A mixed-effects ANOVA indicated a significant main effect of sex (F (1, 61) = 7. 439, p = 0.008) and time (F (6, 366) = 40.66, p < 0.001), such that intake increased across the 7 days for both sexes. Figure 1E displays fentanyl self-administration during the escalation phase. A mixed-effects analysis indicated a significant effect of time (F (20, 952) = 22.68, p < 0.001), such that intake increased across the 21 days. Sex differences seen in acquisition were not observed in the escalation phase (F (1, 56) = 1.80, p = 0.19). Individual data points for reinforced, inactive, and futile presses can be found in Suppl. Fig. 3 . Latent growth curve modeling To evaluate whether sucrose PR performance predicted fentanyl seeking, we used latent growth curve modeling. We quantified a pattern of reinforced lever pressing for fentanyl as the initial number of presses at day 1 of training as well as the change in number of presses through the subsequent training days, with both allowed to depend on sucrose preference. We found that sucrose preference did not predict the overall pattern of reinforced presses during acquisition (p = 0.53; Fig. 1F ) but did predict the overall pattern of reinforced pressing during escalation (p = 0.04; Fig. 1G ). Interestingly, for the escalation phase, it was the initial number of reinforced presses (p = 0.03), rather than the change in reinforced pressing across sessions (p = 0.67), that was associated with antecedent sucrose PR performance. More specifically, latent growth curve modeling predicted 54 initial fentanyl infusions for a rat with average sucrose PR performance and a difference of ~6 reinforced presses for each standard deviation increase in sucrose PR performance at day 1 of escalation. On the last day of escalation, each standard deviation increase in sucrose PR breakpoint predicted only a marginally larger difference of 7 presses, highlighting that the change in lever pressing over the escalation phase was not associated with antecedent sucrose PR performance. We also pursued latent growth curve modeling to test for associations between sucrose PR and inactive or futile lever pressing during the acquisition and escalation phases of fentanyl self-administration. During the acquisition phase, sucrose PR did not predict inactive (p = 0.93) or futile (p = 0.84) lever pressing (Suppl. Fig. 4A, B ). Likewise, during the escalation phase, sucrose PR did not predict inactive (p = 0.61) or futile (p = 0.53) lever pressing (Suppl. Fig. 4C-D ). Thus, these results support a relationship specifically between sucrose PR and reinforced lever pressing resulting in fentanyl infusions during the first session of the escalation phase, but not between sucrose PR and futile presses on the active lever, nor between sucrose PR and inactive lever presses regardless of phase of self-administration training. Group-based trajectory modeling We also investigated whether sucrose PR performance could predict membership in groups of phenotypically similar rats. To this end, we used group-based trajectory modeling to probabilistically cluster rats based on the combination of reinforced, futile, and inactive lever presses. For the acquisition phase, group-based trajectory modeling identified five behavioral phenotypes, exhibited by n = 15 (Group 1), n = 21 (Group 2), n = 6 (Group 3), n = 20 (Group 4), and n = 1 (Group 5) rats respectively. Group 5 contained one rat and was excluded from further analysis over doubts that behavior of a single rat is generalizable, but lever press data from this group is nevertheless presented in Suppl. Fig. 5 . The remaining group numbers reflect fentanyl intake in the ascending order (Fig. 2A–C ). Rats in Group 1 maintained relatively stable levels of reinforced and futile presses, whereas rats in Group 4 reached a plateau of reinforced presses by day 5–6 of fentanyl acquisition with stable futile presses. Rats in Groups 2 and 3 displayed near linear increases in reinforced lever pressing at similar levels but differed on futile presses. While futile press numbers in Group 2 were stable, futile presses in Group 3 approximately tripled over the seven days of acquisition training (Fig. 2C ). As for inactive presses, Group 3 and Group 2 began at high and low levels, respectively, with intermediate (and nearly identical) numbers of inactive presses by Groups 1 and 4. All groups converged on similar numbers of inactive presses by the last day (Fig. 2B ). We found no significant association between sucrose PR z-scores and acquisition group assignment (p = 0.43; Fig. 2D ). Fig. 2. Group based trajectory modeling. Open in a new tab Group assignment for the acquisition phase was based on the following parameters: ( A ) reinforced lever presses, ( B ) inactive lever presses, ( C ) futile lever presses for rats placed in Group 1 (green), Group 2 (yellow), Group 3 (orange), and Group 4 (red). A separate model demonstrating group assignment for the escalation phase based on ( E ) reinforced lever presses, ( F ) inactive lever presses, ( G ) futile lever presses for rats placed in Group 1 (green), Group 2 (orange), Group 3 (red). Panels ( D ) and ( H ) show box plots of sucrose breakpoint z-scores by assigned acquisition and escalation groups. For the escalation phase, group-based trajectory modeling discovered four behavioral phenotypes, expressed by n = 17 (Group 1), n = 24 (Group 2), n = 15 (Group 3), and n = 2 (Group 4) rats, respectively. Again, because of dubious generalizability, Group 4 was largely excluded from further analysis (Fig. 2E–G ). The remaining group numbers were assigned to represent increasing fentanyl intake, such that escalation Groups 1, 2, and 3 can be characterized as low-trajectory intake rats, medium escalators, and high escalators, respectively. Reinforced presses increased nearly linearly for low-trajectory intake rats and medium escalators but plateaued in high escalators by day 16–17 (Fig. 2E ). Medium and high escalators, but not low-trajectory intake rats, increased futile lever pressing (Fig. 2G ). Whereas high escalators maintained a low level of inactive presses, low-trajectory intake rats and medium escalators increased their inactive lever pressing (Fig. 2F ). Assignment to escalation groups was not clearly associated with sucrose PR (p = 0.18; Fig. 2H ). However, the excluded Group 4 rats had the highest sucrose PR performance (Suppl. Fig. 6A ), had extreme numbers of futile presses and displayed parabolic patterns of reinforced, inactive, and futile presses (Suppl. Fig. 6B-D ). Three-dimensional scatterplots to show acquisition and escalation trajectories of reinforced, futile, and inactive pressing, as well as animations to depict these trajectories from different perspectives, can be found in Suppl. Fig. 7 . Finally, we examined whether sex was associated with group assignments generated by group-based trajectory modeling. Acquisition group membership was significantly associated with sex (p = 0.001), such that low fentanyl intake groups were predominantly male, whereas high fentanyl intake groups were predominantly female (Suppl. Fig. 8E ). In contrast, there was no significant association between sex and escalation group membership (p = 0.39) (Suppl. Fig. 8F ). Further detail on the examination of sex and latent growth curve modeling is provided in Suppl. File 1 . Relationship between acquisition and escalation groups Because group-based trajectory models of acquisition and escalation phases were constructed independently, we examined the association between acquisition and escalation group assignments. A Chi-Square test for association between acquisition and escalation group assignments yielded p = 0.14. Therefore, relatively lower or higher number of presses in the escalation phase was not significantly associated with lower or higher presses during acquisition, highlighting a dissociation between behavioral measures of fentanyl seeking under short and extended access conditions. Differential gene expression among escalation groups Bulk RNA sequencing identified 17337 unique genes across the three escalation groups (N = 8 animals/group) established by the group-based trajectory modeling. Genes with 10 or fewer total counts across all 24 samples were filtered out for pyDESEQ analysis, leaving 14277 identified genes. We examined differential gene expression in high escalators vs low-trajectory intake rats (Group 3 vs 1), high vs medium escalators (Group 3 vs 2), and medium escalators vs low-trajectory intake rats (Group 2 vs 1). Traditional sequencing analysis approaches indicated that 10 genes had significantly increased expression (FDR < 0.05 and log2(FC) > 0.5) in Group 3 compared to Group 1 rats (Fig. 3A ), including genes coding for voltage-gated potassium channels and NMDA receptors. An additional 37 genes were identified in Group 3 vs Group 1 with an FDR < 0.05 and |log2(FC)| ≤ 0.5 (32 upregulated and 5 genes down regulated). A single gene, Fgf1 (fibroblast growth factor1), was identified with statistically significant differential expression in Group 3 vs Group 2, (Fig. 3B ), and no genes were differentially expressed in Group 2 vs Group 1 at a statistically significant level (Fig. 3C ). Differential expression data for all identified genes in the three comparisons can be found in Suppl. File 2 . Fig. 3. Differential Gene Expression in Fentanyl Escalating Rats. Open in a new tab A Volcano plot of differently expressed genes (DEGs) in fentanyl escalation Group 3 vs Group 1. B Volcano plot of DEGs in fentanyl escalation Group 3 vs Group 2. C Volcano plot of DEGs in fentanyl escalation Group 2 vs Group 1. Adjusted p values in volcano plots were calculated according to the pyDESeq2 pipeline. Significantly upregulated genes (log2fold change > 0.5 and FDR < 0.05) are colored in red. Blue dotted lines are significance thresholds (-log10(0.05)) and green lines are cut-offs for the log2fold change (log2fold change =0.5 and −0.5). When using permutations to calculate unadjusted p-values, no statistically significant results remained for individual genes after adjustment for multiple comparisons (Suppl. Table 1 ). In comparing the three main escalation groups, there were 1056 genes (7.4%) with unadjusted p-values less than 0.05, suggesting that differential expression may exist (as 7.4% is greater than 5%). However, the smallest adjusted p-value for any individual gene was 0.36. In pairwise comparisons between escalation groups, adjusted p-values for all genes were at least 0.19. When looking at the raw numbers of fentanyl infusions during the last three days as a numeric measure of behavior (instead of groups derived from trajectory analysis), one gene had an adjusted p-value of 0.071, while all other genes had adjusted p-values of at least 0.271. Identification of gene modules associated with high fentanyl escalation We then performed a signed WGCNA to take into account directionality of the correlation in the gene co-expression network, and three modules (Whitesmoke – 282 genes, Brown – 335 genes, and Lightgrey – 144 genes) were significantly correlated to fentanyl escalation group with no modules correlated to sex (Fig. 4A, B ). Specifically, all three modules were positively correlated with Group 3 assignment [Whitesmoke (R = 0.5522, p = 0.005), Brown (R = 0.5418, p = 0.006), and Lightgrey (R = 0.4837, p = 0.017)]. The Brown module was negatively correlated with Group 1 assignment [R = −0.5361, p = 0.007], while the negative correlations of Whitesmoke (R = −0.3884, p = 0.061) and Lightgrey (R = −0.3378, p = 0.106) with Group 1 assignment did not reach statistical significance (Fig. 4B ). The Brown module included 6 of the 10 differentially expressed genes with FDR < 0.05 and log2(FC) > 0.5 and 13 of 37 differentially expressed genes with FDR < 0.05 and |log2(FC)| < 0.5 in the Group 3 rats compared to Group 1. The Whitesmoke included 2 of the 10 differentially expressed genes with FDR < 0.05 and log2(FC) > 0.5 and 14 of 37 differentially expressed genes with FDR < 0.05 and |log2(FC)| < 0.5 in the Group 3 rats compared to Group 1. The Lightgrey module contained no genes with FDR < 0.05 in the Group 3 rats compared to Group 1. Gene lists for each module can be found in Suppl. File 3 , and heat maps showing correlation between all modules and each sample can be found in Suppl. Fig. 11 . Fig. 4. Identification of Gene Modules Associated with High Fentanyl Escalation. Open in a new tab A Graphical representation of the approach of pyDESeq2 and PyWGCNA. B Heatmap displaying the module-trait relationship for each eigengene, color-coded for directionality of expression (downregulated – blue, upregulated – red). The first three rows indicate group-based model assignments of fentanyl escalation phenotypes. The top value in each box is R and the values in parentheses are P values, with significant values (P < 0.05) displayed in white. C Function enrichment analysis result tables of Brown module conducted via enrichr, ( C.i ) GO Biological Processes 2023, ( C.ii ) GO Cellular Component 2023, ( C.iii ) GO Molecular Function 2023. Term columns are the GO terms, Adjusted P-value are the corrected p-values, overlap is the number of genes in the Brown module that appear in the GO term, and genes is the list of genes from the Brown module in the GO term. D Function enrichment analysis bar plot of Brown module conducted via enrichr for the GWAS catalog. The y-axis are the GWAS terms, the x-axis is the -log10 of the adjusted p-value, and the parentheses to the right of the bar plot are the overlap of genes from the Brown module and genes in the GWAS catalog term. BART transcriptional regulator predictions generated via the BART mouse library ( E ) and the BART human library ( F ). In E and F red dots are the top 5 predicted transcriptional regulators. The y-axis represents the log adjusted p-value and the x-axis represents the BART-associated Z-score. We next performed functional enrichment analysis of the three modules significantly correlated to fentanyl escalation using the Biological Process, Cellular Component, and Molecular Function ontologies [ 38 , 39 ], as well as the GWAS catalogue [ 40 ]. For the Brown module, there was a total of 46 statistically significant terms across all three ontologies (Biological Process, Cellular Component, and Molecular Function). Terms containing potassium channel genes were present in all three ontologies and terms associated with glutamatergic signaling were found in the Biological Process and Cellular Component ontologies (Fig. 4C ). The Whitesmoke Module had a total of 17 statistically significant terms across all three ontologies and included potassium channel genes, but not glutamatergic signaling genes (Suppl. File 4 ). In the Lightgrey module, a single Biological Process term (Lipid Oxidation (GO:0034440)) was significant, and no significant terms were identified in the Cellular Component or Molecular Function ontologies (Suppl. File 4 ). Finally, we explored if any of the modules were enriched for genes that have SNPs associated with substance use disorders (SUDs) via the GWAS catalog accessed via the GSEAPY enrichr function. The Brown module had 157 GWAS catalogue terms that were statistically significant, and 3 of the top 5 GWAS catalog terms were strongly associated with SUDs: Smoking Initiation , Schizophrenia , and Drinks Per Week (Fig. 4D ). While the total number of statistically significant GWAS catalog terms for the Brown module is rather large, this was not the case for all modules, as the Whitesmoke module only had 10 GWAS catalogue terms that were statistically significant, the top term being schizophrenia (Suppl. File 4 ). The Lightgrey module had 5 significant GWAS catalogue terms that were mostly related to red blood cell characteristics (Suppl. File 4 ). Furthermore, of the 8 modules identified via the WGCNA that were not correlated to fentanyl escalation group, only one (Black) had any statistically significant GWAS catalogue terms [ 47 ], many of which were associated with immune function (Suppl. File 4 ). Transcriptional drivers of fentanyl escalation To further explore potential mechanisms driving phenotypic differences in fentanyl escalation, we used Binding Analysis for Regulation of Transcription (BART) to predict transcription factors regulating each module based on published mouse and human chromatin immunoprecipitation-sequencing (ChIP-Seq) libraries [ 41 ], since rat libraries are not yet available. For the Brown module, two of the top 5 predicted transcription factors overlapped in mouse and human libraries: enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2) [mouse: absolute rank = 1, p < 0.001; human: absolute rank = 2, p < 0.001] and jumonji and AT-rich interaction domain containing 2 (JARID2) [mouse: absolute rank = 4, p = 0.001; human: absolute rank = 4, p < 0.001] (Fig. 4E, F ). EZH2 is a core subunit of the Polycomb repressive complex 2 (PRC2) and is prominently involved in cancer [ 47 , 48 ]. Consistently, in the Whitesmoke module, EZH2 was also in the top 5 predicted transcription factors [mouse: absolute rank = 5, p = 0.002; human: absolute rank = 5, p < 0.001], while JARID2 was in the top 15 [mouse: absolute rank = 14, p = 0.001; human: absolute rank = 11, p = 0.003] (Suppl. File 5 ). Neither transcription factor reached significance in the Lightgrey module in results from the human or mouse libraries. Modeling contributions of escalation risk genes in biological replicates Finally, we conducted RT-PCR on a select number of genes identified in the Brown module in biological replicates of rats from Groups 1, 2, and 3. Given that the Brown module was positively correlated with Group 3 and negatively correlated with Group 1, we combined the differential expression data of Group 3 versus Group 1 with the Brown module to generate a priority scores for each gene within the Brown module (Suppl. File 6 ). We selected 4 potassium channels ( Kcna2, Kcnma1, Kcnh7, Kcnq3 ) in the top 25 genes given that function enrichment analysis of the Brown module highlighted terms associated with potassium channel activity and prior evidence implicating potassium channel gene dysregulation in SUDs [ 49 ]. We additionally investigated 3 genes in the top 25 associated with glutamatergic signaling ( Grin2a , Grin2b , Slc1a2 ) given the highlighted terms associated with glutamatergic signaling, and 2 genes in the top 25 associated with GABAergic interneurons ( Erbb4 , Tnr ) (Fig. 5A–I ). Fig. 5. Testing of Escalation Risk Genes in Biological Replicates. Open in a new tab A – I Bar graphs show qPCR analysis of select potassium channels, glutamatergic signaling, and GABAergic interneuron associated genes in the three escalation groups. J Linear regression of predicted versus observed total z-score value. Predicted values were generated using gene expression values and combined lever press z-score calculated from the last 7 days of self-administration. Green circles represent Group 1, blue squares represent Group 2, and red triangles represent Group 3. K Coefficients used in the gene expression model to predict the combined z-score values. Kcna2 , Potassium Voltage-Gated Channel Subfamily A Member 2; Kcnq3 , Potassium Voltage-Gated Channel Subfamily Q Member 3. [n = 4 to 7 per group; Fig. 5A–I are mean ± SEM]. Similar to the conclusions of permutations-based differential expression analyses, no individual gene in the PCR data reached statistical significance in a one-way ANOVA treating fentanyl escalation group as the factor. However, linear regression analyses found that expression levels of Kcnq3 and Kcna2 predicted the combined lever press z -score (R 2 = 0.432; F(2,13) = 4.942, p = 0.025, see Methods). The predicted z- score values from the last seven days of self-administration significantly correlated to the observed values (R 2 = 0.4319, F(1,14) = 10.64, p = 0.006) (Fig. 5J, K ). Similar to the z-score approach, linear regression analyses using qPCR data and PC1 of each rat also found expression of Kcnq3 and Kcna2 having statistically significant predictive validity of respective PC1 values (R 2 = 0.435; F(2,13) = 5.014, p = 0.024; Suppl. Fig. 12 ). Discussion A hallmark characteristic of OUD in the DSM-5 emphasizes a user’s overall “loss of control”, which leads to the compulsion to seek and obtain drugs, ultimately resulting in increased consumption [ 50 ]. Preclinical rodent models of escalated drug intake allow rats extended access to fentanyl in SA sessions, resulting in heterogeneity of intake across rats, and individual differences in fentanyl intake may be related to inherent variability in reward sensitivity. First, using latent growth curve models, we found that sucrose breakpoints were not predictive of fentanyl acquisition nor change during the escalation phase, although initial reinforced lever pressing in the escalation phase was predicted by the antecedent sucrose breakpoints. Group-based trajectory modeling was then used to isolate four non-anomalous distinct behavioral phenotypes in the fentanyl acquisition phase and three non-anomalous behavioral phenotypes in the escalation phase. The three-group model following fentanyl escalation was used for transcriptomic analyses for the remainder of the study, and bulk RNA-sequencing of the PFC was utilized to determine differentially expressed genes among the three groups. While a conventional monogenic approach identified a number of genes associated with the endophenotype displaying high escalation of fentanyl intake, the permutation-based monogenic approach revealed no significant associations. Given that complex behavioral responses, such as fentanyl seeking, are likely to develop on the background of cellular plasticity originating from multiple gene loci, a polygenic approach using WGCNA and functional enrichment analysis revealed several transcriptomic targets involving voltage-gated potassium channels, and excitatory glutamate signaling. Further, GWAS catalog functional enrichment results identified modules that were enriched for genes associated with known risk alleles for substance use disorders. Taken together, identification of gene networks that examine the joint contribution of genes may provide insight into druggable targets that underlie the complex endophenotype associated with the ‘loss of control’ phenomenon seen in OUD. A major finding in the current study was that breakpoints in the sucrose PR task did not predict subsequent fentanyl self-administration in acquisition sessions or the temporal change in lever presses over the 21-day escalation phase. It is possible that short duration fentanyl sessions limited the ability of sucrose PR to predict fentanyl seeking during the acquisition phase. That is, while each individual rat may possess the phenotypic features that allow for an association between sucrose PR and fentanyl self-administration, a short testing window precluded us from observing such an association. Although there are reports of chronic opioid exposure leading to increased sugar intake in both preclinical and clinical models (For review, see [ 51 ]) previous assessments of whether responsivity to sucrose can predict subsequent opioid self-administration is scarce. Vasquez, Shah, Re, Laezza and Green [ 13 ] found that PR breakpoints for sucrose were not predictive of PR breakpoints for fentanyl SA, although a frustrative non-reward task was able to predict PR breakpoints for fentanyl SA. Within opioids, assessments of antecedent sensitivity to non-drug reward have been evaluated using several different paradigms, with mixed results. A study by McNamara, Dalley, Robbins, Everitt and Belin [ 12 ] was not able to predict escalation in heroin SA with the 5-choice serial reaction time (5CSRT), a task that measures attention and impulsive action in rodents via the ability to withhold a behavioral response to a reward-associated stimulus [ 52 ]. Contrary to the current findings with fentanyl, unpublished data from our laboratory has shown that antecedent sensitivity to non-drug reward, measured by an identical PR task for sucrose, did predict responding for cocaine SA in a PR task (personal communication, Dr. Pavel Ortinski). Similar results showing that discrete traits predict drug intake have been reported when assessing stimulants, as individual differences in baseline impulsive action on the 5CSRT task have predicted greater rates of SA with nicotine [ 53 ] and cocaine [ 54 ]. Moreover, a study in female rats used a delay discounting task for food reinforcement to screen for high or low impulsivity and found that high impulsivity rats escalated cocaine intake, while low impulsivity rats did not [ 14 ]. Given the mixed results and pharmacological differences between stimulants and opioids, the inability of non-drug reward sensitivity to predict fentanyl escalation may lend support to the hypothesis that chronic opioid SA itself may cause effects on “top-down” control within the PFC [ 19 ]. Moreover, executive control over drug taking becomes progressively weaker as drug experience increases, leading to the “loss of control” phenomenon [ 55 ]. Whereas the present study failed to identify antecedent behavioral markers of risk, the transcriptomic response of opioid intake may provide insight into the identification of biomarkers that govern this PFC-mediated process of escalation of drug intake. While the phenotypic identification of intake trajectories is an artifact of overall operant behavior, there may be genetic or epigenetic differences that influence the emergence of group membership. The failure of a permutation-based monogenic approach to produce candidate genes may have occurred because the ‘loss of control’ phenotype seen in OUD relies on several subtle changes within the transcriptomic landscape triggered by escalated fentanyl use. While GWAS of OUD has generated multiple candidate loci associated with increased risk of OUD ( OPRM1 , DRD2 , KCNG2 , KCNC1 , CNIH3 , etc.) [ 56 , 57 ], previous research has demonstrated that many psychiatric traits including SUDs, initiation of drug use, and frequency of drug use are highly polygenic [ 58 ]. The use of polygenetic risk scores (PRS), which can be used to estimate the risk of disease state development by aggregating the risk of multiple genetic variations from a discovery GWAS [ 59 ], have highlighted the cumulative nature of genetic variation and risk for OUD development. A large GWAS study including ~640,000 individuals produced an OUD PRS that was a significant predictor of OUD diagnosis [ 60 ]. This OUD PRS was also found to interact with education level and post-traumatic stress disorder, as both variables moderate OUD PRS risk for the development of OUD [ 61 ]. Further PRS for other psychiatric traits have been found to be associated with OUD diagnosis, as high PRS for risk-taking and neuroticism were associated with OUD diagnosis [ 62 ]. Thus, a polygenic approach using WGCNA was implemented to examine the joint contribution of genes to identify correlation patterns within the three distinct phenotypes ascertained by our model. Within the modules (ME) identified by WGCNA, the correlation pattern within the Brown ME is of particular interest. The “addiction-prone” high escalators in Group 3 were positively correlated to the Brown ME (0.54), while the low-trajectory intake Group 1 rats were negatively correlated (−0.54), and Group 2 showed virtually no directionality (−0.01). While it is possible that differing levels of overall fentanyl exposure could be contributing to the observed Group 1 and Group 3 transcriptomic differences, the Brown ME showed no correlation to Group 2 (R = −0.01, P = 0.98) despite Groups 2 and 3 having similar exposure levels (ie, similar rates of reinforced lever presses by the end of extended access SA sessions). This suggests that the observed transcriptomic correlations are not solely due to differences in fentanyl exposure, but a core difference in phenotype. However, further study is required to completely rule out the cumulative effect of increased fentanyl exposure across the 21 extended access SA sessions. Interestingly, the Brown ME contained Oprm1 , the gene for the μ-opioid receptor. One of the most consistently found SNPs associated with an increased likelihood of developing OUD is rs1799971 in OPRM1 , with both GWAS and candidate gene studies finding the minor G allele increasing risk of OUD, but these findings are not consistent and often no association is observed in people with non-European ancestry (for review see [ 56 ]). Additionally, it is worth noting that the majority of GWAS studies use non-OUD diagnosis as control, and in GWAS studies that compared OUD diagnosis to opioid exposed non-OUD controls the G allele was associated with decreased likelihood of OUD [ 63 ] or no association with OUD [ 62 ], albeit the later sample size was much smaller than the former. People with the AG/GG genotypes rated low doses of opioid more rewarding and experienced less unpleasant responses to opioid than AA genotypes [ 64 ], however this is not a consistent finding and differences in experimental design and sample size may bias results [ 65 ]. Genetic mouse models with the equivalent SNP (A112G) have been generated and used to investigate OUD-like phenotypes. G/G mutant mice had higher heroin and oxycodone self-administration compared to A/A controls [ 66 , 67 ]. Other SNPs in the OPRM1 gene have been associated with increased risk of OUD in both GWAS and candidate gene approaches. For instances multiple SNPs in OPRM1 were associated with altered mRNA expression in the PFC of postmortem humans, and some of these SNPs were associated with increased likelihood of heroin addiction [ 68 ]. Further, additional SNPs in OPRM1 have been associated with altered response to opioid. Individuals with a minor C allele in the OPRM1 SNP rs6848893 had decreased peak liking rating compared to major T allele carriers. However, the authors reported no impact of OPRM1 SNP rs1799971 on subjective effects of opioids [ 65 ], although sample size was small. While the inclusion of Oprm1 in the Brown ME is noteworthy, connectivity of Oprm1 was 7.22, which is the 194 th most interconnected gene out of the 335. The top 5 hub genes in the Brown ME based on connectivity are Pclo, Grin2a, Ppp1r9a, Slc9a7 , and Lrrc7 . Three of these genes are related to synapse structure and function: Pclo, Ppp1r9a , and Lrrc7 . Pclo is involved in presynaptic structure and function and contributes to calcium channel localization and is involved in synaptic vesicle maintenance [ 69 , 70 ]. Lrrc7 encodes protein Densin-180, which performs a similar synaptic structure and organization function postsynaptically, but has no role in synaptic vesicles [ 71 ]. Ppp1r9a encodes the protein Neurabin-I, which in neurons is associated with filopodia and dendritic spine formation [ 72 , 73 ]. These structural genes can influence strength of synaptic connections through recruitment and stabilization of receptors and/or proteins like CaMKII [ 74 ]. Interestingly the second most connected gene in the Brown ME was Grin2a , an NMDA receptor subunit. Further, Grin2b is also included in the Brown ME and both genes were included in the 10 differentially expressed genes in the Group 3 vs Group 1 comparison. In human postmortem tissue, GluN2B (encoded by GRIN2B ) expression was increased in the PFC of opioid users compared to non-drug user controls [ 75 ]. Injection of an NMDA receptor inhibitor into the prelimbic cortex of rats increased morphine induced conditioned place preference (CPP) at a sub rewarding morphine dose, but injections into the anterior cingulate cortex and infralimbic cortex had no impact on CPP. Further, a GluN2B specific antagonist injected into the prelimbic cortex also elicited morphine CPP at the sub rewarding morphine dose [ 76 ]. While this may appear contrary to our results, this increased CPP was likely attributable to increased rewarding effects of morphine as the NMDA antagonism elicited CPP at a sub rewarding morphine dose, and other studies have found antagonism of NDMA receptors in the prelimbic cortex reduced morphine self-administration [ 77 ]. Potentially structural ( Pclo, Ppp1r9a , and Lrrc7) and functional ( Grin2a and Grin2b ) differences at synapses in the PFC underlie differences in fentanyl reward. Increased NMDA activity, driven by both structural connectivity between synapses and increased NMDA subunits, may result in lower sensitivity to opioid reward which drives the higher responding rate in self-administration. Further experiments are needed to explore the potential role of these hub genes in fentanyl self-administration, and whether they act as concurrent up stream modulators of the Brown ME . While Pclo, Ppp1r9a , and Lrrc7 were not chosen for PCR follow-up, it is worth noting Grin2a and Grin2b were not associated with escalation phenotype in the biological replicate PCR analysis. A functional enrichment analysis of the Brown ME then identified several gene ontologies, including Voltage-Gated Monoatomic Cation Channel Activity (GO:0022843) and Voltage-Gated Potassium Channel Activity (GO:0005249), both of which contained genes ( Kcna2, Kcna3, Cacna1e) that were significantly upregulated when comparing the entire transcriptome of Group 3 and Group 1 using DESeq2 . The prominence of voltage-gated potassium channels genes in our analyses is consistent with their central role in regulation of neuronal firing [ 78 ]. Potassium channels have been linked previously to opioids [ 79 ], as the KCNQ (Kv7) activators retigabine and flupirtine, show antinociceptive effects [ 80 ], as well as enhanced analgesic effects when combined with morphine [ 81 , 82 ]. Moreover, whole-cell patch clamp electrophysiology has shown that fentanyl blocks hERG channels in a voltage dependent manner [ 83 ], an effect also seen with heroin and methadone [ 84 ]. Interestingly, we have previously explored non-opioid effects of potassium channels in high and low sucrose preferring rats. Sequencing data identified 6 differentially expressed potassium channel transcripts in the NAc between high and low sucrose preferring rats, and pharmacological manipulation of the K V 1.4 current increased responding in a sucrose PR task in low rats [ 22 ]. Transcription factor analysis via BART identified two transcription factors as potential mediators of our transcriptional changes, EZH2 and JARID2. EZH2 is a core subunit of the Polycomb repressive complex 2 (PRC2), which is a complex that regulates methylation of histone H3 lysine 2, with core units of EZH2, EED, SUZ12, RBBP7/4, and acts as a major transcriptional repressor [ 47 ]. Given the integral role of EZH2 in PRC2 histone methylation, it serves multiple important biological functions such as embryonic stem cell differentiation [ 85 – 87 ], nervous system development [ 88 , 89 ], immune cell development [ 90 ], skeletal formation [ 91 ], and axon regeneration in the peripheral nervous system [ 92 ]. In adult root ganglia cells, overexpression of EZH2 in vivo causes increased axon regeneration following injury and suppression of ion channel pathways. Interestingly, overexpression of a modified EZH2 with loss of methylation function, ablating the PRC2’s ability to suppress gene expression via HK27me3, causes an upregulation of ion channels pathways, including three differentially expressed genes that overlap with our Brown ME ( Kcnma1, Kcnq3, Grin2b ) [ 92 ]. This suggests a potential shift from EZH2-mediated histone methylation to a non-methylating role of EZH2 in modulating the observed upregulation of ion channel in Group 3 high escalating rats. EZH2 has multiple other non-canonical mechanisms to regulate transcription including complexing with or increasing expression of other transcriptional regulators [ 48 ] and a previous ATAC-seq study examining the striatum of heroin users found enriched enhancer regions of EZH2 [ 93 ]. It is worth mentioning that many of the studies demonstrating co-activation of transcription factors were conducted in cancer cells, as EZH2 is highly implicated in cancer [ 48 , 94 ]. The PRC2 can complex with multiple subunits that influence function, such as JARID2 [ 95 ], which influences embryonic development [ 96 , 97 ] and neurogenesis [ 98 ], and similar to EZH2, is highly implicated in cancer progression [ 99 ]. JARID2 mutations are associated with neurodevelopmental syndrome [ 100 ], and SNPs in JARID2 are associated with affective dysfunction during alcohol withdrawal [ 101 ]. GWAS can provide key insights into interpreting association from a transcriptional data set in a biological and genomic context [ 102 ]. Several significant GWAS terms were generated from the Brown ME , which was positively correlated with the “addiction prone” high escalating rats in Group 3, that may be associated with SUD such as Smoking Initiation, Schizophrenia , and Drinks Per Week . Gelernter, Kranzler, Sherva, Koesterer, Almasy, Zhao and Farrer [ 49 ] conducted a GWAS with 5,697 individuals from the United States with OUD and found the most significant genes were involved in potassium and calcium signaling pathways. Of note, the voltage-gated potassium channels, KCNC1 and KCNG2 , had single-nucleotide polymorphisms (SNPs) that were significantly associated with OUD. Additionally, SNPs in CACNA1C , a voltage-gated calcium channel, were associated with risk for schizophrenia [ 103 ] and bipolar affective disorder [ 104 ], two psychiatric conditions often co-morbid with SUD. Taken together with the transcription factor analyses, this suggests that a narrow group of genes in the Brown ME may be driving susceptibility to escalated fentanyl intake. This supposition is further supported by our PCR data. While the permutation-based monogenic approach identified no single gene that was statistically different between the three fentanyl escalation groups, two K + channels, Kcnq3 and Kcna2 , were found to predict the combined lever press z- score over the last 7 self-administration sessions. While the literature is scant on either Kcnq3 or Kcna2 role in opioid escalation, contrasting to our findings, previous work found that multiple potassium channel mRNAs were down regulated in the PFC of high oxycodone self-administering rats compared to the low self-administering group. However, Kcna2 was not explored, and while Kcnq3 expression was reduced in the high rats compared to the low this did not reach statistical significance [ 105 ]. In the case of Kcna2 specifically, multiple rodent models of neuropathic pain suggest a potential role of Kcna2 [ 106 ], highlighting a gene that may carry risk alleles for OUD development in patients prescribed opioids for neuropathic pain. This supposition is purely speculative, and further investigation is needed for a stronger association between genetic variation in Kcna2 and OUD risk. In conclusion, we sought to determine whether non-drug reward sensitivity could predict fentanyl escalation, and to determine the transcriptomic signatures within the PFC that are associated with the “loss of control” phenomenon seen in OUD. The gene networks associated with fentanyl escalation highlighted in the present study may represent concordant genomic changes essential to the “loss of control” endophenotype present in OUD. Furthermore, exploratory studies assessing the potential of these gene expression changes as drug discovery targets for OUD, namely CNS excitation networks via voltage-gated potassium channels may lead to new classes of therapeutics for substance use disorders. Supplementary information Supplemental File 1 (1.3MB, pdf) Supplemental File 2 (3.7MB, xlsx) Supplemental File 3 (577.5KB, xlsx) Supplemental File 4 (4.7MB, xlsx) Supplemental File 5 (888.6KB, xlsx) Supplemental File 6 (42.7KB, xlsx) Acknowledgements This work is supported by: NIH R01 DA053070; U01 DA051377; F31DA057812; the University of Kentucky SUPRA Faculty Pilot Grant; T32 DA035200; AFPE regional award; the University of Kentucky SUPRA Graduate/Professional Student Pilot Grant; AFPE Pre-Doctoral Fellowship. RNA-sequencing data generated in this study have been deposited in the NCBI Gene Expression Omnibus (GEO) under accession number GSE311548 . The submission includes raw read files, processed count tables, and complete metadata for all samples. Data analyses were conducted using a combination of publicly available software packages (SAS, R, MATLAB, Galaxy, pyDESeq2, PyWGCNA, GSEAPY, and BART; versions specified in Methods). Analysis workflows and parameters are described in the Methods and Supplementary Materials. Custom scripts used for data processing are available from the corresponding author upon reasonable request. This manuscript is currently available as a preprint on bioRxiv. 10.1101/2024.12.06.627259. Author contributions JK performed bioinformatic data analysis, manuscript writing, and manuscript editing. RC designed experiments, performed behavioral and bioinformatic data analysis, and helped write and edit manuscript. JS performed biochemical assays, helped write and edit manuscript. EP performed biochemical assays, manuscript writing, and manuscript editing. MX performed animal surgeries and ran animal behavioral assays. ED performed animal surgeries and ran animal behavioral assays. CB assisted in bioinformatics analysis and edited the manuscript. JM assisted in bioinformatics analysis and edited the manuscript. PO procured funding, designed experiments, and helped write and edit manuscript. MB procured funding, designed experiments, and helped write and edit manuscript. JT procured funding, designed experiments, and helped write and edit manuscript. Competing interests All authors report no biomedical financial interests or potential conflicts of interest. Ethics and approval All experimental procedures were followed in accordance with the NIH Guide for the Care and Use of Laboratory Animals (8th edition, 2011) and approved by the University of Kentucky Institutional Animal Care and Use Committee (protocol # 20221-3779). Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. These authors contributed equally: Jack Keady, Richard Charnigo, Jakob D. Shaykin. Supplementary information The online version contains supplementary material available at 10.1038/s41398-026-04034-1. References 1. Vermeire A, Remon JP, Rosseel MT, Belpaire F, Devulder J, Bogaert MG. Variability of morphine disposition during long-term subcutaneous infusion in terminally ill cancer patients. Eur J Clin Pharmacol. 1998;53:325–30. [ DOI ] [ PubMed ] [ Google Scholar ] 2. Chen ZR, Somogyi AA, Reynolds G, Bochner F. Disposition and metabolism of codeine after single and chronic doses in one poor and seven extensive metabolisers. Br J Clin Pharmacol. 1991;31:381–90. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 3. Säwe J. High-dose morphine and methadone in cancer patients. Clinical pharmacokinetic considerations of oral treatment. Clin Pharmacokinet. 1986;11:87–106. [ DOI ] [ PubMed ] [ Google Scholar ] 4. Hayes CJ, Krebs EE, Hudson T, Brown J, Li C, Martin BC. Impact of opioid dose escalation on the development of substance use disorders, accidents, self-inflicted injuries, opioid overdoses and alcohol and non-opioid drug-related overdoses: a retrospective cohort study. Addiction. 2020;115:1098–112. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 5. Vendruscolo LF, Schlosburg JE, Misra KK, Chen SA, Greenwell TN, Koob GF. Escalation patterns of varying periods of heroin access. Pharmacol Biochem Behav. 2011;98:570–4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 6. Wade CL, Vendruscolo LF, Schlosburg JE, Hernandez DO, Koob GF. Compulsive-like responding for opioid analgesics in rats with extended access. Neuropsychopharmacology. 2015;40:421–8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 7. Harding WM, Zinberg NE. Occasional opiate use. Adv Subst Abuse. 1983;3:27–61. [ Google Scholar ] 8. Lenoir M, Ahmed SH. Heroin-Induced reinstatement is specific to compulsive heroin use and dissociable from heroin reward and sensitization. Neuropsychopharmacology. 2007;32:616–24. [ DOI ] [ PubMed ] [ Google Scholar ] 9. Malone SG, Keller PS, Hammerslag LR, Bardo MT. Escalation and reinstatement of fentanyl self-administration in male and female rats. Psychopharmacology (Berl). 2021;238:2261–73. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 10. Gipson CD, Dunn KE, Bull A, Ulangkaya H, Hossain A. Establishing preclinical withdrawal syndrome symptomatology following heroin self-administration in male and female rats. Exp Clin Psychopharmacol. 2021;29:636–49. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Zernig G, Ahmed SH, Cardinal RN, Morgan D, Acquas E, Foltin RW, et al. Explaining the escalation of drug use in substance dependence: models and appropriate animal laboratory tests. Pharmacology. 2007;80:65–119. [ DOI ] [ PubMed ] [ Google Scholar ] 12. McNamara R, Dalley JW, Robbins TW, Everitt BJ, Belin D. Trait-like impulsivity does not predict escalation of heroin self-administration in the rat. Psychopharmacology (Berl). 2010;212:453–64. [ DOI ] [ PubMed ] [ Google Scholar ] 13. Vasquez TES, Shah P, Re JD, Laezza F, Green TA. Individual differences in frustrative nonreward behavior for sucrose in rats predict motivation for fentanyl under progressive ratio. eNeuro. 2021;8:ENEURO.0136-21.2021. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 14. Anker JJ, Perry JL, Gliddon LA, Carroll ME. Impulsivity predicts the escalation of cocaine self-administration in rats. Pharmacol Biochem Behav. 2009;93:343–8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 15. Perry JL, Larson EB, German JP, Madden GJ, Carroll ME. Impulsivity (delay discounting) as a predictor of acquisition of IV cocaine self-administration in female rats. Psychopharmacology (Berl). 2005;178:193–201. [ DOI ] [ PubMed ] [ Google Scholar ] 16. Swain Y, Waller NG, Gewirtz JC, Harris AC. Individual differences in different measures of opioid self-administration in rats are accounted for by a single latent variable. Front Psychiatry. 2021;12:712163. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 17. Perry JL, Joseph JE, Jiang Y, Zimmerman RS, Kelly TH, Darna M, et al. Prefrontal cortex and drug abuse vulnerability: Translation to prevention and treatment interventions. Brain Res Rev. 2011;65:124–49. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 18. Anastasiades PG, Carter AG. Circuit organization of the rodent medial prefrontal cortex. Trends Neurosci. 2021;44:550–63. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 19. Kalivas PW, Volkow ND. The neural basis of addiction: a pathology of motivation and choice. Am J Psychiatry. 2005;162:1403–13. [ DOI ] [ PubMed ] [ Google Scholar ] 20. Liu SX, Gades MS, Swain Y, Ramakrishnan A, Harris AC, Tran PV, et al. Repeated morphine exposure activates synaptogenesis and other neuroplasticity-related gene networks in the dorsomedial prefrontal cortex of male and female rats. Drug Alcohol Depend. 2021;221:108598. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 21. Browne CJ, Godino A, Salery M, Nestler EJ. Epigenetic Mechanisms of Opioid Addiction. Biol Psychiatry. 2020;87:22–33. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 22. O’Donovan B, Adeluyi A, Anderson EL, Cole RD, Turner JR, Ortinski PI. Altered gating of Kv1.4 in the nucleus accumbens suppresses motivation for reward. eLife. 2019;8:e47870. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 23. Shaykin JD, Olyha LN, Van Doorn CE, Hales JD, Chandler CM, Hopkins DM, et al. Effects of isolation stress and voluntary ethanol exposure during adolescence on ethanol and nicotine co-use in adulthood using male rats. Drug Alcohol Depend Rep. 2024;12:100277. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 24. Testen A, Sepulveda-Orengo MT, Gaines CH, Reissner KJ. Region-Specific reductions in morphometric properties and synaptic colocalization of astrocytes following cocaine self-administration and extinction. Front Cell Neurosci. 2018;12:246. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Bassareo V, Cucca F, Frau R, Di Chiara G. Changes in dopamine transmission in the nucleus accumbens shell and core during ethanol and sucrose self-administration. Front Behav Neurosci. 2017;11:71. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Richardson NR, Roberts DC. Progressive ratio schedules in drug self-administration studies in rats: a method to evaluate reinforcing efficacy. J Neurosci Methods. 1996;66:1–11. [ DOI ] [ PubMed ] [ Google Scholar ] 27. Jones BL, Nagin DS, Roeder K. A SAS procedure based on mixture models for estimating developmental trajectories. Sociol Methods Res. 2001;29:374–93. [ Google Scholar ] 28. Jones BL, Nagin DS. Advances in group-based trajectory modeling and an SAS procedure for estimating them. Sociol Methods Res. 2007;35:542–71. [ Google Scholar ] 29. Ligges UM. Martin Scatterplot3d - an R package for visualizing multivariate data. J Stat Softw. 2003;8:1–20. [ Google Scholar ] 30. Batut B, van den Beek M, Doyle MA, Soranzo N. RNA-Seq data analysis in galaxy. Methods Mol Biol. 2021;2284:367–92. [ DOI ] [ PubMed ] [ Google Scholar ] 31. Hiltemann S, Rasche H, Gladman S, Hotz HR, Lariviere D, Blankenberg D, et al. Galaxy Training: A powerful framework for teaching!. PLoS Comput Biol. 2023;19:e1010752. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 32. Batut B, Hiltemann S, Bagnacani A, Baker D, Bhardwaj V, Blank C, et al. Community-driven data analysis training for biology. Cell Syst. 2018;6:752–8.e751. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 33. Batut MFB, Heydarian M, Erxleben A, Videm P, Blank B, Doyle M, et al. Reference-based RNA-Seq data analysis (Galaxy Training Materials). https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/ref-based/tutorial.html Online; accessed Mon Jan 14 2024. 34. Muzellec B, Telenczuk M, Cabeli V, Andreux M. PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics. 2023;39:btad547. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 35. Rezaie N, Reese F, Mortazavi A. PyWGCNA: a Python package for weighted gene co-expression network analysis. Bioinformatics. 2023;39:btad415. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 36. Fang Z, Liu X, Peltz G. GSEApy: a comprehensive package for performing gene set enrichment analysis in Python. Bioinformatics. 2023;39:btac757. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 37. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44:W90–7. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 38. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 39. C. Gene Ontology, Aleksander SA, Balhoff J, Carbon S, Cherry JM, Drabkin HJ, et al. The Gene Ontology knowledgebase in 2023. Genetics. 2023;224:iyad031. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 40. Sollis E, Mosaku A, Abid A, Buniello A, Cerezo M, Gil L, et al. The NHGRI-EBI GWAS Catalog: knowledgebase and deposition resource. Nucleic Acids Res. 2023;51:D977–D85. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 41. Wang Z, Civelek M, Miller CL, Sheffield NC, Guertin MJ, Zang C. BART: a transcription factor prediction tool with query gene sets or epigenomic profiles. Bioinformatics. 2018;34:2867–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 42. Kumar M, Rainville JR, Williams K, Lile JA, Hodes GE, Vassoler FM, et al. Sexually dimorphic neuroimmune response to chronic opioid treatment and withdrawal. Neuropharmacology. 2021;186:108469. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 43. Adeluyi A, Guerin L, Fisher ML, Galloway A, Cole RD, Chan SSL, et al. Microglia morphology and proinflammatory signaling in the nucleus accumbens during nicotine withdrawal. Sci Adv. 2019;5:eaax7031. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 44. Fisher ML, LeMalefant RM, Zhou L, Huang G, Turner JR. Distinct roles of CREB within the ventral and dorsal hippocampus in mediating nicotine withdrawal phenotypes. Neuropsychopharmacology. 2017;42:1599–609. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 45. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25:402–8. [ DOI ] [ PubMed ] [ Google Scholar ] 46. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 47. Margueron R, Reinberg D. The Polycomb complex PRC2 and its mark in life. Nature. 2011;469:343–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 48. Wang J, Wang GG. No Easy Way Out for EZH2: its pleiotropic, noncanonical effects on gene regulation and cellular function. Int J Mol Sci. 2020;21:9501. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 49. Gelernter J, Kranzler HR, Sherva R, Koesterer R, Almasy L, Zhao H, et al. Genome-Wide association study of opioid dependence: multiple associations mapped to calcium and potassium pathways. Biol Psychiatry. 2014;76:66–74. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 50. Mathis V, Kenny PJ. From controlled to compulsive drug-taking: The role of the habenula in addiction. Neurosci Biobehav Rev. 2019;106:102–11. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 51. Mysels DJ, Sullivan MA. The relationship between opioid and sugar intake: review of evidence and clinical applications. J Opioid Manag. 2010;6:445–52. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 52. SK Asinof, TA Paine. The 5-choice serial reaction time task: a task of attention and impulse control for rodents. J Vis Exp. 2014;90:e51574. [ DOI ] [ PMC free article ] [ PubMed ] 53. Diergaarde L, Pattij T, Poortvliet I, Hogenboom F, de Vries W, Schoffelmeer AN, et al. Impulsive choice and impulsive action predict vulnerability to distinct stages of nicotine seeking in rats. Biol Psychiatry. 2008;63:301–8. [ DOI ] [ PubMed ] [ Google Scholar ] 54. Dalley JW, Fryer TD, Brichard L, Robinson ES, Theobald DE, Lääne K, et al. Nucleus accumbens D2/3 receptors predict trait impulsivity and cocaine reinforcement. Science. 2007;315:1267–70. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 55. Everitt BJ, Robbins TW. Neural systems of reinforcement for drug addiction: from actions to habits to compulsion. Nat Neurosci. 2005;8:1481–9. [ DOI ] [ PubMed ] [ Google Scholar ] 56. Crist RC, Reiner BC, Berrettini WH. A review of opioid addiction genetics. Curr Opin Psychol. 2019;27:31–5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 57. Deak JD, Johnson EC. Genetics of substance use disorders: a review. Psychol Med. 2021;51:2189–200. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 58. Watanabe K, Stringer S, Frei O, Umićević Mirkov M, de Leeuw C, Polderman TJC, et al. A global overview of pleiotropy and genetic architecture in complex traits. Nat Genet. 2019;51:1339–48. [ DOI ] [ PubMed ] [ Google Scholar ] 59. Kember RL, Davis CN, Feuer KL, Kranzler HR. Considerations for the application of polygenic scores to clinical care of individuals with substance use disorders. J Clin Invest. 2024;134:e172882. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 60. Deak JD, Zhou H, Galimberti M, Levey DF, Wendt FR, Sanchez-Roige S, et al. Genome-wide association study in individuals of European and African ancestry and multi-trait analysis of opioid use disorder identifies 19 independent genome-wide significant risk loci. Mol Psychiatry. 2022;27:3970–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 61. Na PJ, Deak JD, Kranzler HR, Pietrzak RH, Gelernter J. Genetic and non-genetic predictors of risk for opioid dependence. Psychol Med. 2024;54:1779–86. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 62. Polimanti R, Walters RK, Johnson EC, McClintick JN, Adkins AE, Adkins DE, et al. Leveraging genome-wide data to investigate differences between opioid use vs. opioid dependence in 41,176 individuals from the Psychiatric Genomics Consortium. Mol Psychiatry. 2020;25:1673–87. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 63. Zhou H, Rentsch CT, Cheng Z, Kember RL, Nunez YZ, Sherva RM, et al. Association of OPRM1 functional coding variant with opioid use disorder: a genome-wide association study. JAMA Psychiatry. 2020;77:1072–80. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 64. Dunn KE, Huhn AS, Finan PH, Mange A, Bergeria CL, Maher BS, et al. Polymorphisms in the A118G SNP of the OPRM1 gene produce different experiences of opioids: A human laboratory phenotype-genotype assessment. Addict Biol. 2024;29:e13355. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 65. Jones JD, Mumtaz M, Manubay JM, Mogali S, Sherwin E, Martinez S, et al. Assessing the contribution of opioid- and dopamine-related genetic polymorphisms to the abuse liability of oxycodone. Pharmacol Biochem Behav. 2019;186:172778. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 66. Zhang Y, Picetti R, Butelman ER, Ho A, Blendy JA, Kreek MJ. Mouse model of the OPRM1 (A118G) polymorphism: differential heroin self-administration behavior compared with wild-type mice. Neuropsychopharmacology. 2015;40:1091–100. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 67. Zhang Y, Randesi M, Blendy JA, Kreek MJ, Butelman ER. Impact of OPRM1 (Mu-opioid Receptor Gene) A112G polymorphism on dual oxycodone and cocaine self-administration behavior in a mouse model. Neuroscience. 2024;539:76–85. [ DOI ] [ PubMed ] [ Google Scholar ] 68. Hancock DB, Levy JL, Gaddis NC, Glasheen C, Saccone NL, Page GP, et al. Cis-Expression quantitative trait loci mapping reveals replicable associations with heroin addiction in OPRM1. Biol Psychiatry. 2015;78:474–84. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 69. Gundelfinger ED, Reissner C, Garner CC. Role of bassoon and piccolo in assembly and molecular organization of the active zone. Front Synaptic Neurosci. 2015;7:19. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 70. Ackermann F, Schink KO, Bruns C, Izsvák Z, Hamra FK, Rosenmund C, et al. Critical role for Piccolo in synaptic vesicle retrieval. eLife. 2019;8:e46629. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 71. Thalhammer A, Trinidad JC, Burlingame AL, Schoepfer R. Densin-180: revised membrane topology, domain structure and phosphorylation status. J Neurochem. 2009;109:297–302. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 72. Terry-Lorenzo RT, Roadcap DW, Otsuka T, Blanpied TA, Zamorano PL, Garner CC, et al. Neurabin/protein phosphatase-1 complex regulates dendritic spine morphogenesis and maturation. Mol Biol Cell. 2005;16:2349–62. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 73. Ryan XP, Alldritt J, Svenningsson P, Allen PB, Wu GY, Nairn AC, et al. The Rho-specific GEF Lfc interacts with neurabin and spinophilin to regulate dendritic spine morphology. Neuron. 2005;47:85–100. [ DOI ] [ PubMed ] [ Google Scholar ] 74. Strack S, Robison AJ, Bass MA, Colbran RJ. Association of calcium/calmodulin-dependent kinase II with developmentally regulated splice variants of the postsynaptic density protein densin-180. J Biol Chem. 2000;275:25061–4. [ DOI ] [ PubMed ] [ Google Scholar ] 75. Daneshparvar H, Sadat-Shirazi MS, Fekri M, Khalifeh S, Ziaie A, Esfahanizadeh N, et al. NMDA receptor subunits change in the prefrontal cortex of pure-opioid and multi-drug abusers: a post-mortem study. Eur Arch Psychiatry Clin Neurosci. 2019;269:309–15. [ DOI ] [ PubMed ] [ Google Scholar ] 76. Bishop SF, Lauzon NM, Bechard M, Gholizadeh S, Laviolette SR. NMDA receptor hypofunction in the prelimbic cortex increases sensitivity to the rewarding properties of opiates via dopaminergic and amygdalar substrates. Cereb Cortex. 2011;21:68–80. [ DOI ] [ PubMed ] [ Google Scholar ] 77. Aboutalebi F, Alaei H, Oryan S, Radahmadi M. Blockade of prelimbic glutamate receptor reduces the reinforcing effect of morphine. Can J Physiol Pharmacol. 2018;96:815–22. [ DOI ] [ PubMed ] [ Google Scholar ] 78. Jan LY, Jan YN. Voltage-gated potassium channels and the diversity of electrical signalling. J Physiol. 2012;590:2591–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 79. McCoy MT, Jayanthi S, Cadet JL. Potassium channels and their potential roles in substance use disorders. Int J Mol Sci. 2021;22:1249. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 80. Dost R, Rostock A, Rundfeldt C. The anti-hyperalgesic activity of retigabine is mediated by KCNQ potassium channel activation. Naunyn Schmiedebergs Arch Pharmacol. 2004;369:382–90. [ DOI ] [ PubMed ] [ Google Scholar ] 81. Goodchild CS, Kolosov A, Tucker AP, Cooke I. Combination therapy with flupirtine and opioid: studies in rat pain models. Pain Med. 2008;9:928–38. [ DOI ] [ PubMed ] [ Google Scholar ] 82. Kolosov A, Goodchild CS, Williams ED, Cooke I. Flupirtine enhances the anti-hyperalgesic effects of morphine in a rat model of prostate bone metastasis. Pain Med. 2012;13:1444–56. [ DOI ] [ PubMed ] [ Google Scholar ] 83. Tschirhart JN, Li W, Guo J, Zhang S. Blockade of the human ether A-Go-Go–related gene (hERG) potassium channel by fentanyl. Mol Pharmacol. 2019;95:386–97. [ DOI ] [ PubMed ] [ Google Scholar ] 84. Zünkler BJ, Wos-Maganga M. Comparison of the effects of methadone and heroin on human ether-à-go-go-Related gene channels. Cardiovasc Toxicol. 2010;10:161–5. [ DOI ] [ PubMed ] [ Google Scholar ] 85. Shan Y, Liang Z, Xing Q, Zhang T, Wang B, Tian S, et al. PRC2 specifies ectoderm lineages and maintains pluripotency in primed but not naive ESCs. Nat Commun. 2017;8:672. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 86. Collinson A, Collier AJ, Morgan NP, Sienerth AR, Chandra T, Andrews S, et al. Deletion of the polycomb-group protein EZH2 leads to compromised self-renewal and differentiation defects in human embryonic stem cells. Cell Rep. 2016;17:2700–14. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 87. Illingworth RS, Holzenspies JJ, Roske FV, Bickmore WA, Brickman JM. Polycomb enables primitive endoderm lineage priming in embryonic stem cells. eLife. 2016;5:e14926. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 88. Zhang M, Zhang Y, Xu Q, Crawford J, Qian C, Wang GH, et al. Neuronal histone methyltransferase EZH2 regulates neuronal morphogenesis, synaptic plasticity, and cognitive behavior in mice. Neurosci Bull. 2023;39:1512–32. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 89. Kim H, Langohr IM, Faisal M, McNulty M, Thorn C, Kim J. Ablation of Ezh2 in neural crest cells leads to aberrant enteric nervous system development in mice. PLoS ONE. 2018;13:e0203391. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 90. Nutt SL, Keenan C, Chopin M, Allan RS. EZH2 function in immune cell development. Biol Chem. 2020;401:933–43. [ DOI ] [ PubMed ] [ Google Scholar ] 91. Dudakovic A, Camilleri ET, Xu F, Riester SM, McGee-Lawrence ME, Bradley EW, et al. Epigenetic control of skeletal development by the histone methyltransferase Ezh2. J Biol Chem. 2015;290:27604–17. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 92. Wang XW, Yang SG, Hu MW, Wang RY, Zhang C, Kosanam AR, et al. Histone methyltransferase Ezh2 coordinates mammalian axon regeneration via regulation of key regenerative pathways. J Clin Invest. 2024;134:e163145. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 93. Egervari G, Akpoyibo D, Rahman T, Fullard JF, Callens JE, Landry JA, et al. Chromatin accessibility mapping of the striatum identifies tyrosine kinase FYN as a therapeutic target for heroin use disorder. Nat Commun. 2020;11:4634. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 94. Gan L, Yang Y, Li Q, Feng Y, Liu T, Guo W. Epigenetic regulation of cancer progression by EZH2: from biological insights to therapeutic potential. Biomark Res. 2018;6:10. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 95. Sanulli S, Justin N, Teissandier A, Ancelin K, Portoso M, Caron M, et al. Jarid2 methylation via the PRC2 complex regulates H3K27me3 deposition during cell differentiation. Mol Cell. 2015;57:769–83. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 96. Loh CH, van Genesen S, Perino M, Bark MR, Veenstra GJC. Loss of PRC2 subunits primes lineage choice during exit of pluripotency. Nat Commun. 2021;12:6985. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 97. Petracovici A, Bonasio R. Distinct PRC2 subunits regulate maintenance and establishment of Polycomb repression during differentiation. Mol Cell. 2021;81:2625–39.e2625. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 98. Zhang J, Roberts JM, Chang F, Schwakopf J, Vetter ML. Jarid2 promotes temporal progression of retinal progenitors via repression of Foxp1. Cell Rep. 2023;42:112237. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 99. Sreeshma B, Devi A. JARID2 and EZH2, the eminent epigenetic drivers in human cancer. Gene. 2023;879:147584. [ DOI ] [ PubMed ] [ Google Scholar ] 100. Verberne EA, Goh S, England J, van Ginkel M, Rafael-Croes L, Maas S, et al. JARID2 haploinsufficiency is associated with a clinically distinct neurodevelopmental syndrome. Genet Med. 2021;23:374–83. [ DOI ] [ PubMed ] [ Google Scholar ] 101. Peng Q, Bizon C, Gizer IR, Wilhelmsen KC, Ehlers CL. Genetic loci for alcohol-related life events and substance-induced affective symptoms: indexing the “dark side” of addiction. Transl Psychiatry. 2019;9:71. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 102. Uffelmann E, Huang QQ, Munung NS, De Vries J, Okada Y, Martin AR, et al. Genome-wide association studies. Nat Rev Methods Primers. 2021;1:1–21. [ Google Scholar ] 103. Nyegaard M, Demontis D, Foldager L, Hedemand A, Flint TJ, Sørensen KM, et al. CACNA1C (rs1006737) is associated with schizophrenia. Mol Psychiatry. 2010;15:119–21. [ DOI ] [ PubMed ] [ Google Scholar ] 104. Sklar P, Smoller JW, Fan J, Ferreira MAR, Perlis RH, Chambert K, et al. Whole-genome association study of bipolar disorder. Mol Psychiatry. 2008;13:558–69. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 105. Wabreha AY, Adjei N, Ladenheim B, Cadet JL, Daiwile AP. Escalated oxycodone self-administration is associated with expression of voltage gated and calcium activated potassium channels in the mesocorticolimbic system in rats. Front Pharmacol. 2025;16:1653356. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 106. Xie C, Kessi M, Yin F, Peng J. Roles of KCNA2 in neurological diseases: from physiology to pathology. Mol Neurobiol. 2024;61:8491–517. [ DOI ] [ PubMed ] [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplemental File 1 (1.3MB, pdf) Supplemental File 2 (3.7MB, xlsx) Supplemental File 3 (577.5KB, xlsx) Supplemental File 4 (4.7MB, xlsx) Supplemental File 5 (888.6KB, xlsx) Supplemental File 6 (42.7KB, xlsx) Articles from Translational Psychiatry are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (2.0 MB) Cite Collections Permalink RESOURCES Similar articles Cited by other articles Links to NCBI Databases Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top

Related documents

• Experimental assessment of intravenous lipid emulsion in the management of acute fentanyl poisoning#190999
• The Donovan Memorial Lecture: 'Fentanyl-Plus': A New Era of Fentanyl Polydrug Combinations in the North American Overdose Crisis.#159622
• Sequence-dependent modulation of hepatorenal biochemical markers following Artemether-Lumefantrine and Sulfadoxine-Pyrimethamine exposure in Wistar rats#383457
• The Illicit Fentanyl Crisis in the United States: Root Causes, Federal Policy Responses, and State-Level Solutions#201198
• Effects of Different Doses of Dexmedetomidine on Fentanyl-induced Cough#350042