A computational strategy for finding novel targets and therapeutic compounds for opioid dependence.

Opioids are widely used for treating different types of pains, but overuse and abuse of prescription opioids have led to opioid epidemic in the United States. Besides analgesic effects, chronic use of opioid can also cause tolerance, dependence, and even addiction. Effective treatment of opioid addiction remains a big challenge today. Studies on addictive effects of opioids focus on striatum, a main component in the brain responsible for drug dependence and addiction. Some transcription regulators have been associated with opioid addiction, but relationship between analgesic effects of opioids and dependence behaviors mediated by them at the molecular level has not been thoroughly investigated. In this paper, we developed a new computational strategy that identifies novel targets and potential therapeutic molecular compounds for opioid dependence and addiction. We employed several statistical and machine learning techniques and identified differentially expressed genes over time which were associated with dependence-related behaviors after exposure to either morphine or heroin, as well as potential transcription regulators that regulate these genes, using time course gene expression data from mouse striatum. Moreover, our findings revealed that some of these dependence-associated genes and transcription regulators are known to play key roles in opioid-mediated analgesia and tolerance, suggesting that an intricate relationship between opioid-induce pain-related pathways and dependence may develop at an early stage during opioid exposure. Finally, we determined small compounds that can potentially target the dependence-associated genes and transcription regulators. These compounds may facilitate development of effective therapy for opioid dependence and addiction. We also built a database (http://daportals.org) for all opioid-induced dependence-associated genes and transcription regulators that we discovered, as well as the small compounds that target those genes and transcription regulators.

Introduction

Opioids such as morphine have long been used as mainstay therapy for treating different types of chronic severe pains such as cancer pain, noncancer-related pain, and neuropathic pain [1]. As a result of relaxation on restriction of prescription opioids for treating chronic noncancer pain and promotion of opioids in treatment by pharmaceutical industry, practitioners and many organizations, non-medical use and abuse of prescription opioids has been increasing rapidly in the United States [2], which in turn leads to opioid epidemic [3].

Despite that opioids have beneficial analgesic effects of alleviating acute and chronic pain, chronic opioid use can lead to adverse side effects including tolerance, hyperalgesia, withdrawal reactions, and even dependence [2, 4, 5]. In order to achieve optimal pain management, many studies have been carried out to elucidate mechanisms underlying both the beneficial as well as the adverse effects of opioids. Plenty of evidence has shown that crosstalks between neuronal signaling, immune responses and chemokines play significant roles in the pain pathways responsive to opioids such as morphine, and these signaling networks can lead to both behavioral and structural changes in the brain [6, 7].

Studies on effects of opioids to relieve pain mainly focus on brain regions such as periaqueductal gray, rostral ventromedial medulla, and dorsal root ganglia [6], while investigations of addictive effects of opioids mostly focus on striatum, since striatum is a main component of the reward system responsible for drug dependence and addiction. Striatum receives several types of neuronal inputs from prefrontal cortex (PFC), ventral tegmental area (VTA), and other areas of the brain [8], and over time, drug dependence and addiction can be reinforced [9]. At the molecular level, many transcription factors (TFs) have been associated with behaviors related to drug abuse and addiction. Some TFs known to play key roles in drug addiction include ΔFosB, cyclic AMP-responsive element binding protein (CREB), NF-κB, and MEF2 [7]. For instance, opioids can reduce Fos expression in the direct pathway striatal neurons [10]. Moreover, drugs of abuse can also alter gene transcription and induce addiction by epigenetic mechanisms [7].

Despite that different drugs of abuse often elicit similar behavioral responses in animals and humans, molecular mechanisms underlying addiction induced by different drugs can be distinctly different [10]. Opioids such as morphine and heroin increase dopamine level in nucleus accumbens, the main component of the ventral striatum, through activation of dopaminergic neurons in VTA [10]. It is popularly believed opioids inhibition of GABAergic neurons in VTA is another contributing factor of disinhibition of dopamine in VTA and increased rewarding effects in nucleus accumbens [11-16]. Opioids can also increase glutamate release in the nucleus accumbens, which results in changes in synaptic plasticity such as decreased dendritic branching and spine density [17].

Despite all these efforts, however, much remains to be known about molecular connections between the genes and pathways activated by opioids in the pain-related processes and those involved in opioid dependence and addiction. Elucidating such connections can not only shed light on the mechanisms which contribute to the opioid epidemic, but may also allow people to identify better candidate targets for therapeutic interventions to prevent opioid dependence and addiction during pain management.

In order to treat opioid addiction, opioid antagonists such as naltrexone have been used to treat opioid addiction for several decades. However, they have shown limited efficacy in relapse prevention [18, 19]. Recently, alternative approaches such as receptor-based therapeutic strategies have been proposed which aim to target receptors including G protein-coupled receptors (GPCRs) such as μ-, δ-, and κ-opioid receptors, chemokine receptors, as well as neuroimmune receptors such as Toll-like receptors 4 (TLR4) [20-22]. However, despite all the efforts and advances in understanding the mechanism of addiction in the past decades, they have not led to development of effective new anti-addiction agents [23]. Therefore, finding novel targets and strategies is needed for treating opioid dependence and addiction.

In this work, we developed a new computational strategy for identification of novel genome-wide targets and potential therapeutic treatments for opioid dependence. In particular, this strategy involves first detecting genes and pathways induced by either morphine or heroin which are associated with dependence, then identification of small compounds which can target the dependence-associated genes responsive to the opioids. We analyzed a data set generated from a previous study by Piechota et al [24], in which, gene expression responses in mouse striatum were investigated after mice were exposed to each of the six drugs of abuse (i.e., morphine, heroin, cocaine, methamphetamine, ethanol, or nicotine). Using the strategy we developed, we identified morphine and heroin-induced genes and pathways associated with dependence-related behaviors such as physical dependence and psychological dependence, as well as transcription regulators which can potentially regulate these genes. Finally, we identified small compounds which can potentially target some of the dependence-related genes and transcription factors. A database for all the dependence associated genes and transcription regulators, along with small compounds for targeting those genes and transcription regulators, is available at http://daportals.org. Our findings can facilitate identification of novel candidate gene targets as well as potential therapeutic interventions for treating opioid dependence and addiction.

Results

Identification of genes and patterns induced by either morphine or heroin

In order to identify genes induced by either morphine or heroin over time, we applied a local regression method to gene expression microarray data collected from mice treated by each opioid at different time points [24]. Using this approach, we found that 423 genes were differentially expressed (DE) after morphine administration in mice, while 608 genes were differentially induced by heroin. Using k-means clustering, we identified 6 expression patterns among genes induced by each opioid (Figs 1 and 2), with genes upregulated and downregulated at 3 phases: Immediate-Early (IE), Middle (M), and Late (L), respectively. Compared to morphine, heroin induced twice as many genes in the IE and L phases, respectively (Figs 1 and 2), indicating that heroin elicits neurobiological responses in mice not only faster but also longer lasting than morphine.

Fig 1

Genes differentially expressed after morphine exposure in mouse striatum.

(A) Patterns of differentially expressed genes induced by morphine. (B-C) These plots show six genes upregulated (B) and downregulated (C) by morphine in the IE, M, and L phase, respectively.

Fig 2

Genes differentially expressed after heroin exposure in mouse striatum.

(A) Patterns of differentially expressed genes induced by heroin. (B-C) These plots show six genes upregulated (B) and downregulated (C) by heroin in the IE, M, and L phase, respectively.

Enriched Gene Oncology (GO) and KEGG terms among differentially expressed genes (DEGs)

Our GO and KEGG analyses revealed that many DEGs induced by either morphine or heroin in the mouse striatum were involved in the immune and neuronal processes and pathways (Tables 1 and 2, S1 and S2 Tables). Many of these biological processes are previously known to play important roles in opioid-mediated pain pathways in brain regions such as dorsal root ganglia [6]. However, despite that the neuroimmune signaling processes induced by either morphine or heroin in mouse striatum share some similarities, it is also apparent that the two opioids elicit distinct neurobiological responses in the animals which we will detail below:

Table 1

Significantly enriched biological processes and pathways induced by morphine which were involved in opioid-mediated pain pathways and corresponding literature support.

In the “Phase” column, Up-IE, Up-M, and Up-L represent upregulated in the IE, M, and L phase, respectively, while that Down-IE, Down-M, and Down-L represent downregulated in the IE, M, and L phase, respectively.

	GO/KEGG term	Phase	Effect (Literature support)	DEGs Associated with Dependence and Other Harmful Effects
Immune system	Pattern recognition receptor signaling pathway	Up-M	Proinflammatory responses, tolerance [6, 22]	Irak1: phys dep; Ptafr: phys dep.
	Activation of innate immune response		Proinflammatory responses, tolerance [6, 22]	Irak1: phys dep; Ptafr: phys dep.
	Toll-like receptor signaling pathway		Allodynia and hyperalgesia [25, 26];NF-κB activation [27–29]	Irak1: phys dep.
	Response to xenobiotic stimulus	Up-L	Tolerance [30, 31]
	Negative regulation of NF-κB TF activity		AnalgesiaNF-κB activation [6]
	MAPK signaling pathway	Down-M	Anti-inflammatory responses Analgesia [32]
	Humoral immune response
	Induction of positive chemotaxis	Down-L	Nociceptive pathways [22]
Neuronal signaling	Cell projection morphogenesis	Up-L	Structural plasticity [33]	Numb: acute, dep, HCC, phys dep, phys harm
	Calcium signaling pathway	Down-IE	Analgesia [22]
	Synaptic transmission, glutamatergic	Down-M	Analgesia [34–36]	Grin1: phys dep, pleasure
	Ensheathment of neurons		Nociceptive pathways [22]	Cldn5 (HM*): phys dep.
	Synaptic transmission, GABAergic		Proinflammatory; tolerance [6, 37]
	Glutamate receptor signaling pathway		Analgesia [6]	Cacng7: dep, phys dep, pleasure
	Sensory organ development		Analgesia [38, 39]	Anp32b: phys dep.
	Negative regulation of neuron apoptotic process		Neuronal apoptosis [6]	Grin1: phys dep, pleasure
Other key processes	Protein dephosphorylation	Up-IE	Alleviating inflammatory hyperalgesia [6, 40]	Dusp12 (HM): phys dep.
	positive regulation of autophagy		Production of proinflammatory cytokines, tolerance [6]	Plekhf1: psycho dep.
	Regulation of programmed cell death		Production of proinflammatory cytokines, tolerance [6]	Dapk1 (HM): psycho dep; Plekhf1 (HM): psycho dep; Pim3 (HM): phys dep.

Note: the full names of the behaviors can be found in S3 Table. HM* indicates that the DEG is induced by both heroin and morphine.

Table 2

Significantly enriched biological processes and pathways induced by heroin which were involved in opioid-mediated pain pathways and corresponding literature support.

All of the abbreviations used in this table can be found in the legend of Table 1.

	GO/KEGG term	Phase	Effect(Literature support)	DEGs Associated with Dependence and Other Harmful Effects
Immune system	Response to biotic stimulus	Up-IE	Nociceptive pathways [6]	Ace: dep, phys dep; Baiap2: pleasure, psycho dep; Stab1 (HM): phys dep
	MyD88-dependent toll-like receptor signaling pathway	Down-L	Anti-inflammatory effect; analgesia [30, 41]
	T-helper 1 type immune response (Il4, Tlr6, Il27)
	Positive regulation of I-κB kinase/NF-κB signaling
	Inflammatory response
	Regulation of MAP kinase activity
	Microglial cell activation		Anti-inflammatory responses; analgesia [6]
	Regulation of granulocyte chemotaxis		Nociceptive pathways [22]
Neuronal signaling	Cell morphogenesis involved in neuron differentiation	Up-IE	Nociceptive pathways [22]	Baiap2: pleasure, psycho dep
	Regulation of nervous system development		Nociceptive pathways [22]	Ace: dep, phys dep; Ncs1: dep; Baiap2: pleasure, psycho dep
	Regulation of excitatory postsynaptic membrane potential		Tolerance and hyperalgesia [22]	Sez6: dep
	Potassium ion transport	Up-M	Proanalgesic effect [22]
	Cation transmembrane transport		Tolerance and hyperalgesia [32]	Slc25a42: dep
	Sodium ion transmembrane transport	Up-L	Allodynia and hyperalgesia [22]
	Regulation of synapse organization		Analgesia [22]
	Regulation of neuron death	Down-IE	Linked to anti-inflammatory response; analgesia [42]	Bag1: dep, pleasure; Tfap2d: dep, pleasure, psycho dep
Other key processes	Protein autophosphorylation	Up-IE	Hyperalgesia [22]	Dapk1 (HM): psycho dep; Pim3 (HM): phys dep
	Negative regulation of receptor recycling	Up-L	Tolerance [22]	Pcsk9: pleasure, psycho dep
	Response to cAMP	Down-M	Analgesia [6]

Immune and neuronal responses to morphine in mouse striatum

Immune responses

Most genes involved in the immune system were induced by morphine in the M phase. As shown in Table 1, some genes participating in anti-inflammatory processes were responsive to morphine, e.g., genes involved in negative regulation of NF-κB transcription factor activity were upregulated, and those involved in MAPK pathway, humoral immune response, and induction of positive chemotaxis were downregulated. Since anti-inflammatory pathways are known to play key roles in opioid-induced pain relief [6], these results are consistent with the fact that morphine induces analgesia effect in animals.

On the other hand, some genes participating in proinflammatory responses were also upregulated by morphine, including, e.g., genes involved in Toll-like receptor signaling pathway, and activation of innate immune response. Notably, previous evidence showed that proinflammatory responses play central roles which contribute to tolerance during chronic opioid exposure [7]. Genes involved in response to xenobiotic stimulus were also upregulated in the L phase, in line with the fact that these genes are essential in sensing that the cells are under the ‘insult’ of the drug.

Together, these results suggest that crosstalks between genes involved in proinflammatory and anti-inflammatory pathways have already initiated in mouse striatum during short-term morphine administration.

Neuronal responses

Our results also showed that genes participating in neuronal responses were induced by morphine (Table 1), which is not surprising, given known neurological effects of morphine in the brain. For example, positive regulation of autophagy and programmed cell death were upregulated in the IE phase, and also genes involved in glutamatergic synaptic transmission, glutamate receptor signaling pathway, and sensory organ development were all downregulated in the M phase. Notably, all these neuronal signaling events have been implicated in the analgesia induced by morphine (Table 1), again supporting the notion that morphine can induce analgesia in the treated mice. Also, genes involved in cell projection morphogenesis were upregulated in the L phase, in line with the evidence that morphine can induce structural changes in mice during chronic exposure [33].

Together, our results suggest that even during short-term morphine exposure, complex crosstalks between genes involved in proinflammatory and anti-inflammatory pathways, neuronal signaling, and the chemokine system have already initiated in mouse striatum, which can contribute to analgesia and/or tolerance effects if exposure of the drug lasts longer.

Immune and neuronal responses to heroin in mouse striatum

Our results showed that distinct from morphine, many immune genes induced by heroin were downregulated in the L phase (Table 2), which included MyD88-dependent Toll-like receptor signaling pathway, microglial cell activation, T-helper 1 type immune response (Il4, Tlr6, Il27), positive regulation of I-κB kinase/NF-κB signaling, regulation of granulocyte chemotaxis, inflammatory response, and regulation of MAP kinase activity. Notably, all these pathways are known to be active in proinflammatory responses during chronic opioid exposure [6]. Therefore, downregulation of these pathways indicates that heroin induces strong anti-inflammatory response and thus elicits strong analgesic effects in the mice.

Genes involved in neuronal activities, such as regulation of nervous system development and excitatory postsynaptic membrane potential (IE phase) and cation transmembrane transport (M phase) were upregulated by heroin, whereas genes involved in regulation of neuron death were downregulated (IE phase) (Table 2). These results also agree with previous findings that neuronal responses play active roles in opioid-related pain process [6].

Also, we noticed that genes involved in cell morphogenesis involved in neuron differentiation, positive regulation of axonogenesis, and regulation of synapse organization were upregulated among IE and L phases after heroin exposure. These results are supported by the previous evidence that drugs of abuse can induce changes in structural plasticity in animals during chronic exposure [43].

Other key biological responses

Our results also showed that heroin induced genes participating in other key biological processes involved in pain-related pathways, e.g., genes involved in protein autophosphorylation are upregulated in the IE phase. Since protein phosphorylation is known to play key roles in desensitization and implicated in opioid-induced hyperalgesia [22], our results suggest that these genes contribute to dependence induced by chronic use of heroin; this speculation is confirmed by our association analysis described below.

Association of morphine- and heroin-induced DEGs with harmful effects of drugs of abuse

In order to find out whether DEGs induced by either morphine or heroin are associated with various harmful effects linked to drugs of abuse, we conducted the association analysis between expression levels of morphine- or heroin-induced DEGs and twelve DA-related harmful effects including dependence, physical dependence, psychological dependence, pleasure, physical harm, social harm, health care cost, and conditioned place preference (See S3 Table and Methods for details).

Using this approach, we detected 44 morphine-induced DEGs and 61 heroin-induced DEGs significantly associated with dependence-related behaviors at the nominal level of significance (p < 0.05) (S4 and S5 Tables). Among these dependence-associated DEGs, 9 were induced by both morphine and heroin, of which 6 were induced in the IE phase. These results can be searched in our database at http://daportals.org.

Next, we investigated whether the dependence-associated DEGs induced by either morphine or heroin were involved in the pain-related neuroimmune pathways mediated by opioids. As shown in Tables 1 and 2, we found that a significant number of the dependence-associated DEGs induced by the opioids were involved in the pain-related neuroimmune pathways. Moreover, some of these genes could be induced by both morphine and heroin, e.g., Dapk1, Plekhf1, Pim3, and Dusp12 were upregulated by both morphine and heroin in the IE phase, and were associated with psychological dependence and physical dependence, respectively.

Detection of potential transcription regulators that regulate the opioid-induced dependence-associated DEGs

In order to find out whether the dependence-associated DEGs induced by each opioid were co-regulated by any TFs or epigenetic factors, we performed the TF and epigenetic factor binding site enrichment test using the ENCODE ChIP-Seq significance tool [44]. Our analysis showed that 17 transcription regulators potentially modulated morphine-responsive dependence-associated DEGs (S6 Table), while that 12 transcription regulators modulated heroin-responsive dependence-associated DEGs (S7 Table). More details about our results are described below.

Transcription regulators detected after morphine exposure

Known TFs associated with dependence and addiction

More than half of the TFs we detected regulating dependence-associated DEGs after morphine exposure are previously known to play important roles in drug dependence and addiction. For example, we found that MEF2A upregulated four DEGs which were associated with physical dependence in the M phase, while that MEF2C upregulated one DEG associated with dependence in the IE phase after morphine exposure, agreeing with the evidence that MEF2 is crucial in inducing behavioral changes after exposure to drugs of abuse [7].

Novel TFs induced by morphine

Our results also showed a few novel TFs activated by morphine. In order to quantify the magnitude of the effects the detected transcription regulators on the dependence-related behaviors induced by the opioids, we developed a scoring metric, called dependence score, based on the total fold changes of the dependence-associated DEGs co-regulated by each regulator (see details in Methods). Using this approach, we found that E2f6, which potentially regulated 13 DEGs associated with physical dependence after morphine exposure, had the highest dependence score of 17.93. Also, we detected ZBTB33 and ZKSCAN1 associated with physical dependence having high dependence scores (11 and 6.8, respectively) after morphine exposure. In particular, ZBTB33 encodes a transcriptional regulator Kaiso which can promote histone deacetylation and decrease expression levels of its target genes [45], consistent with our results showing that ZBTB33 downregulates 8 genes in the M phase. Also, Zkscan1 encodes a member of the Kruppel C2H2-type zinc-finger family of proteins, which has been implicated in regulating the expression of GABA type-A receptors in the brain [46].

Epigenetic factors

Our results showed that about half (8 out of 17) of the detected transcription regulators after morphine exposure were epigenetic factors (S6 Table). Literature search (S8 Table) suggests that these factors including HDAC8 (encoding histone deacetylase 8) and HDAC6 (encoding histone deacetylase 6) play important roles in histone acetylation, histone methylation, DNA methylation, and chromatin remodeling [45, 47–52], in line with the previous evidence that all these epigenetic events have been implicated in the neurobiological responses to drugs of abuse in the brain [7]. Notably, among all the epigenetic factors, SAP30 which encodes a component of the histone deacetylase complex had the highest dependence score of 15.49 and can potentially co-regulate 11 morphine-responsive DEGs associated with physical dependence.

Transcription regulators detected after heroin exposure

DEGs upregulated by EGR1 and CREB1 were associated with psychological dependence in the IE phase after administration of heroin, consistent with the fact that both EGR1 and CREB are key TFs which regulate genes involved in dependence-related behavioral responses during exposure of opioids as well as other drugs of abuse [5, 7]. Moreover, our results showed that Polr2a (encoding the largest subunit of RNA polymerase II), EGR1, and CREB1 were associated with social and physical harm with the highest scores (> 20), suggesting that these TFs play key roles in the biological mechanism that underlies the higher social and physical harm caused by heroin compared to other drugs of abuse [4].

Novel TFs induced by heroin

We found that E2F6 (which showed the highest dependence score in morphine) also had the highest dependence score of 6.52, potentially regulating four DEGs associated with psychological dependence after exposure to heroin.

Epigenetic regulators

Similar to morphine, more than 40% of the detected transcription regulators were epigenetic factors after heroin exposure (S7 Table). All these factors have been known to play major roles in histone and DNA methylation, and chromatin remodeling methylation [50-53]. Notably, 3 epigenetic factors CTCF, EZH2, and SUZ12 were identified during both morphine and heroin exposure.

Taken together, our results suggest that distinct transcriptional regulatory mechanisms are responsive to exposure of morphine and heroin in mouse striatum, and that epigenetic regulation plays major roles after exposure to both opioids. Further investigation is needed to elucidate the roles of the novel TFs (e.g., E2F6) with high dependence scores after exposure of either morphine or heroin.

Finding small compounds which can target the DEGs and the dependence-associated transcription regulators induced by the opioids

To facilitate development of therapeutic interventions for treating morphine or heroin dependence, we developed a strategy which allowed us to identify small compounds that can target the opioid-induced dependence-associated DEGs and their potential transcription regulators (see details in Methods). Our results are shown in S9 and S10 Tables for morphine and heroin, respectively. Among all the small compounds we identified, we found that Calmidazolium could target 8 morphine-induced dependence-associated DEGs (S9A Table), and that Securinine could target E2F6 which had the highest dependence score after morphine exposure (Table 3). Also, we identified that Phenacetin and Buspirone could target 11 and 4 heroin-induced dependence-associated DEGs, respectively (Table 4, S1 Fig), and that Meclofenoxate could target E2F6 which had the highest dependence score also after heroin exposure (S10B Table).

Table 3

Small compounds negatively correlated with the potential transcription regulators after morphine administration.

All of the abbreviations used in this table can be found in the legend of Table 1.

Compound	Transcription Regulators	Fold Change	Phase	Associated Harmful Effects
Benserazide	TAF1	2.9	Up-IE	Dep, pleasure
	MEF2A	5.34	Up-M	Phys dep
	ZKSCAN1	-6.81	Down-M	Phys dep
Piperacetazine	SUZ12	1.2	Up-IE	Dep
	ZKSCAN1	-6.81	Down-M	Phys dep
Securinine	MEF2C	1.7	Up-IE	Dep
	BRF2	1.39	Up-M	Phys dep
	E2F6	-17.93	Down-M	Phys dep
Isocorydine	MEF2A	5.34	Up-M	Phys dep
Gabapentin	CTCF	2.97	Up-IE	Chronic, dep, phys harm
	MEF2A	5.34	Up-M	Phys dep
	HDAC6	-4.84	Down-M	Acute, phys dep

Table 4

Small compounds negatively correlated with dependence-associated DEGs after heroin administration.

All of the abbreviations used in this table can be found in the legend of Table 1.

Compound	DEGs	Fold Change	Phase	Harmful Effects Associated with DEGs
Prilocaine	DAPK1	1.63	Up-IE	Acute, intox, psycho dep
	ACE	1.34	Up-IE	Acute, dep, HCC, phys dep, phys harm, soc
	CEP350	1.29	Up-IE	Dep, HCC, soc harm
	FAAH	-1.42	Down-IE	Acute, dep, phys harm, pleasure, soc harm, soc
	MPDU1	-1.39	Down-L	Dep, HCC, phys dep, phys harm, psycho dep
Phenacetin	PLEKHF1	1.64	Up-IE	Psycho dep
	RPL7L1	1.56	Up-IE	Dep, HCC, phys dep, phys harm, soc harm, soc
	EXTL1	1.24	Up-IE	Dep, soc harm, soc
	SPON2	1.23	Up-M	Acute, dep, HCC, phys dep
	TMTC4	1.19	Up-M	Dep, HCC, phys dep, phys harm
	PCSK9	1.52	Up-L	Chronic, CPP, pleasure, psycho dep
	FBRSL1	-1.17	Down-IE	Psycho dep
	ACSL4	-1.34	Down-M	CPP, dep, HCC, phys dep
	FEM1C	-1.32	Down-M	Hcc, phys dep, soc harm
	TOR3A	-1.21	Down-M	Acute, dep, HCC, phys dep, phys harm, soc harm, soc
Procyclidine	KLF2	2.31	Up-IE	CPP, dep, HCC, phys dep, phys harm
	FEM1C	-1.32	Down-M	Hcc, phys dep, soc harm
Spiradoline	PLEKHF1	1.64	Up-IE	Psycho dep
	KIF23	1.36	Up-IE	Dep, HCC, phys harm, soc harm, soc
	CEP350	1.29	Up-IE	Dep, HCC, soc harm
	PTBP3	-1.19	Down-L	Chronic, CPP, intox, phys dep, phys harm, psycho dep
Buspirone	NCS1	1.34	Up-IE	Dep, HCC, soc harm
	TBC1D2B	1.31	Up-IE	Dep, HCC, phys dep, phys harm, soc harm, soc
	SLC35D1	-1.29	Down-M	Acute, phys dep
	TOR3A	-1.21	Down-M	Acute, dep, HCC, phys dep, phys harm, soc harm, soc

Comparison of our approach and findings with those in a previous work by Piechota et al. [24]

Since the data set we analyzed in this work was generated from a previous study in [24], it is worthwhile to compare our approach and results with those reported in [24]. There are two main differences between the two studies. First, Piechota et al. employed two-way ANOVA to identify genes differentially expressed in mouse striatum responsive to some of the six drugs of abuse. The limitation of using two-way ANOVA to analyze time-course data is that time is considered as a factor in the ANOVA analysis and thus the trends of expression levels of genes (over time) are ignored. We, instead, employed a local regression smoothing technique to identify differentially expressed genes responsive to the drug exposure, which was able to capture the trends of the expression levels of the genes over time. Second, Piechota et al. focused on identification and characterization of gene expression patterns induced by multiple drugs of abuse, whereas we were interested in the genes and patterns induced by either morphine or heroin, and particularly when the genes were involved in pain pathways and associated with dependence.

Using our approach, we identified 423 and 608 DE genes induced by morphine and heroin, respectively, in mouse striatum (Figs 1 and 2). Since only the two-way ANOVA was used to identify the 42 DE genes in [24], it is unknown which of the six drug(s) induced the genes. Regardless, we compared the DE genes in the two studies, and found that 11 morphine-induced genes and 17 heroin-induced genes in our work overlapped with the previous study. We also compared the enriched functional groups of genes identified by the two studies, and found that the majority of enriched GO groups identified in the previous study were also discovered in our work (see S11 Table), e.g., protein phosphatase activity, apoptosis, anatomical structure morphogenesis, calcium ion binding, GTPase mediated signal transduction, and transmembrane transporter activity. However, as shown in S1 and S2 Tables, we also identified many more significantly enriched GO and KEGG groups of genes than those in [24], such as those involved in the neuroimmune signaling processes (Tables 1 and 2).

Discussion

Many studies have been conducted which intend to delineate pain-related pathways induced by opioids. However, much remains to be known about the molecular connection between these opioid-mediated pain pathways and those playing key roles in drug dependence and addiction. Dissecting these pathways can facilitate identification of candidate targets for developing effective therapeutic interventions which ideally can target opioid tolerance and dependence while preserving opioid analgesic effect.

In this study, we developed a computational strategy to identify candidate dependence-associated DEGs induced by either morphine or heroin, as well as to find small compounds which could target these genes for treating dependence of the opioids. Using this strategy, we analyzed a time-course gene expression microarray data set generated previously to investigate gene expression patterns responsive to various drugs of abuse in mouse striatum [24]. In particular, we first employed a local regression technique to detect genes differentially expressed over 8 hours of time in mouse striatum after either morphine or heroin exposure. Then, we performed correlation analysis to identify morphine or heroin-induced DEGs which were associated with twelve harmful effects including dependence commonly linked to drugs of abuse. Furthermore, we detected potential transcription regulators including TFs and epigenetic factors that regulated the dependence-associated DEGs using an ENCODE enrichment tool. Finally, to facilitate the identification of candidate targets and development of effective therapy for morphine and heroin-induced dependence, we identified small compounds which could potentially target against some of the detected dependence-associated DEGs and transcription regulators.

Using the approach described above, we found that a significant number of the DEGs responsive to either morphine or heroin in mouse striatum were involved in the neuroimmune signaling pathways, which are typically activated in the pain-related pathways during chronic opioid use previously identified in other brain areas including periaqueductal gray, rostral ventromedial medulla, and dorsal root ganglia [6]. Using correlation analysis, we found that a considerable portion of the pain pathway-related DEGs, previously known to play active roles in opioid analgesia, tolerance, hyperalgesia, and allodynia, were associated with the harmful effects (such as dependence) linked to morphine and heroin as well as many other drugs of abuse, e.g., Irak1 (encoding interleukin-1 receptor-associated kinase 1) in the enriched Toll-like receptor signaling pathway (Table 1) induced by morphine was correlated with physical dependence at a nominal level of significance (p < 0.05). Toll-like receptor signaling pathway has been known to play crucial roles in proinflammatory signaling and tolerance to opioid analgesia. We also noticed that some dependence-associated DEGs could be induced by both morphine and heroin, e.g., among DEGs upregulated in the IE phase, Dapk1 and Plekhf1 were correlated with psychological dependence, while Dusp12 and Pim3 were associated with physical dependence after exposure to both morphine and heroin. It is unclear what roles these genes (induced by both opioids) play when mice were first exposed to morphine, then switched to heroin later on, a scenario commonly seen among human drug abusers.

Despite the similarities in gene expression responses induced by both morphine and heroin in mouse striatum, differences between the two opioids are also obvious. For example, a large number of the DEGs involved in immune signaling were downregulated in the L phase after heroin exposure, as opposed to morphine exposure, suggesting that heroin elicited strong anti-inflammatory responses in the L phase and thus induced acute analgesic effects in the mice. Considering that heroin is diamorphine which is rapidly converted into several psychoactive metabolites including 6-mono-acetylmorphine (6MAM) and morphine in humans and mice [54, 55], it is not surprising that morphine and heroin elicit similar responses in vivo. However, several factors may account for the differences between the two opioids. First, evidence suggests that heroin and its metabolite 6MAM can elicit different responses in mice than morphine, e.g., the former, but not the latter, can elicit acute toxic effect on locomotor activity, particularly at high doses [55], suggesting that 6MAM may contribute to the differential responses in mice induced by heroin and morphine. Second, since both heroin and 6MAM have much higher degrees of lipophilicity than morphine, after intraperitoneal (i.p.) injection as done in [24], heroin and 6MAM can cross the blood-brain barrier much faster than morphine [56] and thus some of the late gene expression response to heroin (but not morphine) exposure could be captured by the gene expression profiling which measured gene expression levels at no later than 8 hours after the opioid injection. Finally, since different dosages of morphine (20 mg/kg) and heroin (10 mg/kg) were used for i.p. injection of the mice, it is difficult to assess to what extend this dosage effect affects the gene expression responses in mouse striatum induced by morphine and heroin.

Among detected transcription regulators that potentially regulate the dependence-associated DEGs, MEF2A induced by morphine as well as EGR1 and CREB1 by heroin are known to play crucial roles in drug addiction. We also found that more than 40% of the detected transcription regulators are epigenetic factors after both morphine and heroin exposure, including HDAC8 and HDAC6 activated by morphine, which supports the previous notion that epigenetic factors are important for addiction. Furthermore, using the dependence score, a metric we developed for measuring the extent the detected transcription regulators affect the dependence-associated DEGs, we found that E2F6 has the highest dependence scores after exposure of both morphine and heroin.

In summary, our work here intent to elucidate molecular connections between the analgesic and tolerance-related pain pathways and harmful side effects of opioid use during pain treatment. Despite the general belief that morphine is safe for managing patients with pain, our results suggest that morphine may induce tolerance to analgesia and dependence on the drug in the patients in the very early stage, which may increase the possibility of the same patients to abuse heroin thereafter, since heroin may further induce acute analgesic effects as suggested by our results. Moreover, because heroin can cause both structural and behavioral changes among patients, abusing heroin after morphine may lead to more potent dependence on the drugs among the patients.

Furthermore, we found several small compounds which could potentially target some of the dependence-associated DEGs and the detected transcription regulators induced by the opioids. In particular, we identified Securinine and Meclofenoxate which could target E2F6 in humans after exposure to morphine and heroin, respectively. These compounds can facilitate future development of effective therapeutic interventions which can target the adverse side effects of morphine and heroin, while preserving their analgesic effects.

We also compared the approach and findings from this work with those from the previous study by Piechota et al. [24] from which the data set was generated. The approaches employed by the two studies were different but complement with each other. Piechota et al. aimed to identify and characterize genes differentially induced by any of the six drugs using the two-way ANOVA analysis, while we focused on finding and investigating genes induced by either morphine or heroin (using a local regression technique) which were also associated with dependence and previously known to affect pain pathways. Comparison of the enriched functional (GO and KEGG) groups of genes identified in the two studies suggest that the majority of the GO groups discovered in the previous study were also identified in our work and that our results revealed many more biologically meaningful functional groups of genes involved in neuroimmune signaling pathways.

The limitations of our work include the following. The gene expression microarray data we analyzed spanned only 8 hours after administration of morphine and heroin, which limited our ability to discover chronic effects of the drugs. In the future study, we intend to employ the same strategy to investigate long-term effects of morphine and heroin, and to compare them with the acute effects we discovered in this work. Also, our results in this work were generated based on analyzing a single gene expression profiling data set, further biological validation of some of the genes differentially induced over time in mouse striatum by morphine and heroin should be conducted to verify our findings here. Despite these limitations, we found that our results agree well with the previously known evidence about drug abuse and addiction, suggesting our findings are valid and worth further in-depth investigation. Moreover, our work provides insight into the molecular connections between the opioid-induced pain-related pathways and the adverse harmful effects associated with morphine and heroin. Understanding such connections may facilitate development of effective therapies which allow people to target dependence-associated genes and transcription regulators at an early stage of opioid use while preserving analgesic effects of opioids.

Methods

Dataset

The gene expression microarray data set we analyzed in this work was obtained from the NCBI Gene Expression Omnibus (GEO) database under the accession number [GEO:GSE15774]. This data set was generated from a previous work described in [24], in which, gene expression alterations in mouse striatum were investigated after the mice were treated by various drugs of abuse, including morphine, heroin, methamphetamine, cocaine, nicotine and alcohol. Detailed description of the data set can be found in [24]. Briefly, after a single dose of drug administration, gene expression was obtained from the mouse striatum at 1, 2, 4, 8 hours afterwards. Meanwhile, samples from saline- and naïve-treated control group were collected at 0, 1, 2, 4, 8 hours as controls. There were three biological replicates for each drug group and each time point.

Identification of genes differentially expressed over time in mouse striatum after exposure of either morphine or heroin

In order to identify genes differentially expressed over time in mouse striatum after administration of either morphine or heroin, we employed a local regression smoothing technique [57] to estimate the smoothed time course gene expression data for each opioid. The detailed description of the strategy can be found in [58]. For each opioid, expression values of each gene (i.e., transcript) were available for 1, 2, 4, and 8 hours, and expression values for time point 0 for the corresponding genes from the naïve group were used to represent the control time point (i.e., 0 hour) for each opioid. In particular, expression values for each gene over different time points were first fitted using a local polynomial quadratic (degree = 2) model with the bandwidth optimally estimated using a leave-one-out cross validation procedure [57]. To determine whether a gene is differentially expressed over time with respect to the control time point, we calculated the simultaneous 95% confidence intervals for the fitted (or expected) intensity values using a method due to Sun and Loader [59]. The p-values were adjusted using the Bonferroni correction to account for multiple hypothesis testing. We determined a gene as differentially expressed if its expression value at any time point T relative to the control time point satisfied: 1) adjusted p-value < 0.05, and 2) fold change ≥ 1.2.

Identification of temporal patterns for DEGs induced by either morphine or heroin using cluster analysis

In order to identify temporal patterns for DEGs responsive to either morphine or heroin exposure, we applied a k-means clustering algorithm proposed by Hartigan and Wong [60] to the temporal expression values of the DEGs. The Euclidean distance was used to measure dissimilarities between different genes. A thousand iterations were performed to find an optimal partition of K clusters where K is pre-assigned. To determine an optimal number of the clusters for the DEGs, we employed the average silhouette width (ASW) as described in [61], and when K = 6, ASW is the largest for the DEGs induced by both morphine and heroin.

GO and KEGG pathway enrichment analysis

We performed GO and KEGG pathway enrichment analysis to identify biological processes and pathways that were overrepresented among DEGs in each cluster after exposure of either morphine or heroin. We performed the enrichment analysis with the GOstats R software package [62], which finds enriched functional groups using the hypergeometric test with the aid of the functional terms in the GO and KEGG databases. GO terms and KEGG pathways were considered as significantly enriched if their p-values < 0.05.

Our association analysis aimed to determine whether morphine- and heroin-induced DEGs were associated with any harmful effects of drugs of abuse. The scores which assessed the magnitudes of the twelve harmful effects of various drugs of abuse were taken from Nutt et. al. [4] and can be found in S3 Table. In particular, the scores for the harmful effects encompassing three categories, including physical harm (overall, acute, chronic), dependence (pleasure, psychological, physical), social harm (overall, health-care costs), and conditioned place preference for the drugs including morphine, heroin, cocaine, methamphetamine, ethanol, and nicotine were used to calculate the association of each harmful effect and the DEGs induced by either morphine or heroin. Specifically, let S denote a vector of the scores for harmful effect i corresponding to drugs D, where D = [morphine, heroin, cocaine, methamphetamine, ethanol, and nicotine]. Let G denote a vector of expression values of gene G corresponding to drugs D at time point j; G is a DEG induced by either morphine or heroin at time point j, but the gene is not required to be differentially induced by the other drugs at the same time point. The expression values of gene G for the drugs including cocaine, methamphetamine, ethanol, and nicotine were estimated for each drug by using the same local regression smoothing techniques as described above for morphine and heroin. Finally, we calculated the correlation between G and S using both the Pearson correlation and a quadratic polynomial regression; if the resulting p-value from any of the methods was less than 0.05, G was considered as significantly associated with S.

Identification of human transcription and epigenetic factors that potentially regulate the dependence-associated DEGs induced by each opioid

To identify potential transcription and epigenetic factors that regulate dependence-associated DEGs responsive to each opioid in mouse striatum, we employed the ENCODE ChIP-Seq significance tool [44] to identify human TFs and epigenetic factors whose binding sites were significantly enriched among the DEGs associated with dependence (i.e., dependence, psychological, and/or physical dependence) in each of the six identified clusters after either morphine or heroin exposure. The ENCODE ChIP-Seq significance tool calculates enrichment scores of the transcription regulators using the hypergeometric test, and the resulting p-values were corrected by an FDR procedure to account for multiple hypothesis testing. A 1000-base pair (bp) window upstream of the transcription start site (TSS) and downstream of the transcription termination site (TTS) were considered for each DEG. A TF or an epigenetic factor was considered as significantly enriched if its FDR p-value < 0.05.

Furthermore, to facilitate ranking the significantly enriched TFs and epigenetic factors in terms of their impact on dependence, we developed a scoring metric called the ‘dependence score’ as follows. For each transcription regulator R activated in a certain phase P, we assume that R regulates a number N of morphine or heroin-induced DEGs associated with a dependence-related harmful effect (such as dependence, psychological, or physical dependence) within phase P. As shown in Figs 1 and 2, the IE, M, and L phases correspond to 0–2 hours, 2–4 hours, and 4–8 hours after exposure to either morphine or heroin. The ‘dependence score’ for R in phase P was then defined as , where FC represents the maximum absolute fold change of the expression values of a DEG G within phase P, relative to that of the control time point, and G is regulated by R. The higher the dependence score, the more impact the transcription regulator R can have on dependence.

Using a similar concept as the dependence score, we also assigned association scores to transcription regulator R if the DEGs it regulates were also associated with other harmful effects of drugs (as shown in S3 Table). Specifically, an association score between a transcription regulator R and a harmful effect H was defined as the sum of the maximum absolute fold change of the DEGs regulated by R in phase P that were associated with H.

Finding small molecular compounds to target the opioid-induced dependence-associated DEGs and the transcription regulators

Gene expression patterns in cells can change during treatment by small-molecule drugs or compounds [63]. If a small compound has an opposite effect on transcription than opioids, the small compound has potential to reverse the gene signature induced by opioids and hence the subsequent harmful effects caused by opioids. With the availability of gene expression profiles of small molecular compounds, we were able to compare them with those of morphine and heroin, and identify small compounds with the potential for treating dependence and addiction induced by each opioid.

Specifically, we employed the following two-step strategy to find the small compounds:

Finding DEGs induced by small compounds

First, a commercial Illumina BaseSpace (former Nextbio™) software (Santa Clara, CA, USA, http://www.nextbio.com) were used to obtain the DEGs induced by small compounds in cells. In BaseSpace, most of the raw gene expression datasets involving perturbations by small compounds were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Only genes with the p-values < 0.05 and absolute fold changes >1.2 were considered as DEGs induced by small compounds. To identify top compounds that have gene expression profiles most correlated with the dependence-associated DEG or the transcription regulators induced by morphine or heroin, we searched the gene expression profiles of small compounds stored in BaseSpace through BaseSpace integrated Pharmaco Atlas search. Then, the correlation between the DEGs induced by small compounds and the opioid-induced DEGs or transcription regulators was calculated as described below.

Calculation of the correlation between the DEGs induced by small compounds and the opioid-induced DEGs or transcription regulators

For each opioid-induced dependence-associated DEG, we used its maximum absolute fold change over the measured time points (i.e., 1, 2, 4, and 8 hours) (which was defined as the ratio of the highest absolute expression value of the gene relative to that at the control time point) to represent its fold change. For each transcription regulator, we used its dependence score to represent its fold change value.

The correlation between the DEGs induced by each small compound and the opioid-induced dependence-associated DEGs or transcription regulators was calculated using the BaseSpace software. This software provided a modified form of the rank-based enrichment statistics to compare the two sets of the DEGs [64, 65]. BaseSpace pre-processed gene expression data with biomedical ontologies to enable comparison among heterogeneous datasets from different species. It also used meta-analyses to provide consistent predictions from multiple instances of similar perturbations, e.g., genes expression profiles from different cell lines induced by the same compounds [66]. All analyses using the BaseSpace software were performed with the default parameters.

S1 Fig shows an example of the significant negative correlation between the (61) heroin-induced dependence-associated DEGs and the buspirone-induced DEGs (p-value = 0.0277). Four genes were regulated by both heroin and buspirone, but in opposite directions.

Significant negative correlation between the 61 heroin-induced dependence-associated DEGs and the buspirone-induced DEGs (p-value = 0.0277).

Four genes were regulated by both heroin and buspirone, but in opposite directions.

(TIF)

Click here for additional data file.

Significantly enriched GO terms and KEGG pathways induced by morphine in different phases and their association with harmful effects of drugs.

(A) Significantly enriched GO terms induced by morphine. (B) Significantly enriched KEGG pathways induced by morphine.

(XLSX)

Click here for additional data file.

Significantly enriched GO terms and KEGG pathways induced by heroin in different phases and their association with harmful effects of drugs.

(A) Significantly enriched GO terms induced by heroin. (B) Significantly enriched KEGG pathways induced by heroin.

(XLSX)

Click here for additional data file.

The scores of the harmful effects associated with drugs of abuse (adapted from Nutt, et. al., Lancet 2007; 369: 1047–53).

(XLSX)

Click here for additional data file.

44 dependence-associated DEGs induced after morphine exposure.

(XLSX)

Click here for additional data file.

61 dependence-associated DEGs induced after heroin exposure.

(XLSX)

Click here for additional data file.

Significantly enriched transcription regulators and associated harmful effects after morphine exposure.

(DOCX)

Click here for additional data file.

Significantly enriched transcription regulators and associated harmful effects after heroin exposure.

(DOCX)

Click here for additional data file.

Supporting evidence of involvement of transcription regulators in drug dependence and addiction from literature.

(DOCX)

Click here for additional data file.

Small compounds negatively correlated with dependence-associated DEGs and the potential transcription regulators after morphine administration.

(A) Small compounds negatively correlated with dependence-associated DEGs induced by morphine. (B) Small compounds negatively correlated with the transcription regulators induced by morphine.

(XLSX)

Click here for additional data file.

Small compounds negatively correlated with dependence-associated DEGs and the potential transcription regulators after heroin administration.

(A) Small compounds negatively correlated with dependence-associated DEGs induced by heroin. (B) Small compounds negatively correlated with the transcription regulators induced by heroin.

(XLSX)

Click here for additional data file.

Comparison of the enriched functional terms in this work with those in Piechota, et al. (2011).

(XLSX)

Click here for additional data file.