Introduction

Methamphetamine (Meth) addiction is a major public health concern and currently there are no approved pharmacological treatments [1]. Relapse during abstinence is a key challenge for treating Meth addiction [2]. In rats, drug seeking across several drug classes, including Meth, progressively increases after prolonged withdrawal from extended access to drug self-administration [3]. This incubation phenomenon has also been translated to a clinical setting. Wang et al. [4] showed that cue-induced drug craving in Meth-dependent individuals increased for up to 3 months of abstinence. These clinical data highlight the translational value of neural mechanistic studies of this incubation phenomenon in rats.

In previous studies, we and others have identified several mechanisms underlying incubation of Meth craving in rats, including: a critical role of dopamine signaling [5] and histone deacetylase 5 [6] in dorsal striatum; projections from lateral anterior intralaminar nucleus of thalamus to dorsomedial striatum; [7] and calcium-permeable AMPA receptors in nucleus accumbens (NAc) [8]. However, whether neural mechanisms at the transcriptional level [9] contribute to this incubation is unknown. During our initial exploration, we used a candidate-gene approach and demonstrated that incubated Meth seeking is associated with selective gene regulation in activated Fos-expressing dorsal striatal neurons [5].

In the current study, we used a genome-wide approach and characterized the transcriptional profiles of central amygdala (CeA) and orbitofrontal cortex (OFC) during incubation of Meth craving. We chose CeA because we demonstrated that reversible inactivation of CeA by GABAA + GABAB receptor agonists (muscimol + baclofen) decreases incubated Meth seeking after 1-month withdrawal, but has no effect on non-incubated Meth seeking after 2-day withdrawal [10]. The critical role of CeA neural activity in incubation of craving generalizes across different drug classes (cocaine [11, 12], morphine [13], and nicotine [14]), as well as a natural reward (sucrose [15]). Thus, characterization of the CeA transcription profile during incubation of Meth craving may uncover common transcriptional regulation underlying incubation of craving.

We chose OFC based on its critical role in relapse to cocaine [16] and heroin seeking [17], as well as its important role in the cognitive processes (e.g., decision making, learning, and memory) related to drug addiction [18, 19]. Additionally, several clinical studies have demonstrated disrupted OFC function in Meth-dependent individuals [20,21,22]. By contrast, we recently showed that intra-OFC injection of muscimol + baclofen in rats has no effect on incubated Meth seeking [10]. Thus, characterization of OFC transcription profiles may reveal transcriptional mechanisms potentially involved in drug seeking beyond acute neuronal activity in OFC during the relapse/incubation tests.

We used RNA-sequencing (RNA-seq) to perform genome-wide transcriptomic profiling on CeA and OFC tissue collected after 2 days or 35 days of withdrawal from extended access of Meth (or saline, as the control condition) self-administration training. Our findings demonstrate that CeA and OFC exhibited distinct transcriptional patterns during incubation of Meth craving, and provide novel transcriptional targets for future mechanistic studies of the incubation phenomenon.

Material and methods

Subjects. See Supplementary Online Material

Intravenous surgery. See Supplementary Online Material

Apparatus. See Supplementary Online Material

Drug self-administration training. See Supplementary Online Material

Withdrawal phase. See Supplementary Online Material

RNA isolation. See Supplementary Online Material

Library preparation and RNA-seq. See Supplementary Online Material

cDNA synthesis and qPCR. See Supplementary Online Material

Exp. 1 Genome-wide transcriptional profiling of CeA during incubation of Meth craving

As shown in Fig. 1a, we performed intravenous surgeries on two groups of rats (total n = 26) and trained them to self-administer either saline (n = 12) or Meth (n = 14) as described in Supplementary Online Material in 2 independent runs. We performed live decapitation on withdrawal days 2 and 35, and collected CeA tissue (as in Fig. 1c) for mRNA preparation. We used the extracted mRNA for library preparation and RNA-seq. We pooled tissue from two rats as one biological replicate. The number of biological replicates in each group was: Day 2: Saline = 3, Meth = 4; Day 35: Saline = 3, Meth = 3.

Fig. 1
figure 1

Tissue collection for subsequent RNA-seq analysis. a Experimental timeline. b Meth and saline self-administration training. Data are mean ± SEM number of Meth (0.1 mg/kg/infusion) infusions and active and inactive lever presses during the ten 9-h daily self-administration sessions for Exp. 1 (total n = 26). During training, active lever presses were reinforced on an FR1 20-s timeout reinforcement schedule, and Meth or saline infusions were paired with a 5-s tone-light cue. c Diagram for orbitofrontal cortex (OFC) and central amygdala (CeA) tissue collection

Full size image

Exp. 2 Genome-wide transcriptional profiling of OFC during incubation of Meth craving

As above, we performed intravenous surgeries on two groups of rats (total n = 32) and trained them to self-administer saline (n = 16) or Meth (n = 16) in two independent runs. We performed live decapitation on withdrawal days 2 and 35, and collected OFC tissue (as in Fig. 1c) for mRNA preparation. Note: the rats in Exp. 1 and 2 are from different cohorts. We used the extracted mRNA either for library preparation and RNA-seq or for cDNA synthesis and qPCR. We pooled tissue from two rats as one biological replicate. The number of biological replicates in each group was: Day 2: Saline = 4, Meth = 4; Day 35: Saline = 4, Meth = 4.

Statistical analyses of behavioral and qPCR data

See Supplemental Online Material

Statistical and bioinformatic data analysis of RNA-seq data

Differential expression analyses: We aligned the sequencing reads to the rat rn4 reference transcriptome with Tophat, and counted the reads against Ensembl Rnor 6.0 annotation with HTSeq Python package [23]. Next, we normalized raw counts of each gene based on both the sequencing depth and mRNA length to generate reads per kilo-base per million mapped reads (RPKM). All the RNA-seq processing steps were done in parallel using an in-house pipeline called SPEctRA [24]. Subsequently, we performed pair-wise differential expression comparisons using Voom Limma on the filtered gene lists, which met the criterion that the majority of the samples within at least one experimental condition have RPKM of 1.0 or greater [25]. We set the significance threshold of nominal p < 0.05 and log (Fold Change) > 1.3. We performed analyses of gene expression between the saline and Meth groups on either withdrawal day 2 or day 35 separately, because of differences observed in the gene expression patterns of the Saline day 2 and Saline day 35 groups. The reasons for these differences are unknown, but they may be related to either age differences (day 35 rats were older than day 2 rats) or the fact that the brains of these two groups were extracted on different days.

Enrichment analyses: We used the Fisher’s exact test (FET) to evaluate the overlap of differential expression lists between withdrawal day 2 and day 35 [26]. For identifying gene ontology terminology (GO term), we performed PANTHER overrepresentation test (Version 13.1) [27] on the differentially expressed genes (DEGs) in CeA on withdrawal day 35 and on the overlapped genes identified by Rank–rank hypergeometric overlap (RRHO) in OFC (upregulated on day 2 and downregulated on day 35). The number of DEGs in other groups was too small to observe statistically significant overrepresentation.

RRHO analyses: See Supplementary Online Material

Upstream regulator analysis: See Supplementary Online Material

Motif analysis: See Supplementary Online Material

Accession numbers

Raw and processed RNA-seq gene expression data are available via the Gene Expression Omnibus (GEO) database (GSE111243). Gene lists used by HOMER are also available via GEO (GSE21512).

Results

Meth self-administration training

Meth-trained but not saline-trained rats increased their number of infusions over the course of training (Group × Day interaction, F9,216 = 8.2, p < 0.01) (Fig 1b). Moreover, active lever presses, but not inactive lever presses, increased throughout training in the Meth-trained but not saline-trained rats (Lever × Day × Group interaction, F9,216 = 8.8, p < 0.01). The training data presented here were from Exp. 1. The training data from Exp. 2 have been published previously [5] and were similar to the training data of Exp.1.

Exp. 1: Genome-wide transcriptional profiling of CeA during incubation of Meth craving

The goal of Exp.1 was to evaluate transcriptome-wide changes in gene expression in CeA associated with incubation of Meth craving. We performed RNA-seq of CeA after 2-day or 35-day withdrawal from Meth self-administration. We used rats that underwent saline self-administration as the control groups for each timepoint.

We analyzed the data from day 2 and day 35 independently and used the between-subjects factor of group (Meth, saline). We found there was a much greater number of DEGs on withdrawal day 35 when compared with day 2 (METH vs. SAL—Day 2: 215 genes; Day 35: 2217 genes; Table S2). The union heatmap shows highly divergent patterns of DEG induction (Fig. 2a). The Venn diagrams confirmed minimal overlap between upregulated or downregulated DEGs at the two withdrawal timepoints (upregulated: 4 genes; downregulated: 3 genes; Fig. 2b). FET analyses confirmed that there was no significant overlap of DEG lists between the two timepoints (Fig. 2c). In contrast, RRHO analyses showed that there was a significant signal in the upper-left quadrant (Day-2-UP vs. Day-35-DOWN). The overlap was centralized on “day-2-axis” but not on “day-35-axis” (Fig. 2e), indicating that a subset of upregulated genes on day 2 returned to baseline on day 35. However, 1255 out of 1258 genes were not considered DEGs when applying the threshold, suggesting that this phenomenon only occurred in genes not robustly regulated on day 2. Based on these observations, we did not follow up on this pattern.

Fig. 2
figure 2

Greater regulation of genes in central amygdala (CeA) during incubation of Meth craving. a Heatmaps show the union of withdrawal day 2 (Meth vs. Saline) and day 35 (Meth vs. Saline) DEGs in CeA rank ordered by log2(fold change) of withdrawal day 2. b Venn diagrams represent the number of the common and unique downregulated (left panel) and upregulated (right panel) DEGs on withdrawal day 2 and withdrawal day 35 in CeA. c Table shows the odds ratio, Jaccard Index and p-value for Fisher’s exact test for enrichment of DEGs between withdrawal day 2 and withdrawal day 35. N.S., nonsignificant. d Five of the gene ontology terms significantly enriched in downregulated (top) and upregulated (bottom) DEGs on withdrawal day 35 in CeA. Red line indicates p = 0.05. e RRHO map compares all expressed genes between withdrawal day 2 and withdrawal day 35. Each pixel represents the overlap between downregulated and upregulated genes in a transcriptome-wide threshold-free manner on the two withdrawal days. The significance of overlap (-log10(p-value) of a hypergeometric test) is color coded, with warmer colors indicated increasing significance. f Transcription factor motifs enriched in downregulated (left) and upregulated (right) DEGs on withdrawal day 35 in CeA

Full size image

GO-term analysis in CeA on withdrawal day 35 showed that the significantly overrepresented biological processes in the downregulated DEGs included: NADH dehydrogenase complex assembly (fold enrichment (FE) = 23.32, FDR = 7.14 × 10−7); Mitochondrion organization (FE = 5.11, FDR = 4.09 × 10−5); ATP metabolic process (FE = 7.24, FDR = 2.59 × 10−4); Translation (FE = 4.11, FDR = 3.35 × 10−4); and RNA processing (FE = 2.76, FDR = 0.0155). The significantly overrepresented biological processes in the upregulated DEGs included: Protein ubiquitination (FE = 2.41, FDR = 1.15 × 10−8); Regulation of dendritic spine morphogenesis (FE = 5.67, FDR = 9.01 × 10−5); Regulation of post-synapse organization (FE = 4.03, FDR = 0.00115); Activation of mitogen-activated protein kinase (MAPK) (FE = 2.64, FDR = 0.0318); and Histone lysine methylation (FE = 3.26, FDR = 0.0379) (Fig. 2d, Table S3).

We performed IPA analysis and found that on day 35 the majority of upstream regulators were miRNAs, all of which were negatively enriched (Table S4). This prediction is likely driven by robust upregulation of those DEGs, which are potential miRNA downstream targets. HOMER analysis on DEGs from withdrawal day 35 showed 22 enriched transcription factor motifs in the downregulated DEGs and 30 in the upregulated DEGs. Of these, 3 are common: interferon regulatory factor 4 (IRF4), myeloid zinc finger 1 (MZF1), and YY1 transcription factor (YY1) (Fig. 2f).

In summary, the data in Exp. 1 demonstrate that there were 10-fold more DEGs in CeA on withdrawal day 35 than day 2, indicating a positive association between transcriptional regulation in CeA and incubation of Meth craving. Among these DEGs, there were significant enrichment of genes involved in several relevant biological processes. Furthermore, there were several upstream regulators predicted for each group of DEGs, most of which are unique to the direction of gene regulation observed.

Exp. 2: Genome-wide transcriptional profiling of OFC during incubation of Meth craving

The goal of Exp. 2 was to evaluate the transcriptome in OFC associated with incubation of Meth craving. We performed RNA-seq of OFC on withdrawal days 2 and 35 after Meth self-administration. We also performed qPCR on three genes to validate the RNA-seq data.

We analyzed the data from day 2 and day 35 independently and used the between-subjects factor of group (Meth, saline). We found a greater number of DEGs on withdrawal day 2 compared with withdrawal day 35 (METH vs. SAL–Day 2: 118 genes; Day 35: 55 genes; Table S2). Like the CeA, heatmaps of the two timepoints demonstrate visually dissimilar patterns of differential expression (Fig. 3a), and Venn diagrams show minimal overlap of upregulated or downregulated DEGs at the two withdrawal timepoints (upregulated: 1 gene; downregulated: 0 genes; Fig. 3b), which was further confirmed by FET (Fig. 3c).

Fig. 3
figure 3

Opposing regulation of genes in orbitofrontal cortex (OFC) during incubation of Meth craving. a Heatmaps show the union of withdrawal day 2 (Meth vs. Saline) and day 35 (Meth vs. Saline) DEGs in OFC rank ordered by log2(fold change) of withdrawal day 2. b Venn diagrams represent the number of the common and unique downregulated (left panel) and upregulated (right panel) DEGs on withdrawal day 2 and withdrawal day 35 in OFC. c Table shows the odds ratio, Jaccard Index, and p-value for Fisher’s exact test for enrichment of DEGs between withdrawal day 2 and withdrawal day 35. N.S. nonsignificant. *p < 0.05. d RRHO map compares all expressed genes between withdrawal day 2 and withdrawal day 35. Each pixel represents the overlap between downregulated and upregulated genes in a transcriptome-wide threshold-free manner on the two withdrawal days. The significance of overlap (-log10(p-value) of a hypergeometric test) is color coded, with warmer colors indicated increasing significance. Upregulated-Day 2 and Downregulated-Day 35 showed a robust overlap and we performed further analysis in these overlapping genes (e-h). e Five of the gene ontology terms significantly enriched in the overlapping genes from RRHO analysis. Red line indicates p = 0.05. f Transcription factor motifs enriched in the overlapping DEGs from FET analysis. g Upstream regulators predicted in the OFC lists with activation scores indicated for withdrawal days 2 and 35. Negative activation score indicates predicted inhibition of the regulator, and positive activation score indicates predicted activation. h Partial network of known interactions from transcription factor motifs and upstream regulators (IPA). Solid line—direct interaction; dashed line—indirect interaction; arrow—observed activating interaction; blunted line—observed inhibitory interaction. Panel is colorized based on direction of interactions and expression in the Day 2 expression list. Orange—predicted activation; Blue—predicted inhibition; Yellow—increased expression. The top row are upstream regulators predicted by IPA, second row are intermediary regulators, third row are HOMER-predicted transcription factors, and the bottom row are DEGs that are overlapping in the Day 2 Upregulated and Day 35 Downregulated lists

Full size image

Interestingly, we observed a significant overlap between upregulated DEGs on withdrawal day 2 and downregulated DEGs on withdrawal day 35 (Day-2-UP vs. Day-35-DOWN: Odds Ratio = 46.457, p = 4 × 10−9; Fig. 3c). RRHO analysis showed a significant overlap in genes with opposite expression patterns on day 2 and 35, similar to the results from the FET analysis (Fig. 3d). In addition to the large overlap between Day-2-UP and Day-35-DOWN, we also observed a weak overlap between Day-2-DOWN and Day-35-UP. Based on these findings, we focused on the overlapping genes between Day-2-UP and Day-35 DOWN for GO term (RRHO list) and HOMER analyses (DEGs).

GO-term analysis showed several significantly overrepresented biological processes including: sensory perception of chemical stimuli (FE = < 0.01 (negatively enriched), FDR = 1.36 × 10−11); Regulation of cell adhesion (FE = 2.74, FDR = 3.97 × 10−6); Regulation of cell communication (FE = 1.6, FDR = 3.4 × 10−5); G-protein coupled receptor signaling (FE = 0.46 (negatively enriched), FDR = 3.17 × 10−3); and Neurogenesis (FE = 1.65, FDR = 0.00671) (Fig. 3e, Table S5).

HOMER analysis showed 26 predicted transcription factor motifs (Fig. 3f). Notably, one of these (IRF4) was also predicted in CeA day 35 analyses (Fig. 2f). IPA revealed several predicted upstream regulators, the top 20 of which are listed in Fig. 3g with the direction of activation noted (Table S6).

We used IPA software to create a partial network of known interactions that connect a representative list of IPA-predicted upstream regulators to some of the HOMER-predicted transcription regulators. We also included direct and indirect connections from some of these predicted regulators to 2 representative DEGs that belong to the list of overlapping genes between Day-2-UP and Day-35-DOWN – Dcn: decorin and Col3a1: collagen alpha-1(III) chain (Fig. 3h).

We created visualizations of approximate locations of HOMER-predicted transcription factor motifs on Dcn and Col3a1 using known binding data from MatInspector (Fig. 4a). We show the representative expression patterns of Dcn and Col3a1 in Fig. 4b. According to the sequencing dataset, these two genes are increased in the Meth group on withdrawal day 2 and decreased on withdrawal day 35 when compared to the respective saline groups.

Fig. 4
figure 4

Validation of RNA-seq by qPCR. a Approximate locations of predicted transcription motifs on Dcn and Col3a1 (MatInspector). b Wiggle plots of Dcn and Col3a1 from RNA-seq. Each respective gene is illustrated below its expression track (blue). Red line indicates the greatest number of reads for each track. c Dcn, Col3a1, and Henmt1 mRNA expression in OFC. Data are presented as fold change of mean values in the saline group. Error bars indicate SEM. *p < 0.05; n = 4 per group

Full size image

We validated the expression pattern of three genes by qPCR based on RNA-seq data—Dcn, Col3a1, and Henmt1 (a small RNA 2′-O-methyltransferase not changed in RNA-seq) (Fig. 4c). We analyzed the data from days 2 and 35 independently and used the between-subjects factor of group (Meth, saline). Analyses of qPCR data showed that expression of Dcn and Col3a1 mRNA in OFC is increased on withdrawal day 2 (Dcn: t6 = 4.1, p = 0.006; Col3a1: t6 = 2.8, p = 0.030) and decreased on withdrawal day 35 (Dcn: t6 = 2.5, p = 0.048; Col3a1: t6 = 2.5, p = 0.045). Lastly, there was no change of Henmt1 mRNA expression in OFC on either withdrawal day 2 or 35 (p > 0.05). These results by qPCR are consistent with RNA-seq data.

In summary, the data in Exp. 2 demonstrate that there were more DEGs in OFC on withdrawal day 2 than day 35, and a significant overlap between upregulated DEGs on withdrawal day 2 and downregulated DEGs on withdrawal day 35. RRHO further confirmed overlap in a transcriptome-wide threshold-free manner. Among the genes found to overlap by RRHO, there was significant enrichment of several relevant biological processes. Furthermore, predicted upstream regulators of the overlapped DEGs were found to be connected with some predicted transcription factor motifs in known pathways. Two DEGs, validated by qPCR, also contain binding motifs of many of the predicted transcriptional regulators.

A final question is whether DEGs overlap across brain regions. We applied FET and found that there was no significant overlap of DEGs in the two brain regions at either timepoint (Fig. 5a). RRHO analyses showed that there was no overlap between CeA and OFC on withdrawal day 2, and a minimal transcriptome overlap, but centered on the least significantly changed genes, in the same direction on day 35 (Fig. 5b). Taken together, these data suggest that there was no or minimal overlap of DEGs between CeA and OFC.

Fig. 5
figure 5

No overlapping genes between CeA and OFC during incubation of Meth craving. a Table shows the odds ratio, Jaccard Indexm, and p-value for Fisher’s exact test for enrichment of DEGs between CeA and OFC on withdrawal day 2 and withdrawal day 35, respectively. N.S., nonsignificant. b RRHO map compares all expressed genes between CeA and OFC on withdrawal day 2 (left) and 35 (right). Each pixel represents the overlap between downregulated and upregulated genes in a transcriptome-wide threshold-free manner on the two withdrawal days. The significance of overlap (-log10(p-value) of a hypergeometric test) is color coded, with warmer colors indicated increasing significance

Full size image

Discussion

We used RNA-seq and performed transcriptional profiling of two brain regions, CeA and OFC, during incubation of Meth craving. We generated a unique source of genome-wide gene expression data (publicly available in GSE111243) in CeA and OFC after 2-day and 35-day withdrawal from extended access Meth self-administration training in rats. In CeA, we identified a positive association between the number of DEGs and incubation of Meth craving, but did not observe significant overlap of genes among the DEGs lists. While GO-term analysis of downregulated genes in CeA on withdrawal day 35 showed an enrichment of genes mainly involved in energy and metabolism, GO-term analysis of upregulated genes showed an enrichment of genes involved in protein ubiquitination, spine morphology, synaptic function, MAPK activity, and histone methylation. Additionally, we predicted upstream regulators for each DEG list and found that most of these upstream regulators are unique to the direction of regulation.

In OFC, we identified a negative association between the numbers of DEGs in OFC and incubation of Meth craving, with a significant overlap between upregulated DEGs on withdrawal day 2 and downregulated DEGs on day 35. Among these overlapped genes, there was a significant enrichment of genes involved in biological processes such as sensory perception of chemical stimuli and regulation of cell adhesion. Furthermore, we predicted upstream regulators and transcription factor motifs of the overlapped genes and found that they are connected partially in known pathways. These findings not only provide new insight into the transcriptional mechanisms during incubation of Meth craving, but also provide a direct comparison between two discrete brain regions previously examined in a reversible inactivation study of this incubation [10].

Distinct transcriptional regulation of CeA and OFC during incubation of Meth craving

In the reversible inactivation study mentioned above, inactivation of the CeA but not OFC decreased incubation of Meth craving. Similarly, our data suggest that at the transcriptional level, the CeA plays a more important role in relapse to Meth seeking than the OFC. First, there were more overall DEGs in CeA than in OFC on both withdrawal day 2 (~2 fold more) and day 35 (~20 fold more). Second, the greater number of DEGs in CeA on withdrawal day 35 than day 2 suggests that transcriptional mechanisms in CeA display a form of “incubation” as seen behaviorally. In contrast, the greater number of DEGs in OFC on withdrawal day 2 than day 35 suggests that OFC may play a role in non-incubated Meth seeking during early withdrawal.

We identified a greater number of upregulated DEGs than downregulated DEGs in OFC on withdrawal day 2, indicating a general activation of gene transcription during early withdrawal. Additionally, there was a significant overlap between upregulated DEGs on withdrawal day 2 and downregulated DEGs on day 35, which was absent in CeA. This analysis suggests that several upregulated DEGs on day 2 in OFC underwent a time-dependent suppression after prolonged withdrawal. One interpretation is that upregulated genes in OFC might suppress Meth seeking during early withdrawal, and loss of these genes leads to the elevated Meth seeking after prolonged withdrawal. The transcriptional profile in OFC also correlates well with previous functional imaging studies showing that OFC is hypermetabolic during early withdrawal and hypometabolic after prolonged withdrawal in cocaine [28, 29] or Meth [30,31,32,33] users.

In contrast, in CeA, we identified a greater number of upregulated DEGs than downregulated DEGs on both withdrawal day 2 and 35, but there was no overlap among the two sets of DEGs. These data suggest a general activation of transcription but that distinct sets of genes are regulated at the two different timepoints of CeA during incubation of Meth craving. Additional analyses suggest that the negatively enriched miRNAs underline the massive gene upregulation in CeA on withdrawal day 35, which identify a direction for future studies.

Overall, these patterns of data are consistent with our previous finding showing that reversible inactivation of CeA by GABA receptor agonists decreases incubated Meth seeking, but have no effect on Meth seeking during early withdrawal [10]. Additionally, global OFC inactivation has no effect on incubated Meth seeking [10], but a potential role of OFC in Meth seeking during early withdrawal is supported by a recent clinical study demonstrating an association between Meth seeking and deficits in dopamine signaling in lateral OFC of recently abstinent Meth users [20]. However, whether OFC plays a role in Meth seeking during early withdrawal in rats has not been examined.

Finally, we found no or minimal overlapping DEGs between CeA and OFC during incubation of Meth craving. However, direct comparison of DEGs from CeA and OFC should be made with caution here, because of the different numbers of biological replicates and independent cohorts for tissue collection for CeA and OFC.

Candidate-gene targets in CeA during incubation of Meth craving

Our data also provide candidate-gene targets for future studies on the role of CeA transcriptional mechanisms in incubated Meth seeking after prolonged withdrawal. For example, upregulated DEGs were significantly enriched in the MAPK signaling pathway in CeA on withdrawal day 35. This finding suggests that this pathway in CeA, which has been implicated in incubation of cocaine [11] and opiate [13] craving, also contribute to incubation of Meth craving. Another interesting biological process enriched among the upregulated DEGs in CeA on withdrawal day 35 is protein ubiquitination. The ubiquitin-proteasome system (UPS), which mediates protein turnover, has been previously implicated in synaptic plasticity and learning and memory [34,35,36]. Additionally, UPS function in the NAc has been shown to play a role in memories associated with the rewarding effects of cocaine [37, 38] and morphine [39], as assessed by conditioned place preference (CPP) and self-administration procedures.

Most notably, expression of several E3 ubiquitin ligases, which catalyze ubiquitin transfer to a specific substrate before its degradation in proteasomes [40], was selectively elevated in CeA on withdrawal day 35 in the Meth but not Saline group. This finding is novel because very few studies to our knowledge have determined the role of E3 ubiquitin ligases in the behavioral effects of addictive drugs [41, 42]. Jiao et al. [41] reported that expression of Synovial apoptosis inhibitor 1, an E3 ubiquitin ligase, in dorsal striatum increases after Meth CPP, which leads to degradation of the GABAA1 receptor subunit. Caffino, et al. [42] reported that exposure to cocaine during adolescence is associated with increased ubiquitin-protein ligase E3A in medial prefrontal cortex, which leads to Arc protein degradation in adulthood.

Among the E3 ubiquitin ligases observed in our study, two of them are of interest. One is March7, which has been implicated in protein degradation after activating D1-receptor signaling in vitro [43]. This is of potential relevance because our recent study demonstrates a critical role of D1-receptor signaling in CeA in relapse to Meth seeking after food choice-induced voluntary abstinence [44]. The other is Rnf2: Kumari et al. [45] reported that loss of Rnf2 in cultured hippocampal neurons contributes to generation of silent synapses, which within the NAc play a critical role in incubation of cocaine craving [46].

Additionally, we found that expression of several histone lysine methyltransferases in CeA were elevated on withdrawal day 35. This finding extended previous studies which demonstrated extensive epigenetic adaptations after Meth exposure in cortical and striatal areas [47], as well as the role of striatal histone methylation and acetylation in Meth-associated memories (assessed by CPP procedure) [48] and incubated Meth seeking [6], respectively. One of the candidate histone lysine methyltransferases identified here is Kmt2a (also known as Mll1), which is also elevated in activated dorsal striatal neurons during incubated Meth seeking [5]. Knockdown of Kmt2a in NAc has been shown to decrease Meth CPP [48].

A final note is that we could not validate our RNA-seq data from CeA by qPCR, due to unexpected mRNA degradation during sample storage. However, validation of RNA-seq data from OFC by qPCR (Fig. 4) demonstrates the general reliability of the DEGs reported in our study. Nonetheless, the candidate genes identified here need further validation both at the mRNA and protein levels, as well as at the behavioral level in future studies.

Conclusions

Our genome-wide study has identified distinct transcriptional profiles in CeA and OFC during incubation of Meth craving. Our findings highlight the CeA as a key site for transcriptional regulation in incubated Meth seeking, which is consistent with our previous anatomical mapping study showing that CeA plays a critical role in incubation of Meth craving [10] and relapse to Meth seeking after food choice-induced voluntary abstinence [44]. In contrast, transcriptional mechanisms in OFC may be associated with Meth seeking during early withdrawal. Finally, the present data reveal numerous candidate-gene targets and potential upstream regulators for examining the transcriptional mechanisms in CeA in relapse to Meth seeking, as well as relapse to other addictive drugs.