|Home | About | Journals | Submit | Contact Us | Français|
Chronic morphine administration may alter the expression of hundreds to thousands of genes. However, only a subset of these genes is likely involved in analgesic tolerance. In this report we utilized a behavior-genetics strategy to identify candidate genes specifically linked to the development of morphine tolerance. Two inbred genotypes (C57BL/6J (B6), DBA2/J (D2)) and two reciprocal congenic genotypes (B6D2, D2B6) with the proximal region of chromosome 10 (Chr10) introgressed into opposing backgrounds served as the behavior genetic filter. Tolerance following therapeutically-relevant doses of morphine developed most rapidly in the B6 followed by the B6D2 genotype and did not develop in the D2 mice and only slightly in the D2B6 animals indicating a strong influence of the proximal region of Chr10 in the development of tolerance. Gene expression profiling and pattern matching identified 64, 53, 86 and 123 predisposition genes and 81, 96, 106, and 82 tolerance genes in the periaqueductal gray (PAG), prefrontal cortex, temporal lobe and ventral striatum, respectively. A potential gene network was identified in the PAG in which 19 of the 34 genes were strongly associated with tolerance. Eleven of the network genes were found to reside in QTLs previously associated with morphine-related behaviors while 7 were predictive of tolerance (morphine-naïve condition). Overall, the genes modified by chronic morphine administration show a strong presence in canonical pathways representative of neuroadaptation. A potentially significant role for the microRNA and epigenetic mechanisms in response to chronic administration of pharmacologically-relevant doses of morphine was highlighted by candidate genes Dicer and H19.
The therapeutic value of opioids for pain relief is well established. Repeated opioid administration however, can lead to a progressive decline in analgesic efficacy. Extensive research into the mechanisms involved in analgesic tolerance has provided substantial insight into its neurobiology (Williams et al. 2001). However, despite recent advances, the cellular and molecular alterations downstream of the primary effects at the mu-opioid receptor (μOR) are complex and not well understood (Law et al., 2000).
One strategy to discover the mechanisms of morphine’s behavioral effects is to identify regions of the genome associated with these traits through linkage analysis. Behavior genetic investigations of morphine tolerance demonstrate a high degree of heritability (Liang et al., 2006; Mas et al., 2000, Kest et al., 2002). Behavioral Quantitative Trait Locus (bQTL) mapping assigns a significant portion of the genotype-dependent variance in morphine consumption and acute analgesia to the Mop2 locus on the proximal region of Chr10 (Bergeson et al., 2001; Belknap et al., 1995; Berrettini et al., 1994; Doyle et al., 2006). This is a promising candidate region since the gene encoding the μOR (Oprm) is located in Mop2. Reciprocal congenic genotypes in which this region has been introgressed from the opposing background (C57BL/6J (B6) or DBA/2J D2)) have been developed to further aid in the isolation of genetic factors (Ferraro et al., 2005; Bergeson et al., 2001). These studies have established a basic strategy with which to pursue specific genomic regions and candidate genes.
Large-scale gene expression profiling studies have also been utilized as a strategy to find candidate genes involved in tolerance. For instance, gene expression profiling in four opioid-naïve (Korotstynski et al., 2006) and morphine treated mouse strains have identified basal and morphine-induced gene expression differences in the striatum (Korotstynski et al., 2007) and nucleus accumbens (Grice et al., 2007) that may contribute to opioid-related phenotypes. However, none of the aforementioned studies anchored their gene expression results to within subjects derived behavioral data or investigated multiple brain regions especially to include those more directly involved in opioid-induced analgesia (i.e. periaqueductal gray (PAG)) and to assess region by strain and drug interactions. Characterizing several brain regions and anchoring expression data to the same experimental protocol and dosing schedule may improve the likelihood of focusing on genes directly relevant to the phenotype in question.
Reciprocal congenic genotypes have been an excellent tool for confirming regions of the genome hypothesized to be involved in a particular trait (Bergeson et al., 2001; Ferraro et al., 2005). The purpose of our study was to utilize reciprocal congenics to investigate the promising role of proximal Chr10 in nociception and the acute and chronic effects of pharmacologically relevant doses of morphine. A behavior genetics approach (Letwin et al., 2006) was employed to identify changes in gene expression directly related to morphine analgesic tolerance. In this manner morphine-induced changes in gene expression that are irrelevant to this trait were filtered from the pool of genes. Only those genes specifically correlated to the trait of interest were identified as candidate genes significantly influencing the development of tolerance.
Adult male B6 (C57BL/6J) and D2 (DBA/J) (Jackson Laboratories, Bar Harbor, ME), D2B6 (Mop2 locus from B6 introgressed into D2 background), and B6D2 (Mop2 locus from D2 introgressed into B6 background) mice 60–120 days old and weighing approximately 21–28 g at the start of the experiment were used. The D2B6 and B6D2 mice were shipped to the MPRC from the University of Pennsylvania where brother x sister mating continued to propagate the reciprocal strains. All animals were experimentally naive, housed in a temperature-controlled room (21° C; 0700–1900 lights on), and given free access to chow and tap water during the entire experimental procedure.
The reciprocal congenic genotypes were generated in the following manner. The portion of the chromosome thought to account for a significant portion of variance in response to morphine was reciprocally introgressed between the parental genotypes, B6 and D2. The introgress protocol began with a B6 × D2 intercross followed by backcross of the resulting F1 generation progeny to both parental genotypes. The backcross consisted of F1 males and females being mated with progenitor mice. The N2 backcross progeny were genotyped using DNA extracted from tail snips at the time of weaning. Mice were then selected for subsequent rounds of backcross breeding based upon heterozygosity at specific set of microsatellite DNA markers (B6D2: D10Mit75 (2.0cM) proximally to D10Mit124 (15.0cM) distally; D2B6: D10Mit75 (2.0cM) proximally to D10Mit61 (32.0cM) distally; Ferraro et al., 2005). Ten serial backcross generations were created in a reciprocal manner. At generation N10, heterozygous mice were intercrossed and DNA from the N10 F1 progeny was analyzed to identify homozygous congenic mice D2B6 and reciprocal congenic genotype B6D2 (Supplemental Figure 1).
Baseline sensitivity to pain and morphine-induced analgesia was measured by the hot-plate test (Harvard Hot-plate Analgesia Meter; Holliston, MA). The hot plate was maintained at a constant temperature of 55 °C for all test conditions. Latency to paw lick (front or hind) or jump was used as the dependent measure. Under all experimental conditions a cut-off time of four times each genotype’s saline control value was used to avoid tissue damage.
To determine the baseline sensitivity for each genotype, we examined the latencies of each genotype naïve to morphine. Each genotype received one injection of saline and was tested 15 minutes later on the hotplate (n = 6–9animals/genotype). Data were analyzed using a one-way ANOVA with genotype as a main effect.
The subject was placed on the hot-plate 15 min post s.c. injection of 0, 0.3, 1.0, 3.0, 10.0 or 30.0 mg/kg morphine to determine the relative potency of morphine to produce analgesia. Each dose for each genotype consisted of 6–9 animals and was tested for only one day. A two-way analysis of variance (ANOVA) was performed for dose-response curves across genotype. Data were analyzed as the percent of maximal analgesic response (%MPE) as determined by the following formula: 100 * [(test latency) - (saline baseline latency)] ÷ [(4 * baseline latency) - (baseline latency)]. The ED90 for each genotype was derived from the regression analysis of the linear portion of each dose-response curve.
Analgesic tolerance to morphine was determined in three distinct experiments. In each experiment the dose of morphine required to produce approximately 90% (ED90 as determined above) of the maximal possible analgesic effect (%MPE) was used to induce tolerance in each genotype. The ED90 dose of morphine for the DBA, D2B6, C57, and B6D2 mice was as follows: 7.7, 7.7, 31.2, and 31.2 mg/kg, respectively. Mice were tested on the hot-plate 15 min post injection. Saline or the ED90 dose of morphine was administered s.c. between the hours of 8:00–11:00am.
The three experimental conditions were as follows: 1) Three Injection. For each genotype, two groups of mice (Saline and Morphine, N = 12–16/group) were tested once every other day during a five day period (i.e. Monday, Wednesday, Friday). Brains were harvested for gene expression analysis only in this experimental condition, 2) Five injection, phenotypic confirmation. This test was performed as a complimentary means to confirm the results of the three-injection schedule and to push the dosing to determine if D2 mice would develop tolerance using a shorter injection interval and more injections. For each genotype, two groups of mice (Saline and Morphine, N = 6–9/group) were tested once every day during a five day period. 3) Learning Assessment. A subset of genotypes were run in a manner that would determine to what degree, if any, a genotypes response to morphine was influenced by repeated exposure to the hot plate. In this experiment, the ED90 dose of morphine was given in a manner exactly the same as the Three dose group (once every other day; Monday, Wednesday, Friday), however, in one group the hot plate was turned off until the final day of testing while the hot plate was turned on (55°C) for the other group as described above. Only B6D2 and B6 mice were run in this experiment since they were the two genotypes that develop tolerance to a degree in which a learning component, if present, might be discerned (n=6–9/group per genotype). The mice were tested on the hotplate 15 minutes post injection.
In the Three Injection and Five Injection Conditions, a one-way ANOVA was used to determine if equi-analgesic effects of morphine were produced on the first day of drug administration. A two-way (genotype X day) repeated measures ANOVA was used to determine genetic differences in the analgesic effects of morphine across repeated administration. Morphine-treated mice at each respective dosing interval were compared to saline-treated mice tested under the same procedure. Data are presented as %MPE. The Three Injection treatment was compared to the Five Injection phenotype confirmation group using a three-way repeated measures ANOVA (genotype X treatment paradigm X injection number). Two comparisons were made; one using the first 3 tests of each paradigm and second using the first, third and fifth tests of the Five Injection treatment paradigm (compared to the three tests of the Three Injection paradigm). In the Learning Assessment paradigm a one-way ANOVA (heated vs. room temperature hot plate) was conducted on the last test session to determine if learning occurred due to repeated exposure to the hot plate.
Brain tissue from the Three Injection tolerance study was harvested from all four genotypes and treatment conditions (n = 12–16 mice/group) 3 hours following the end of tolerance testing on the last test day. Each mouse was euthanized by CO2 asphyxiation followed by cervical dislocation and decapitation. A plexiglass brain mold (David Kopf Instruments, Tujunga, CA) was used to slice the fresh whole brain into coronal slices appropriate for each section desired (see below).
Using the atlas of Franklin and Paxinos (Paxinos and Franklin, 2001) as a reference the following sections were taken: Prefrontal cortex (PFC): Interaurally 5.90-3.70 mm, from the dorsal most point to 2.25 mm ventrally at a medial lateral width of 1 mm centered on midline. Ventral striatum (VS): Interaurally 5.90-3.70 mm, from a dorsal-ventral position of 3.75 mm ventrally to just short of the ventral surface (not including olfactory tubercle) and from a medial lateral position 2.0 mm to midline. This region included regions of interest such as the nucleus accumbens, ventral pallidum, bed nucleus of the stria terminalis and substantia innominata. Temporal lobe (TL): Interaurally 2.74-1.74 mm from the dorsal most point at 3.80 mm to the ventral surface from a medial point of 1.80 mm laterally to the lateral surface. Periaqueductal gray (PAG) containing dorsal raphe: Interaurally 0.00 – (−1.16) mm from the dorsal point 1.80 mm to 3.00 mm ventrally at a medial lateral width of 1.5 mm centered on midline.
The PAG has consistently been implicated in nociception, antinociception, and tolerance (Tortorici et al., 2001)). In addition, the VS (Manning et al., 1994), TL (Nandigama and Borszcz, 2003; Pavlovic et al., 1996) and PF (Hardy and Haigler, 1985) have also been shown to play a role in nociception and the analgesic effects of opioids (Chudler and Dong, 1995; Zubieta et al., 2001). In addition to nociception and analgesia these areas were of interest in terms of their role in conditioned reinforcement and the reinforcing properties of drugs of abuse (VS, PFC) (Everett et al., 1999; Shippenberg and Elmer, 1998; Wise, 1998).
Tissue from the four brain regions was harvested for gene expression analysis on a 27,648 element mouse cDNA microarray employing a common reference design (Letwin et al., 2006; Lee and Saeed, 2007). The microarray was comprised of the National Institute on Aging Ko set of 15,247 cDNA clones (Tanaka et al., 2000) and National Institute of Mental Health Brain Molecular Anatomy Project (BMAP) set of 11,136 cDNA clones, altogether representing ~21,000 distinct mouse genes. Each brain region from 2–3 subjects was pooled to obtain adequate tissue samples, thereby forming 3–5 independent replicates. Total RNA was prepared using Trizol reagent (Invitrogen, Carlsbad, CA) and the RNAeasy mini kit (Qiagen, Valencia, CA) as per manufacturers’ instructions. Each pooled RNA sample was further divided into two equal aliquots to enable a technical replication known as a dye-swap hybridization (Lee and Saeed, 2007). Thus, expression analysis of a brain region for a particular genotype and treatment was derived from 3–5 independent hybridizations and 3–5 corresponding dye-swap hybridizations. For statistical analysis, a biological replicate experiment was defined as the average of an independent hybridization and its corresponding dye-swap hybridization (i.e. 5 independent and 5 corresponding dye-swap hybridizations yield 5 biological replicates). A hybridization experiment consisted of Cy-5 labeled cDNA that was reverse transcribed from 15 micrograms of pooled total RNA and co-hybridized with Cy3-labeled cDNA synthesized from an equal amount of the Stratagene Universal Mouse Reference RNA. Dye swap hybridizations were performed by reversing the dyes for each of the RNA samples. Hybridizations were performed for 18–24 h at 42°C followed by washing in decreasing concentrations of SSC at room temperature and spun dry. Microarray image scanning, fluorescence intensity measurements, LOWESS data normalization, normalization across replicate experiments, experimental noise determination and cluster analysis were performed as described previously (Letwin et al., 2006; Teramoto et al., 2005).
Gene expression values from saline and morphine treated subjects in four genotypes in four brain regions can be defined as cgsr or mgsr, where c is saline control, m is morphine, g is gene, s is genotype and r is region. These data yield two base datasets: A) control gene expression without morphine intervention (cgsr), and B) gene expression following morphine administration (mgsr). From these two data sets we derive an additional data set C), the change in gene expression relative to the mean value of the genotype and region control expression (mgsr - xcgsr, or delta). Since the morphine expression data is important mostly in relation to the control levels, the delta dataset was emphasized over mgsr in future analyses and discussion. Statistical analysis for genotype dependent differences in datasets cgsr and delta was accomplished by an analysis of variance (ANOVA) with 10% FDR to control for type I error resulting from multiple comparisons (Reiner et al., 2003). The implementation of a 10% FDR corresponded to an effective alpha value of ≤ 0.00787 (data not shown). Subsequently, the significant genes were filtered for a minimum 1.5-fold change across genotypes based on ‘self versus self’ hybridizations in which aliquots of the same RNA sample were labeled separately with Cy3 and Cy5 dyes and co-hybridized onto a microarray. Results from these control hybridizations (data not shown) and previous published studies (Yang et al., 2002; Lee et al., 2007; Liang, et al., 2008) demonstrate that a Z-score of approximately 2 to 3 corresponds to a 1.4-fold ‘difference’, providing a metric for inherent assay noise. Significant genes displaying a ≤1.4-fold change in the strain comparisons exhibited a low validation success rate (~50%) as assessed by quantitative real time RT-PCR (data not shown), whereas genes above this cut-off displayed a validation rate (>90%) in line with the chosen FDR. We will refer to these statistically significant and filtered datasets as: 1) ‘genotype-signature genes’, cgsr genes that are differentially expressed under control conditions as a function of genotype in each brain region, and 2) ‘morphine-signature genes’, delta genes that are significantly different as a function of genotype in each brain region.
The phenotypic endpoint used to define tolerance was set as [100 – (%MPE on day 3)] for each genotype. This number reflects the diminution in analgesic efficacy by the end of the chronic treatment regiment and assigns a higher value to the genotypes showing the most tolerance by the end of the treatment. A two-tier analysis strategy was employed to define candidate genes associated with the tolerance phenotype. First, genes had to exhibit significant genotype-dependent expression differences (i.e. genotype-signature and morphine-signature genes). Second, the genotype-signature or morphine-tolerance genes must be significantly correlated (p<0.05) with the tolerance phenotype as defined by pattern matching/feature selection (Letwin et al., 2006; Pavlidis and Noble, 2001). The subset of genotype-signature genes that were significantly associated with tolerance are defined as candidate Predisposition genes. The subset of morphine-signature genes that were significantly associated with tolerance are defined as candidate Tolerance genes.
Biological themes associated with morphine-signature genes were identified by using the three gene ontology (GO) categories of molecular function, biological process, and cellular component in the Expression Analysis Systematic Explorer (EASE) application (Hosack et al., 2003), which is executable using TIGR Multi Experiment Viewer (TMEV; available at www.tigr.org/softlab).
Identification of candidate genes that are highly interconnected or over-represent a biological pathway may provide insight into molecular events that predispose or are causally related to the development of tolerance. In order to provide a functional framework to assist interpretation and narrow the candidate gene list to high value targets, gene networks and canonical pathways were explored using Ingenuity Pathway Analysis (IPA, Redwood City, CA), Pathway Studio (Rockville, MD), the GO Consortium (www.geneontology.org), GenMAPP (www.genmapp.org) and Kyoto Encyclopedia of Genes and Genomes (www.genome.jp/kegg).
Ingenuity Systems is an online-based bioinformatics resource that algorithmically identifies gene networks and canonical pathways based upon the gene lists. The gene lists uploaded into the database consisted of unique gene identifier (e.g. RefSeq, GenBank accession, Affymetrix ID), expression and correlation values. In order to identify gene networks, IPA rank orders the gene list in order of interconnectedness (triangle connectivity), then adds genes to the top focus gene based upon specific connectivity (how much a new gene’s connections overlap with the seed gene), then combines smaller networks through the addition of linker genes found either in the uploaded list or IPA’s Global Molecular Network. A p-value, the number of focus molecules (f) and top functional potential is given for each network. The p-value is defined as the probability of finding f or more Focus Genes in a set of n genes randomly selected from the Global Molecular Network database using a right-tailed Fisher’s exact test. All significant biological themes associated with these networks were subsequently cross-validated using Pathway Studio, GenMAPP, GO and KEGG pathways. In order to identify canonical pathways, IPA quantifies the likelihood that the number of uploaded genes that overlap with any one of 163 established canonical pathways occurs at a level greater than chance using a right-tailed Fisher’s exact test. The p-value and ratio value (number of molecules in a given pathway divided by total number of molecules that make up the pathway) are provided for each pathway.
Predisposition and Tolerance genes in the PAG, PFC, TL and VS were associated with expression QTLs (eQTLs) using the WebQTL module available at www.genenetwork.org/home.html. INIA Brain mRNA M430 (Jun06) RMA, VCU BXD PFC Sal M430 2.0 (Dec06) RMA, INIA Brain mRNA M430 (Jun06) RMA, HBP Rosen Striatum M430V2 (Apr05) RMA served as the databases for eQTL searches in the PAG, PFC, TL and VS, respectively. Significant linkage of eQTL markers to correlated genes was defined by marker regression plots with 1000 permutations. The results from the permutation tests provide likelihood ratio statistics (LRS) scores that are suggestive, significant or highly significant.
Expression levels of a portion of the genes were confirmed by real time RT-PCR on an ABI Prizm 7700 Sequence Detection System as previously described (Malek et al., 2002). Total RNA from PFC, TL, VS, and PAG in all genotypes was reverse transcribed using random primers following the manufacturer’s instructions. The resulting cDNA was diluted and used as a template for quantitative PCR. PCR primers were selected for specificity by the National Center for Biotechnology Information Basic Local Alignment Search Tool (NCBI BLAST) of the mouse genome, and amplicon specificity was verified by first-derivative melting curve analysis with the use of software provided by PerkinElmer (Emeryville, CA) and Applied Biosystems. Quantitation and normalization of relative gene expression were accomplished by using the comparative threshold cycle method as described previously (Joe et al., 2005). The following genes were identified by microarray analysis as differentially expressed across genotypes and subjected to RT-PCR validation: protein kinase C gamma (NM_011102), protein kinase A (NM_008854), cannabinoid receptor 1 (NM_007726), adenosine receptor 1 (NM_001008533), and dicer1 (NM_148948). The expression of the “housekeeping” genes carnitine acetyltransferase (NM_007760), peptidylprolyl isomerase-like 4 (NM_026141) and succinate-CoA ligase alpha subunit (NM_019879) were used for normalization as these genes did not exhibit differential expression in our microarray assays. Sequence of primers used for the real-time RT-PCR is given in the Supplemental Table 1.
There was a significant main effect of genotype (F(genotype) df3,47 = 8.67, p<.0001) on baseline nociception (Figure 1A). The B6 and B6D2 genotypes were significantly more sensitive than the D2 and D2B6 mice. The B6 and B6D2 genotypes did not differ from each other nor did the D2 and D2B6 genotypes. Therefore, baseline pain sensitivity was not affected by the introgression of the proximal region of Chr10.
There was a significant main effect of genotype on acute sensitivity to the analgesic effects of morphine (F(genotype) df3,23 = 35.40, p<.0001); the B6 and B6D2 genotypes were significantly less sensitive than the D2 and D2B6 mice (Figure 1B). ED50 values for morphine were 4.2 and 5.1 mg/kg for the D2 and D2B6 mice and 11.8 and 12.2 mg/kg for the D2 and D2B6 mice, respectively. As was true for baseline nociception, introgression of the proximal Chr10 region did not alter acute sensitivity to the analgesic effects of morphine; the B6 and B6D2 genotypes did not differ from each other nor did the D2 and D2B6 genotypes. This was a surprising finding considering bQTL analyses have mapped the proximal chromosome region of 10 to morphine analgesia (Belknap et al., 1995) and other morphine related-behaviors (Berrettini et al., 1994; Ferraro et al., 2005)
Figure 2A (Three Injection Experiment) shows the analgesic effects of morphine given every other day for five days as a percent of analgesia produced by the first injection. There was no significance difference in the analgesic effects of morphine following the first injection; therefore tolerance assessment was performed at approximately equi-analgesic doses. Genotype significantly influenced the development of tolerance (F(genotype) df3,50 = 17.25, p<.0005). Tolerance developed most rapidly in the B6 followed by the B6D2 strain. Tolerance did not develop in the D2 mice and only slightly in the D2B6 animals. Pair-wise comparisons for statistical significance following the last injection are given in the table insert.
Figure 2B (Five Injection, Phenotype Confirmation) shows the analgesic effects of morphine given everyday for five days as a percent of that produced by the first injection. Again, there was no significant difference in the analgesic effects of morphine following the first injection; therefore tolerance assessment was performed at approximately equi-analgesic doses. Genotype significantly influenced the development of tolerance ((F(genotype) df3,23 = 19.74, p<.0001). Tolerance developed most rapidly in the B6 followed by the B6D2 and D2B6 mice. Tolerance did not develop in the D2 mice. Pair-wise comparisons (Genotype x Day) for statistical significance following the last injection are given in the table insert. Overall, the shorter dose interval increased the rate of tolerance that developed after 3 injections (F(treatment schedule) df1,71 = 7.19, p<.009) as did the additional injections of morphine (F(treatment schedule) df1,71 = 14.27, p<.0003) (Figure 2A vs 2B). However, regardless of the test paradigm, the D2 genotype did not develop tolerance to morphine and there was no genotype by treatment schedule interaction (p<.31;ns); decreasing the inter-injection interval and increasing the number of injections produced a parallel change in the development of tolerance across genotypes. Thus, the Five Injection phenotype confirmation experiment confirms the results of the Three Injection schedule and extends the generality to schedules using a shorter injection interval and more injections.
The degree of tolerance was not influenced by learning factors associated with repeated hot-plate testing. There was no difference between the experimental groups exposed to a heated plate versus those exposed to room temperature plate in the B6D2 or B6 genotypes (p<0.86, p<0.34, respectively, both ns; Supplemental Figure 2).
The expression profiles of ~21,000 cgsr genes were measured in each of four brain regions (PAG, PFC, TL, VS) across four saline-treated strains (B6, D2, B6D2, D2B6), and a one-way ANOVA (main factor genotype) at 10% FDR was performed per brain region. Altogether, a total of 1070, 2550, 1758 and 2684 genotype-signature genes were identified (cgsr genes exhibiting significant strain-dependent differences in expression of at least 1.5-fold or greater) in the PAG, PFC, TL and VS, respectively (Figure 3A; Supplementary Table 2), representing 4562 unique genes across the 4 brain regions. In order to define genotype-dependent differences in morphine-induced changes in gene expression, the delta value for each gene was computed from the morphine and saline treatment data and subjected to a one-way ANOVA at 10% FDR per brain region. A total of 1182, 1642, 2041 and 2074 morphine-signature genes were detected in the PAG, PFC, TL, VS, respectively (Figure 3B, Supplementary Table 2), representing 3942 unique genes across the 4 brain regions. The number of unique morphine-signature genes is not unexpected given the well documented effects of morphine on the expression and/or activity of transcription factors such as ΔFosB and CREB (Zachariou et al., 2006; Han et al., 2006). The maximum delta/minimum delta ratio for each gene across the four genotypes ranged from 1.5 to 6, 1.5 to 8, 1.5 to 11 and 1.5 to 10-fold in the PAG, PFC, TL and VS, respectively (data not shown).
Morphine-signature genes (delta genes) were subjected to EASE analysis in order to identify potential biological themes associated with the four brain regions. Genes were classified by using GO categories belonging to the ontologies of molecular function, biological process, and cellular component. A total of 73, 166, 139 and 144 GO categories were significantly over-represented (Fisher’s exact test, p < 0.05) by morphine-signature genes in the PAG, PFC, TL and VS, respectively. Particularly noteworthy were the identification of epigenetic-related categories such as chromatin remodeling complex (GO:0016585) in the PFC, RNA interference (GO:0016246) and regulation of cell differentiation (GO:0045595) in the TL, and plasticity-related categories such as neuron development (GO:0048666) and cell projection organization and biogenesis (GO:0030030) in the VS.
An associative analysis was performed between signature genes and the behavioral endpoint tolerance for candidate gene discovery. Genotype-signature genes that correlate with tolerance are postulated to predispose animals to morphine-induced tolerance (such gene associations are termed ‘Predisposition genes’). A similar approach has been employed to identify “naïve” genes predisposing inbred mice to ethanol-related behaviors (Letwin et al., 2006). The tolerance value [100 – (%MPE on day 3)] for each genotype was used as the template to search for genotype-signature genes (Figure 3A); the expression value of a gene across the four genotypes was required to match the tolerance template across the four genotypes in order to be identified. This was done in each of the four brain regions. Employing pattern matching as a correlative search tool (also known as feature selection; Pavlidis and Noble, 2001), 64, 53, 86 and 123 Predisposition genes (positive and negative) were identified in the PAG, PFC, TL and VS, respectively (pattern match p-value < 0.05; Figure 4A; Supplemental Table 3). Tspyl4 and 6330407J23Rik were the only Predisposition genes physically located within the congenic region demarcated by markers D10Mit75 and D10Mit61. Using the GeneNetwork expression QTL (eQTL) resource (http://www.genenetwork.org), Tspyl4 was found to be regulated in trans by an eQTL locus associated with marker rs3682996 on chromosome 1 while 6330407J23Rik was determined to be regulated in cis (Supplemental Table 3).
Next, the tolerance value for each genotype was used as the template to search for ‘Tolerance genes’ (Figure 3B) within the morphine-signature dataset; again, the expression value of a gene across the four genotypes was required to match the tolerance template across the four genotypes in order to be identified as a candidate gene. Using feature selection, 81, 96, 106 and 82 Tolerance genes (positive and negative) were identified in the PAG, PFC, TL and VS, respectively (pattern match p-value < 0.05; Figure 4B; Supplemental Table 3). These genes are postulated to represent potential mediators of tolerance whose expression is modulated by morphine.
We also performed a two-way ANOVA (genotype and treatment serving as the main factors) with 10% FDR correction in each brain region and subjected the resulting significant genotype x treatment interactions to pattern matching in order to identify the tolerance genes. Despite the more conservative nature of the two-way ANOVA approach compared to the one-way ANOVA strategy used above, a substantial overlap of pattern matched genes was found between the two approaches. For example, 79% of the PFC tolerance genes originally identified by the one-way ANOVA approach were also captured in the two-way ANOVA approach. Overall, the overlap between the two approaches across the 4 brain regions ranged from 77% in the PAG to 80% in the VS. Given the paucity of information pertaining to genes associated with tolerance behavior, all subsequent analyses will be based on the results from the one-way ANOVA.
A phenotype-gene association (PGA) clustergram was generated by hierarchical clustering of the 369 Tolerance genes across the four brain regions (Figure 5). This plot is analogous to hierarchical clustering of gene expression values typically seen in microarray experiments (Eisen et al., 1998); however, rather than clustering one-dimensional expression values, PGA analysis clusters biological information pertaining to the direction and significance level of gene-phenotype correlations (Reiner-Benaim et al., 2007). Tolerance genes exhibiting significant positive correlations are indicated in red, significant negative correlations in green and non-significant correlations in black. Noteworthy in Figure 5 is the paucity of genes that are correlated in more than one brain region (<2%). A similar finding was noted for the predisposition set (<2%; see Supplemental Table 3). These observations suggest that networks comprised of different combinations of predisposition and/or tolerance genes in each brain region may be working in concert to promote tolerance susceptibility and acquisition. Tolerance genes could be grouped into a number of interesting categories such as gene expression, 2nd messenger mediated signaling, neuronal synaptic plasticity, response to oxidative stress, signal transduction/G-protein coupled receptor and transport (Figure 5).
The mapping of correlated genes to behavioral intervals provides additional supportive evidence when associating a gene to a behavior. Physical map locations of our correlated genes were obtained from ENSEMBL (http://www.ensembl.org) and morphine-related behavioral QTLs (bQTL) information was acquired from the literature. Provisional bQTLs mapping to chromosomes 1, 6 and 10 for morphine preference (Berrettini et al., 1994), chromosomes 3, 7, 8, 9, 10, 12 and 18 for morphine consumption (Belknap and Crabbe, 1992), chromosomes 2, 3, 4, 7, 9, 10, 11 and 19 for morphine analgesia (Belknap and Crabbe, 1992; Belknap et al., 1995; Bergeson et al., 2001), and chromosomes 1, 5 and 10 for morphine withdrawal (Kest et al., 2004) have been identified through crosses of the B6 and D2 strains, representing ~15% of the mouse genome. For our analysis, a bQTL interval was defined as two LOD units around the linkage peak (Rapp, 2000), provided that information on multiple genetic markers in the vicinity of the peak was given (see Berrettini et al., 1994). Otherwise, the interval was defined by 10–15 megabases on either side of a lone genetic marker defining the bQTL. Of the 369 Tolerance genes, 316 had unambiguous mapping information and of these mapped genes 73 (23%) were found to reside in bQTLs that influence morphine-related behaviors (Supplemental Table 3). Similar results were obtained in Predisposition gene set where 70 out of 258 mapped genes (27%) were found to reside in bQTLs (Supplemental Table 3). In order to ascertain whether or not the mapping of correlated genes in bQTLs may have occurred by chance, we compared the number of correlated and non-correlated genes residing in bQTLs for each brain region. There was a significant occurrence not due by chance alone of tolerance genes mapping to bQTLs for the temporal lobe (χ2, P=0.0069) and periaquaductal gray (χ2, P=0.0216); and predisposition genes mapping to bQTLs for the ventral striatum (χ2, P=0.0059) and prefrontal cortex (χ2, P=0.0221). A number of interesting Tolerance genes residing in bQTLs associated with morphine analgesia included the glycine receptor beta subunit (positively correlated in PAG), signal transducer and activator of transcription 5 (positively correlated in PFC), ADP-ribosylation factor interacting protein 1 (negatively correlated in TL) and tumor suppressor candidate 2 (negatively correlated in VS).
Of the Tolerance and Predisposition genes residing inside bQTLs but outside of the congenic region, 13 and 6 genes, respectively, were regulated in trans by eQTLs located in the congenic region (Supplemental Table 3). Moreover, 39 Tolerance and 25 Predisposition genes residing outside bQTLs were found to be regulated in trans by eQTLs in the congenic region (Supplemental Table 3). There was a significant difference in the occurrence of tolerance versus non-correlated genes mapping inside eQTLs for the temporal lobe (χ2, P= 0.0469), ventral striatum (χ2, P=0.0381) and periaquaductal gray (χ2, P=0.0053); and predisposition versus non-correlated genes for the temporal lobe (χ2, P= 0.0392) and periaquaductal gray (χ2, P=0.0032). Notable trans-regulated examples included synaptic vesicle membrane proteins (synaptotagmin XI, syntaxin 3), transcription factors (Stat5a) and extracellular matrix/cell adhesion molecules (laminin B1, integrin beta 3 binding protein).
Promoter analysis of the 3 network types using the TRANSFAC motif database and MAST (Malek et al., 2006, Lee et al., 2007) did not identify any statistically significant over-represented transcription factor binding motifs (data not shown). To address whether the correlated genes are potentially interlinked into a gene regulatory network, all Predisposition and Tolerance genes were computationally analyzed using microarray bioinformatics resources available through Ingenuity Pathway Analysis, Pathway Studio, the GO Consortium (www.geneontology.org), GenMAPP (www.genmapp.org) and Kyoto Encyclopedia of Genes and Genomes (www.genome.jp/kegg). Three network ‘types’ can be generated in a region-specific manner: a network of genes based solely from the predisposition set, a network derived solely from the tolerance set, and a network derived from the combined consideration of predisposition and tolerance genes. The top 3 networks in each network type for each brain region are provided in Supplemental Table 4. Each network type affords a slightly different viewpoint. The predisposition network may interact in such a way as to provide a genotype-dependent starting ‘set point’ that predisposes the subject to tolerance. These genes may not necessarily change expression in response to morphine but set the stage for morphine’s effects on other genes. The tolerance networks may interact in such a way as to cause (or prevent via compensatory mechanisms) tolerance. These genes do not necessarily differ at baseline rather they change expression in response to morphine administration in a genotype-dependent manner. The combined consideration networks define potential preexisting interactions of predisposition and tolerance genes during the acquisition of tolerance.
Figure 6, Panel B, shows an example network derived from the Tolerance gene list in the PAG. Nineteen of the 34 genes in the PAG tolerance network were either positively (12) or negatively (7) associated with the development of tolerance. Eleven of the Tolerance genes reside in bQTLs previously associated with morphine-related behaviors (see Supplemental Table 3). The probability of having 11 of 19 network tolerance genes residing in morphine-related bQTLs was significantly greater than would be expected by chance alone (χ2, P<0.0001). Despite having less than 2% overlap overall between the Predisposition and Tolerance genes (across all brain regions), 7 genes in the PGA tolerance network had corresponding predisposition gene counterparts. As a means to consider the potential alterations in network associations that may occur following chronic drug treatment, Figure 6, Panel A, shows a simplified network comprised of the 7 Predisposition gene counterparts (H19, Ttk, Hba2, B4galt1, Hmox1, Casp3, cFos). A comparison of Panels A and B was made in order to demonstrate how prevailing gene networks may undergo alterations in network associations following chronic drug treatment. Some Predisposition gene counterparts undergo a reversal in direction of their association with tolerance following morphine treatment (H19, Ttk, Hba2, B4galt1). For instance, Ttk was positively correlated in saline treated animals but negatively correlated in animals injected with morphine. Hba2 in the same network was negatively correlated in saline treated animals but oppositely correlated upon treatment with morphine. Overall, the genes modified by chronic morphine administration show a strong presence in canonical pathways representative of neuroadaptation (see next section) and illustrate potential alterations in network associations following chronic drug administration.
A complementary analysis to the de novo network construction using candidate genes is to determine their representation in molecular pathways with known function. Supplemental Table 5 presents the top 15 canonical pathways associated with genotype-signature, morphine-signature, Predisposition, and Tolerance genes in each brain region along with the number of genes found in each pathway and its corresponding p-value. As noted previously, the pathways identified using the genotype- and morphine-signature lists show a strong presence in canonical pathways representative of neuroadaptation. In this regard synaptic long-term potentiation was identified in the top 5 genotype- and morphine-signature lists in the PAG and TL, and PAG and PFC (Supplemental Figure 3). Axonal guidance signaling was identified in the top 5 canonical pathways in all brain regions within genotype-signature list and in 4 of the 5 brain regions in the morphine-signature list. Additional pathways within the morphine-signature list suggesting a major role of synaptic plasticity include the ephrin receptor, Huntington’s disease and neuregulin signaling. The smaller number of genes in the Predisposition and Tolerance lists limit the statistical probability that they will populate a canonical pathway. Nevertheless, there were several pathways in the Tolerance list suggestive of neural plasticity (neuregulin, notch and ERK/MAPK signaling) and cell death (NRF-2 mediated oxidative stress response, GM-CSF and p53 signaling).
Five Tolerance genes (protein kinase C gamma, protein kinase A alpha subunit, cannabinoid receptor 1, adenosine A1 receptor, dicer-1 protein) involved in 6 correlations were selected for microarray validation. These genes exhibited the smallest fold-change (1.5 to 2-fold) of maximum delta/minimum delta across the four strains, presumably representing the lower limits of our analysis to define significant correlations. Protein kinase C gamma, protein kinase A catalytic alpha subunit, cannabinoid receptor 1 and adenosine A1 receptor have been implicated in morphine-related behaviors including tolerance, albeit these studies are typically performed in a single mouse strain (Kaplan and Leitte-Morris, 1997; Narita et al., 2001; Smith et al., 2007; Trang, Sutak, and Jhamandas, 2007). It should be noted that none of the aforementioned genes had, prior to this study, been correlated with tolerance across multiple genetic strains. The fifth gene chosen for validation, dicer-1, has not previously been linked to drug tolerance. Finally, protein kinase A catalytic alpha subunit and adenosine A1 receptor genes were of interest because they mapped to bQTLs associated with morphine consumption and preference, respectively (Supplemental Table 3). Quantitative real time RT-PCR validation of the five genes in six correlations is shown in Figure 7. A significant correlation was found between the microarray and real time RT-PCR results (p<0.05) and between the RT-PCR results and tolerance behavior for all genes, supporting the reliability of our gene expression measurements.
One of our goals was to investigate the role of the proximal region of Chr10 in the acute and chronic effects of morphine. Previous bQTL studies suggest that this region between markers D10Mit28-D10Mit3 (4–21cM) accounts for a significant amount of variance in morphine two-bottle choice preference and analgesia (Berrettini et al., 1994; Belknap et al., 1995; Bergeson et al., 2001)(Supplemental Figure 1). A comparison of the B6 and D2 strains in the Mouse Genome Database (http://www.informatics.jax.org) reveals the presence of 21 non-synonymous coding and 207 non-coding single nucleotide polymorphisms (SNPs) between D10Mit28-D10Mit3, which were thought to account for the observed variance in behavioral responses. Oprm, the gene responsible for morphine’s primary receptor target (μOR), is located within this region at 8.0cM. Reciprocal introgression of this promising region however did not alter baseline nociception or the acute analgesic effects of morphine in our study. These results are in contrast to those found by Bergeson et al., (2001) using a similar congenic approach in which the Chr10 region between markers D10Mit28 and D10Mit3 (4–21 cM) from the donor B6 strain was introgressed into the D2 background (Supplemental Figure 1). In that case, acute sensitivity to 10 mg/kg morphine (i.p.) was altered in the reciprocal congenic strains; baseline nociception was not reported nor was tolerance assessed. Differences in stimulus intensity, route of administration, endpoint or alterations in baseline sensitivity (Elmer et al., 1998) could contribute to the discrepancy. Equally important is the fact that the introgressed regions differ substantially in size between our congenics and the ones studied by Bergeson et al. (2001). Specifically, our D2B6 strain carries an additional 11 cM portion or 37.4 Mb of Chr10 derived from the B6 strain (Supplemental Figure 1B). Consequently, genetic interactions (i.e. networks) likely diverge between the two D2B6 lines due to the presence of 162 SNPs residing within 73 genes (many encoding transcription factors) found in this 11 cM region. In fact, 16 of the 162 SNPs have been identified as eQTLs predicted to trans-regulate a subset of our Tolerance and Predisposition genes. The consequences of different introgressed genetic footprints are analogous to the ‘congenic footprint’ phenomenon which can lead to profound differences in gene expression and behavioral responses (Schalkwyk et al., 2007; Lee et al., 2007). Also noteworthy is our finding that while reciprocal introgression was without effect on acute morphine analgesia, presumably due to the genetic footprint, profound effects were observed for chronic morphine analgesia. While differences in environmental context and protocol may ultimately explain a large part of the differences between studies the difference in introgressed footprint may fortuitously have explanatory significance worth pursuing.
Tolerance to the analgesic effects of morphine develops at a rate directly related to the magnitude of the initial pharmacological effect (Cox and Tiffany, 1997). If the same dose of morphine is used to induce tolerance in subjects that differ significantly in morphine’s acute analgesic potency, tolerance could appear to develop more rapidly in sensitive genotypes given the greater magnitude of the initial pharmacological effect. This confound would affect the tolerance outcome and interpretation. The current study accounts for this by employing genotype-specific doses within a ‘therapeutic’ range (ED90). Under these conditions, reciprocal introgression of proximal Chr10 altered the development of tolerance. Thus, the acute analgesic effects of morphine are determined in part by genetically separable mechanisms from those that are involved following chronic administration. This conclusion is supported by previous behavior genetic (Oliverio and Casellano, 1974; Kest et al., 2002; Liang et al., 2006) molecular and pharmacological investigations (Ocana et al., 08; Kaneto et al, 1985). A potential confound to this approach is that the genotype-specific ED90 dose is sufficient to produce analgesia acutely in all genotypes but does not reach a threshold for inducing the molecular mechanisms involved in tolerance in all genotypes. If the morphine dose was increased sufficiently to produce tolerance in the D2 mice (i.e. ~ten-fold dose increase produces tolerance in D2; Kest et al., 2002) a different set of genes may be revealed to produce tolerance. However, this approach also runs the risk of confound by non-specific neuropharmacological mechanisms engaged at supramaximal dosing.
Given the strong bQTL evidence, a straightforward explanation for differential response to chronic morphine in the B6 and D2 mice would be found in the variations in Oprm expression or protein regulation. In this regard, B6 and D2 mice differ qualitatively in VS μOR regulation following chronic high dose morphine administration; μOR concentration decreases (−50%) in the B6 mice and increases (+ 20%) in the D2 mice (Petruzzi et al., 1997). There are two promoter regions in the Oprm gene, one upstream (5’ of exon 11) and one downstream (5’ of exon 1) (Liang and Carr, 1997; Pan, 2002). Five SNPs have been identified within a 1.7 kilobase (kb) region of the Oprm downstream promoter region between the B6 and D2 strains; three were near or within potential transcription binding sites (Doyle et al., 2006). A difference in the basal Oprm promoter activity was seen in vitro between the B6 and D2 constructs suggesting a role of the polymorphisms in gene transcription (Doyle et al., 2006). However, our current study and those utilizing B6 and D2 constructs in μOR-positive BE(2) cell lines (Doyle et al., 2006) suggest other factors are involved in tolerance besides Oprm expression since μOR was not differentially expressed in the cell lines nor did we find differential expression of the μOR or differential expression of genes known to promote or repress Oprm (Law et al., 2004; see also Contet et al., 2008). Given the complex and often-times controversial evidence surrounding μOR regulation in tolerance, systematic large-scale gene expression analysis may yield alternative molecular mechanisms not found by investigating prototypical candidate mechanisms.
Large-scale gene expression studies were performed in conjunction with the behavioral studies to identify candidate genes associated with tolerance. A stumbling block in expression studies has been the lack of a systematic method to iteratively focus on the most promising candidates. We employed canonical pathway and network analyses to help identify high value candidate predisposition and tolerance genes. Canonical pathway analysis of the morphine-responsive genes suggests a significant role for neuroadaptive processes in response to chronic morphine administration (i.e. LTP, axonal guidance, ephrin, neuregulin pathways; see Supplemental Figure 3 for LTP example). Network analysis was used to identify highly interconnected genes in order to provide insight into molecular events that predispose or are causally related to the development of tolerance. It is possible to analyze and view gene networks from several different vantage points. One viewpoint would be to determine network associations comprised of tolerance genes (Figure 6, right panel) and determine the correlation status of the same genes prior to morphine administration (Figure 6, left panel). An alternative permutation would consist of predisposition genes (Supplemental Figure 3, left panel and their status after chronic morphine administration (Supplemental Figure 3, right panel). Changes in the correlation status of genes transitioning from baseline to morphine treatment in the network permutations may be of significant interest. This is nicely exemplified by the top network involving tolerance genes in the PAG. Nineteen of the 34 genes in this network were either positively (12) or negatively (7) associated with the development of tolerance. Protein kinase C, cannabinoid receptor 1, and calmodulin are known to play a prominent role in the current theories of morphine tolerance and neuroadaptation ((Law et al., 2000; Williams et al., 2001; Trang et al., 2007; McClung and Nestler, 2008). Out of the 19 tolerance genes contributing to this network, more than half (11) remarkably resided in bQTLs previously associated with morphine-related behaviors. The occurrence of these correlated genes in bQTLs was significantly greater than would be expected by chance alone, presumably reflecting their importance in tolerance acquisition. Prior to morphine administration, 7 of the 34 genes making up this network were found to be associated with tolerance. Three of these genes (c-Fos, Casp, Pld) were not correlated after tolerance developed which could be traced to a parallel morphine response (down-regulation) across the genotypes. Strikingly, 4 genes (B4galt1, Hba2, Ttk, H19) exhibited a reversal in the direction of correlation following chronic morphine administration. For example, high H19 expression is predictive of tolerance development (direction of predisposition association) and a negative change in expression (i.e. down-regulation) is associated with tolerance (direction of tolerance association). H19 expression appears to decrease in the genotypes that become tolerant to morphine and increase in those that do not.
Several recent studies provide intriguing insight into the potential role of H19 in response to chronic morphine administration. First, it has been suggested that H19 is a primary microRNA precursor (to miR-675) (Cai and Cullen, 2007). Coincidently, our data reveals that Dicer1, encoding the pre-microRNA processing enzyme (Bernstein et al., 2001), is associated with the development of tolerance in the PFC, and EASE analysis of morphine-responsive genes identified the RNA interference pathway as playing a role in the TL. Second, analysis of 12,000 genes and ESTs between control undifferentiated and de-differentiation derived cell clones establish only one gene as differentially regulated, H19 (Scott et al., 2005). In a separate study, Ayesh et al., (2002) identified 48 differentially expressed genes modulated by H19 in a manner that may promote metastasis; interestingly we found that morphine administration altered the expression of 50% of these genes but not all were directly associated with the development of tolerance. Third, the H19 imprinting region has been implicated in dopamine neuron differentiation (Freed et al., 2008) and gene regulation via interchromosomal interactions involving the CCCTC-binding factor (Ling et al., 2006). Interestingly, our findings reveal the CCCTC-binding factor to be a genotype-signature gene as well. Moreover, five dopamine-related genes (dopamine D3 and D4 receptor, dopamine transporter, DARPP-32 and calcyon) were all found to be differentially regulated by morphine in this study. While none of these genes were found to be correlated with analgesic tolerance, they may be involved in traits such as the rewarding effects of morphine. Overall, our findings suggest that a form of transdifferentiation or metaplasia-like activity may be an important component in the neuroadaptive response to chronic morphine. Thus, while this imprinting region has been recognized as important in neurodevelopment and cellular differentiation, our study suggests that it may play a significant role in response to chronic administration of pharmacologically-relevant doses of morphine and warrants future investigation.
In summary, epigenetic and miRNA mechanisms and neuroadaptation appear to play a significant role in the nervous systems response to chronic morphine administration as corroborated by converging lines of evidence from EASE/GO, network and canonical pathway analyses. Specific components of the neuroadaptive response are likely related to specific behavioral endpoints.
Funded by NIDA/NIH Grant# DA015087 (GIE, NHL). We would also like to acknowledge Wade H. Berrettini, MD, Ph.D. and Thomas N. Ferraro, Ph.D. for their valuable contribution of mice and consultation.