Adult male (10-week-old) C57BL/6 J mice (purchased from HFK Bioscience Co., Ltd, Beijing, China) were used in this study. Animals were habituated for at least 1 week in the new environment before experimentation. They were housed 5 per cage, on a 12 h light-dark cycle, and in a constant temperature environment (21 ± 1°C), with access to food and water ad libitum.

All mice were randomly divided into cages, with five mice per cage. Cocaine was dissolved in 0.9% saline. Repeated-cocaine treatments is a popular paradigm to explore the mechanisms of cocaine addiction. Consecutive intraperitoneal injections for 7 or 14 days are commonly used in previous studies, which had been demonstrated neuroadaptational changes associated with cocaine addition (Feng et al., 2014; Puig et al., 2014). In the present study, animals received intraperitoneal injections of cocaine (Qinghai Pharmacoceutical, China) at a dose of 20 mg/kg/d during 9:00 to 11:00 a.m. for 14 consecutive days. All control groups received saline injections for 14 consecutive days. Mice were asphyxiated using dry ice, and then the bilateral PFC (0–1.5 mm from the tip of brain) was excised from each animal individually at 2 h, 24 h, or 7 days after the final injection of cocaine and then pooled randomly. The cortical subregions in the excised tissue included frontal association cortex, prelimbic cortex and orbital cortex. To account for inter-animal variations, we obtained 2 biological replicates for each treatment group, with each replicate representing the PFCs pooled from 5 animals. The tissue was stored in liquid nitrogen immediately and then transferred to −80°C refrigerator. All animal protocols were approved by the Review Board of the Institute of Psychology, Chinese Academy of Sciences and were performed strictly in accordance with the Guideline for Care and Use of Laboratory Animals of the Chinese Academy of Sciences.

To characterize the transcriptome of the PFC, we used RNA sequencing (RNA-seq) to measure the expression levels of all polyA-containing transcripts in the PFC of mice treated chronically with cocaine or saline (control).

RNA was extracted from PFC tissues using Trizol (Invitrogen, Cat. 15596026) following the manufacturer's instructions. Briefly, 50–100 mg tissue was homogenized in 1 ml of TRIzol using a homogenizer. RNA was separated by adding 0.2 ml of chloroform per 1 ml of TRIzol reagent and then precipitated with isopropanol. RNA quantified using a NanoDrop 2000 (Thermo-Fisher Scientific, Waltham, MA, USA), and assessed with an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). The RNA integrity number (RIN) of each sample was >8.

Total RNA (240–380 ng/ul, 2 ug) was used to prepare cDNA libraries according to the manufacturer's protocol (TruSeq™ RNA Sample Prep Kit v2 -Set A, Illumina, San Diego, CA, USA). The brief procedure of library preparation is as follows: enrichment of mRNA by poly A, mRNA fragmentation, first strand synthesis, second strand synthesis, end repair, A tailing, adaptor ligation and PCR (13 cycle, 22.4–26.6 ng/ul × 30 ul).

Paired-end 100 nt sequencing was performed using the Illumina HiSeq2000 (TruSeq PE Cluster Kit v3–cBot - HS and TruSeq SBS Kit v3–HS). These 12 raw sequencing data sets were deposited in the Gene Expression Omnibus of the NCBI (https://www.ncbi.nlm.nih.gov/geo) under accession number GSE89572.

After removal of the adaptor sequences and trimming of low quality sequences using cutadapt (Martin, 2011), all reads were aligned to the UCSC mm9 version of the mouse genome assembly using the TopHat2 short read alignment program(Trapnell et al., 2009). Only uniquely mapped reads were used in this analysis. Fragments Per Kilobase of sequence per Million mapped reads (FPKM) were used to normalize the number of aligned reads by the size of the gene and the total number of mapped reads. The changes in gene expression between different groups were analyzed using the Cufflinks 2.0 package (Trapnell et al., 2012). Cufflinks can estimate the relative abundances of these transcripts based on how many reads support each one, taking into account biases in library preparation protocols (http://cole-trapnell-lab.github.io/cufflinks/manual/). Difference in expression of genes with FDR (False Discovery Rate) value < 0.05 were considered as differentially expressed genes.

We identified statistically enriched functional terms from the public databases of the KEGG Pathways (http://www.genome.jp/kegg/pathway.html). Functional terms with FDR-value < 0.05 was considered as statistically significant enrichment.

Time-series clustering was performed to identify statistically significant expression patterns and gene groups among differentially expressed genes (DEGs) using the short time-series expression miner (STEM; Ernst and Bar-Joseph, 2006). STEM uses a method of analysis that takes advantage of the number of genes being large and the number of time points being few to identify statistically significant temporal expression profiles and the genes associated with these profiles. User-defined parameters were set to 35 profiles, maximizing the number of clustering genes, with a minimum correlation threshold of 0.7. Five statistically significant clusters were identified (p < 0.05).

cDNA was synthesized from 1 μg of total RNA (template) using Superscript III (Invitrogen, CA, USA). The validation of DEGs were determined by real-time quantitative PCR using SYBR Green Mix (CWBIO, Beijing, China). PCR conditions were: 95°C for 10 min, 40 cycles at 95°C for 10 s, 60°C for 60 s, and signal detection for 10 s; dissociation curve analysis was performed at the end of each run, and relative mRNA expression was determined by the standard ΔΔCt method using glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as a normalization control. The primers are listed in Table S1.

Statistical analyses were performed with GraphPad Prism 6.0 software. Means of relative gene expression between two groups were compared with a two-tailed Student's t-test. Statistical significance was set at p < 0.05.

We obtained 36, 913, 048 to 50, 325, 962 short reads of 100 bp from each replicate from both the cocaine treated and control animals. Of these, 84–87.08% were successfully aligned to a reference gene database (Ensembl: Mus musculus, NCBI) by TopHat (Trapnell et al., 2009).

For our initial analysis, we used stringent false discovery rate (FDR) cutoffs of < 0.05. As showed in Figure 1, we obtained 463, 14, and 535 DEGs from samples obtained 2 h, 24 h, and 7 days after the final administration of cocaine, respectively (corrected p < 0.05). Among of them, there were 62/401 (2 h), 3/11 (24 h), and 311/224 (7 d) DEGs with decreased/increased expression level, respectively (Tables S2, S3 and S4). Only 4 DEGs were shared by all 3 time points, namely Dusp6, Fos, Nr4a1, and Egr1 (Figure 1).

Through functional enrichments of DEGs, several signaling pathways were enriched at each of the time points (Table 1). Two categories of signaling pathways were enriched (corrected p < 0.05) at 2 h. The first category comprised the signaling pathways regulating cell division, cell differentiation, and cell death, e.g., the Wnt signaling pathway, mitogen activated protein kinase (MAPK) signaling pathway, and Jak-STAT signaling pathway. The second enriched pathway was related to the circadian rhythm. Even though the number of DEGs was greatly decreased after 24 h of withdrawal, the circadian rhythm and MAPK signaling pathway were significantly enriched. At the 7 days post-withdrawal time-point, the number of DEGs was increased; however, most of these genes were significantly enriched (FDR < 0.01) in metabolic pathways, for example, the electron transport chain, ribosomes, proteasomes, citrate cycle, and RNA polymerase. These results indicate that distinct pathways are involved in the acute, sub-acute, and long-lasting stages of withdrawal after chronic cocaine treatment.

To further utilize the dynamic expression information, we performed time-series correlation analysis on the total DEGs of all 3 time points using the STEM method, which produced 5 clusters of DEGs with significant serial correlation of expression changes (r > 0.7; p < 0.05). Among them, the expression levels of DEGs in Cluster 1 gradually decreased with time, while those in Cluster 2 first increased and then decreased or returned to the baseline level (the same level as those in the saline group). Expression changes in Cluster 3 and Cluster 4 shared a similar pattern of fluctuation, with an initial increase followed by a decrease and then finally another increase. Whereas, expression levels of DEGs in Cluster 5 initially decreased, then increased, and decreased again at the end (Figure 2).

There were 3 sub-clusters in Cluster 1 that were consistent with the general tendency of gradual or fluctuant decrease, as shown in Figures 2A, 3A. The most highly enriched pathways in Cluster 1 were oxidative phosphorylation, and neurodegenerative disorders such as Huntington's disease, Parkinson's disease, and Alzheimer's disease. In addition, the oxidative phosphorylation chain and proteasome were significantly enriched in both Clusters 1 and 5 (Figures 2E, 3E).

There were 3 sub-clusters in Cluster 2 as showed in Figures 2B, 3B. KEGG pathway analysis of genes within Cluster 2 showed that the MAPK signaling pathway, the p53 pathway, TGF-β signaling pathway, insulin-signaling pathway, and circadian pathway were significantly enriched.

Clusters 3 and 4 (Figures 2C,D, 3C,D) shared a similar pattern of changes in the expression of DEGs, an initial increase followed by a decrease and then finally an increase. KEGG pathway analysis of Clusters 3 showed that they were mainly enriched in plasticity-related pathways, such as long-term potentiation and the regulation of the actin cytoskeleton, while axon guidance was enriched in both Clusters 3 and 4.

We validated these RNA-seq data independently from tissue samples using qRT-PCR. As shown in Table 2, we identified the expression of 10 genes at 3 different withdrawal time points, with 83.3% of genes/time points showing consistent results. Some of validated genes have been widely studied on their effects on cocaine exposure, while others have known functions, which may contribute to drug addiction. Fosb encoded ΔFosB, a truncated and stable protein, which is a well-characterized transcription factor induced by chronic exposure to virtually all drugs of abuse (Wallace et al., 2008). Nr4a1 increases in frontal cortex (Freeman et al., 2002) and Nr4a3 increases in the medial caudate putamen and cingulate cortex (Werme et al., 2000) after chronic cocaine administration. Nr4a2 plays an important role in development and differentiation of midbrain DA neurons and Nr4a2 gene expression decreased in human post-mortem midbrain of cocaine abusers (Leo et al., 2007). Our results further demonstrated that Nr4a1 persistently increased while Nr4a2 and Nr4a3 showed transient changes after cocaine withdrawal, indicating their different roles in different stages of cocaine addiction. The roles of Wnt pathway in the NAc in long-term cocaine-induced neuroplasticity have been reported (Cuesta et al., 2017). Notch signaling is an important regulator of neuronal homeostasis and Notch1 activation is involved in cocaine-mediated regulation of PDGF-B expression (Yao et al., 2011). In the present study, we found both Wnt4 and Notch2 expression increased 2 h after withdrawal and their roles in cocaine addiction should be further investigated. Cocaine-induced histone methylation at many specific candidate genes is already known (11). We found Kdm6b, a histone demethylase, increased 7d after cocaine withdrawal, indicating its potential role in cocaine addiction for the first time.

Mitochondrial energy production, which occurs through oxidative phosphorylation, requires the action of various respiratory complexes together known as the electron transport chain enzymes (ETC). There are 5 ETC complexes: (i) Complex I (NADH dehydrogenase), (ii) Complex II (succinate dehydrogenase), (iii) Complex III (cytochrome bc1 complex), (iv) Complex IV [cytochrome c oxidase (COX)], and (v) Complex V (ATP synthase). Our results found that over 40 genes in the PFC involved in each ETC complex were gradual or fluctuant decrease after cocaine-withdrawal, indicating that energy metabolism impairment induced by repeated-cocaine treatment showed a time-dependent increase following withdrawal. Studies in humans, nonhuman primates, and mice have demonstrated cocaine-associated reductions in metabolic activity in the brain. By testing the regional cerebral metabolic rate for glucose, it was found that 10 or 14 days of repeated daily treatment decreased metabolic activity in the mesolimbic dopaminergic pathway, including the PFC (Hammer and Cooke, 1994). The glucose metabolism rate and blood flow significantly decreased in the PFC of cocaine addicts, indicating lower metabolism and activity in brain of cocaine abusers (Goldstein and Volkow, 2002; Volkow et al., 2004). At the molecular level, both metabolic proteins in cocaine addicts and mitochondrial proteins in animal model of cocaine abuse are decreased (Lull et al., 2008). Our results further showed a dynamic metabolic decrease after cocaine withdrawal and suggest that this gradual decrease could represent an adaption of the metabolic response to cocaine in the PFC. A recent study found that compulsive drug-seeking behavior induced by repeated cocaine treatment resulted from reduced pyramidal cell excitability mediated by PFC metabolic dysfunction (Chen et al., 2013). Here, we revealed changes of the energy metabolism network underlying repeated cocaine treatment and believe that an adaption of the metabolic response to cocaine in the PFC may be another important mechanism of relapse during protracted abstinence.

Dynamically decreased genes were also enriched in the proteasome pathway, including the subunits composing the core particle of proteasomes during ubiquitin-mediated proteolysis. The ubiquitin proteasome system (UPS) controls the degradation of misfolded, newly synthesized proteins, as well as the turnover of specific target proteins. Consistent with our results, many transcriptomic and proteomic studies have described the effects of treatment with drugs of abuse on proteasome subunits or proteins involved in the ubiquitination process (Massaly et al., 2014). In the case of cocaine treatment, ubiquitin-conjugating enzyme E2N and proteasome subunit-alpha type 2 were involved in cocaine-conditioned place preference (Guan and Guan, 2013). Interestingly, it is suggested that synaptic and signaling proteins, especially transcription factors, which mediate enduring neuronal plasticity, could be regulated directly or indirectly by a UPS-dependent process (Massaly et al., 2014). Combined with lasting changes in transcription factors such as Nr4a1, Fos, and Egr1 in our results, we speculate that UPS-dependent regulation plays an important role in neuronal plasticity associated with the long-term effect of cocaine.

The MAPK pathway is one of the most conserved signaling pathways in eukaryotes. Multiple proteins involved in this pathway transfer extracellular biological, chemical, or physical signals into the cell by adding or removing adjacent phosphate groups of the protein molecule(Roux and Blenis, 2004). The MAPK signal pathway is divided into 3 sub-pathways: the ERK (extracellular signal-regulated kinase), JNK (c-Jun NH₂-terminal kinase), and p38 signal pathways (Roux and Blenis, 2004). Four genes, phosphorylase Dusp6, and transcription factors Egr1, Fos, and Nr4a1, which all belong to the ERK signaling pathway, showed a persistent expression change after cocaine withdrawal (Table 3).

Several lines of evidence have demonstrated the important roles of the ERK pathway in cocaine addiction. Both acute and repeated injections of cocaine activated ERK in the mesolimbic reward system, including the VTA, PFC, bed nucleus of the stria terminalis, and amygdala (Lu et al., 2006). In addition, ERK contributes to various behavioral effects of cocaine, such as psychomotor sensitization, conditioned place preference, and reconsolidation of memories for cocaine cues (Lu et al., 2006). Importantly, long-lasting changes of the ERK pathway induced by repeated cocaine treatment have been supported by the association of long-term ERK-mediated effects and ERK-dependent activation of FOS and EGR1 (Radwanska et al., 2005). High rates of relapse to cocaine use in humans after prolonged abstinence indicated that the drug-induced adaptation in cellular mechanisms was persistent during abstinence (Lu et al., 2004). A time-dependent increase in cocaine seeking induced by cocaine cues involves activation of ERK signaling, suggesting this pathway mediates the incubation of cocaine craving (Lu et al., 2005). Together, we suggest that long-lasting changes of the ERK signaling pathway in the PFC induced by cocaine treatment may be one of the most important mechanisms underlying cocaine addiction.

In the present study, we found that long-term potentiation (LTP), the regulation of the actin cytoskeleton, and the axon guidance pathway demonstrated fluctuant changes after withdrawal of repeated cocaine treatment, consistent with the current opinion that during abstinence, neuroplasticity at many levels, ranging from cellular alterations to morphological and structural changes, plays an essential role in addictive behaviors (Kasanetz et al., 2012; Nyberg, 2014). Repeated cocaine exposure modifies the induction of tetanus-induced LTP in vitro and in vivo, which was associated with dysfunction of glutamate receptors (Yao et al., 2004; Goto and Grace, 2005). The actin cytoskeleton has been shown to promote presynaptic vesicle movement, postsynaptic glutamate receptor trafficking, and morphogenesis of dendritic spines (Nyberg, 2014). In addition, changes in dendritic spine diameter and density due to chronic cocaine administration were associated with a deteriorating actin cytoskeleton and reduction in glutamate signaling-related proteins (Shen et al., 2009). However, alterations of drug-induced synaptic plasticity were guided by molecular cues, such as axon guidance molecules, that form projection-target connections during rearrangement of the synapse (Bahi and Dreyer, 2005). Together with these results, our findings indicate that potential associations among these varied levels of neuroplasticity play an important role in cortical reorganization after withdrawal of repeated cocaine treatment.

Previous studies suggest that circadian genes, including Per1, Per2, and Npas2 in NAc and hippocampus play roles in modulating cocaine reward or sensitization (Abarca et al., 2002; Imbesi et al., 2009). In dorsal striatum, >20% (29 genes total) were identified as genes associated with the circadian system in a total of 139 genes, which were regulated by cocaine self-administration (Lynch et al., 2008). Together with our findings of circadian pathway enrichment in the PFC after repeated cocaine treatment, it is indicated that circadian genes play essential roles in cocaine addiction regardless of treatment paradigm or brain regions.

As mentioned above, MAPK pathway also plays important roles in cocaine addiction across different brain regions. There is evidence indicating that MAPKs have roles in the circadian biological clock. MAPK signaling pathways and circadian clocks affect similar biological processes and defects in either pathway lead to many of the same types of human diseases (Goldsmith and Bell-Pedersen, 2013). Interestingly, our further analysis showed that MAPK and circadian pathway had similar dynamic changes after cocaine withdrawal (r > 0.70, p < 0.05), indicating that interaction between these two pathways mediated cocaine-induced neuroadaptations.

On another hand, a study compared changes of genes expression with cocaine abstinence between the mPFC and NAc (Freeman et al., 2010). The results showed that only a limited number of changes were observed in both brain regions. Another study showed that genes involved in cholinergic, glutamatergic, and GABAergic signaling transmissions in NAc were in response to chronic cocaine exposure and withdrawal (Eipper-Mains et al., 2013). In our study, only a few neurotransmitter receptors such as Grin2a and Grin2b, were changed in the PFC 7d after cocaine withdrawal. These inconsistent results may reflect distinct roles of genes in different brain regions. In addition, different cocaine treatment paradigms may contribute to these differences. For example, repeated cocaine treatment resulted in different changes of gene network with or without locomotor sensitization test (Eipper-Mains et al., 2013; Feng et al., 2014).

Repeated exposure to cocaine leads to sensitization, and may persist for weeks to years after cessation of drug taking, which may contribute to drug relapse (Hyman et al., 2006). Our finding indicates that maladaptive neural plasticity after repeated cocaine treatment is an ongoing, degenerative process with dynamic changes of the gene network at different stages of withdrawal. The long-lasting adverse effects on neuropsychiatric impairments that persist following abstinence lead to poor treatment outcomes. Understanding the dynamic profile of gene expression after abstinence could be helpful in developing new therapeutic approaches according to different periods of abstinence.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.